SymArpackSOE.cppGo to the documentation of this file.00001 // Written: Jun Peng 00002 // Created: 12/98 00003 // Revision: A 00004 // 00005 // Description: This file contains the class definition for 00006 // LawSparseinSolver. It solves the SymArpackSOEobject by calling 00007 // some "C" functions. The solver used here is generalized sparse 00008 // solver. The user can choose three different ordering schema. 00009 // 00010 // What: "@(#) SymArpackSOE.C, revA" 00011 00012 #include "SymArpackSOE.h" 00013 #include "SymArpackSolver.h" 00014 #include <Matrix.h> 00015 #include <Graph.h> 00016 #include <Vertex.h> 00017 #include <VertexIter.h> 00018 #include <math.h> 00019 #include <Channel.h> 00020 #include <FEM_ObjectBroker.h> 00021 #include <AnalysisModel.h> 00022 #include <Vector.h> 00023 00024 SymArpackSOE::SymArpackSOE(SymArpackSolver &the_Solver, AnalysisModel &aModel, 00025 double theShift) 00026 :EigenSOE(the_Solver, EigenSOE_TAGS_SymArpackSOE), 00027 size(0), nnz(0), colA(0), rowStartA(0), 00028 factored(false), shift(theShift), theModel(&aModel), 00029 nblks(0), xblk(0), invp(0), diag(0), penv(0), rowblks(0), 00030 begblk(0), first(0) 00031 { 00032 the_Solver.setEigenSOE(*this); 00033 } 00034 00035 00036 SymArpackSOE::~SymArpackSOE() 00037 { 00038 /* 00039 if (invp != 0) free((void *)invp); 00040 if (xblk != 0) free((void *)xblk); 00041 if (penv != 0) { 00042 free((void *)penv[0]); 00043 free((void *)penv); 00044 } 00045 if (diag != 0) free((void *)diag); 00046 if (rowblks != 0) free((void *)rowblks); 00047 if (begblk != 0) free((void *)begblk); 00048 if (first != 0) free((void *)first); 00049 00050 if (rowStartA != 0) delete [] rowStartA; 00051 if (colA != 0) delete [] colA; 00052 */ 00053 } 00054 00055 00056 int 00057 SymArpackSOE::getNumEqn(void) const 00058 { 00059 return size; 00060 } 00061 00062 // extern "C" int symFactorization(int *fxadj, int *adjncy, int neq, int LSPARSE); 00063 extern "C" int symFactorization(int *fxadj, int *adjncy, int neq, int LSPARSE, 00064 int **xblkMY, int **invpMY, int **rowblksMY, 00065 OFFDBLK ***begblkMY, OFFDBLK **firstMY, 00066 double ***penvMY, double **diagMY); 00067 00068 int 00069 SymArpackSOE::setSize(Graph &theGraph) 00070 { 00071 int result = 0; 00072 int oldSize = size; 00073 size = theGraph.getNumVertex(); 00074 00075 // fist itearte through the vertices of the graph to get nnz 00076 Vertex *theVertex; 00077 int newNNZ = 0; 00078 VertexIter &theVertices = theGraph.getVertices(); 00079 while ((theVertex = theVertices()) != 0) { 00080 const ID &theAdjacency = theVertex->getAdjacency(); 00081 newNNZ += theAdjacency.Size(); 00082 } 00083 nnz = newNNZ; 00084 00085 if (colA != 0) 00086 delete [] colA; 00087 00088 colA = new int[newNNZ]; 00089 if (colA == 0) { 00090 opserr << "WARNING SymArpackSOE::SymArpackSOE :"; 00091 opserr << " ran out of memory for colA with nnz = "; 00092 opserr << newNNZ << " \n"; 00093 size = 0; nnz = 0; 00094 result = -1; 00095 } 00096 00097 factored = false; 00098 00099 if (rowStartA != 0) 00100 delete [] rowStartA; 00101 rowStartA = new int[size+1]; 00102 if (rowStartA == 0) { 00103 opserr << "SymArpackSOE::ran out of memory for rowStartA." << endln; 00104 result = -1; 00105 } 00106 00107 // fill in rowStartA and colA 00108 if (size != 0) { 00109 rowStartA[0] = 0; 00110 int startLoc = 0; 00111 int lastLoc = 0; 00112 00113 for (int a=0; a<size; a++) { 00114 theVertex = theGraph.getVertexPtr(a); 00115 if (theVertex == 0) { 00116 opserr << "WARNING:SymArpackSOE::setSize :"; 00117 opserr << " vertex " << a << " not in graph! - size set to 0\n"; 00118 size = 0; 00119 return -1; 00120 } 00121 00122 const ID &theAdjacency = theVertex->getAdjacency(); 00123 int idSize = theAdjacency.Size(); 00124 00125 // now we have to place the entries in the ID into order in colA 00126 for (int i=0; i<idSize; i++) { 00127 int row = theAdjacency(i); 00128 bool foundPlace = false; 00129 00130 for (int j=startLoc; j<lastLoc; j++) 00131 if (colA[j] > row) { 00132 // move the entries already there one further on 00133 // and place col in current location 00134 for (int k=lastLoc; k>j; k--) 00135 colA[k] = colA[k-1]; 00136 colA[j] = row; 00137 foundPlace = true; 00138 j = lastLoc; 00139 } 00140 00141 if (foundPlace == false) // put in at the end 00142 colA[lastLoc] = row; 00143 00144 lastLoc++; 00145 } 00146 rowStartA[a+1] = lastLoc;; 00147 startLoc = lastLoc; 00148 } 00149 } 00150 00151 // begin to choose different ordering schema. 00152 // cout << "Enter DOF Numberer Type: \n"; 00153 // cout << "[1] Minimum Degree, [2] Nested Dissection, [3] RCM: "; 00154 int LSPARSE = 1; 00155 // cin >> LSPARSE; 00156 00157 // call "C" function to form elimination tree and do the symbolic factorization. 00158 // nblks = symFactorization(rowStartA, colA, size, LSPARSE); 00159 nblks = symFactorization(rowStartA, colA, size, LSPARSE, 00160 &xblk, &invp, &rowblks, &begblk, &first, &penv, &diag); 00161 00162 // invoke setSize() on the Solver 00163 EigenSolver *theSolvr = this->getSolver(); 00164 int solverOK = theSolvr->setSize(); 00165 if (solverOK < 0) { 00166 opserr << "WARNING:BandArpackSOE::setSize :"; 00167 opserr << " solver failed setSize()\n"; 00168 return solverOK; 00169 } 00170 00171 00172 return result; 00173 } 00174 00175 00176 int 00177 SymArpackSOE::addA(const Matrix &m, const ID &id, double fact) 00178 { 00179 // check for a quick return 00180 if (fact == 0.0) 00181 return 0; 00182 00183 int idSize = id.Size(); 00184 00185 // check that m and id are of similar size 00186 if (idSize != m.noRows() && idSize != m.noCols()) { 00187 opserr << "SymArpackSOE::addA() "; 00188 opserr << " - Matrix and ID not of similar sizes\n"; 00189 return -1; 00190 } 00191 00192 int *newID = new int [idSize]; 00193 int *isort = new int[idSize]; 00194 if (newID == 0 || isort ==0) { 00195 opserr << "WARNING SymArpackSOE::addA() :"; 00196 opserr << " ran out of memory for vectors (newID, isort)"; 00197 return -1; 00198 } 00199 00200 int i; 00201 for (i=0; i<idSize; i++) { 00202 newID[i] = id(i); 00203 if (newID[i] >= 0) 00204 newID[i] = invp[newID[i]]; 00205 } 00206 00207 long int i_eq, j_eq, iadd; 00208 int j, nee, lnee; 00209 int k, ipos, jpos; 00210 int it, jt; 00211 int iblk; 00212 int jblk; 00213 OFFDBLK *ptr; 00214 OFFDBLK *saveblk; 00215 double *fpt, *iloc, *loc; 00216 00217 nee = idSize; 00218 lnee = nee; 00219 00220 /* initialize isort */ 00221 for(i = 0, k = 0; i < lnee ; i++ ) 00222 { 00223 if( newID[i] >= 0 ) { 00224 isort[k] = i; 00225 k++; 00226 } 00227 } 00228 00229 lnee = k; 00230 00231 i = k - 1; 00232 do 00233 { 00234 k = 0 ; 00235 for (j = 0 ; j < i ; j++) 00236 { 00237 if ( newID[isort[j]] > newID[isort[j+1]]) { 00238 int temp = isort[j]; 00239 isort[j] = isort[j+1]; 00240 isort[j+1] = temp; 00241 00242 /* isort[j] ^= isort[j+1] ; 00243 isort[j+1] ^= isort[j] ; 00244 isort[j] ^= isort[j+1] ; 00245 */ 00246 k = j ; 00247 } 00248 } 00249 i = k ; 00250 } while ( k > 0) ; 00251 00252 i = 0 ; 00253 ipos = isort[i] ; 00254 k = rowblks[newID[ipos]] ; 00255 saveblk = begblk[k] ; 00256 00257 for (i=0; i<lnee; i++) 00258 { 00259 ipos = isort[i] ; 00260 i_eq = newID[ipos] ; 00261 iblk = rowblks[i_eq] ; 00262 iloc = penv[i_eq +1] - i_eq ; 00263 if (k < iblk) 00264 while (saveblk->row != i_eq) 00265 saveblk = saveblk->bnext ; 00266 00267 ptr = saveblk ; 00268 for (j=0; j< i ; j++) 00269 { 00270 jpos = isort[j] ; 00271 j_eq = newID[jpos] ; 00272 00273 if (ipos > jpos) { 00274 jt = ipos; 00275 it = jpos; 00276 } else { 00277 it = ipos; 00278 jt = jpos; 00279 } 00280 00281 if (j_eq >= xblk[iblk]) /* diagonal block */ 00282 { loc = iloc + j_eq ; 00283 *loc += m(it, jt) * fact; 00284 } else /* row segment */ 00285 { while((j_eq >= (ptr->next)->beg) && ((ptr->next)->row == i_eq)) 00286 ptr = ptr->next ; 00287 fpt = ptr->nz ; 00288 fpt[j_eq - ptr->beg] += m(it,jt) * fact; 00289 } 00290 } 00291 diag[i_eq] += m(ipos, ipos) * fact; 00292 } 00293 00294 delete [] newID; 00295 delete [] isort; 00296 00297 return 0; 00298 } 00299 00300 00301 int 00302 SymArpackSOE::addM(const Matrix &m, const ID &id, double fact) 00303 { 00304 return this->addA(m, id, -shift); 00305 } 00306 00307 00308 double 00309 SymArpackSOE::getShift(void) 00310 { 00311 return shift; 00312 } 00313 00314 void 00315 SymArpackSOE::zeroA(void) 00316 { 00317 factored = false; 00318 } 00319 00320 void 00321 SymArpackSOE::zeroM(void) 00322 { 00323 00324 } 00325 00326 int 00327 SymArpackSOE::sendSelf(int cTag, Channel &theChannel) 00328 { 00329 return 0; 00330 } 00331 00332 int 00333 SymArpackSOE::recvSelf(int cTag, Channel &theChannel, 00334 FEM_ObjectBroker &theBroker) 00335 { 00336 return 0; 00337 } 00338 00339 00340 00341 00342 00343 00344 00345 00346 00347 00348 00349
Generated on Mon Oct 23 15:05:27 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents