SymArpackSolver.cppGo to the documentation of this file.00001 // File: ~/system_of_eqn/linearSOE/LawSolver/SymArpackSolver.C 00002 // 00003 // Written: Jun Peng 00004 // Created: 2/1999 00005 // Revision: A 00006 // 00007 // Description: This file contains the class definition for 00008 // SymArpackSolver. It solves the SymArpackSOE object by calling 00009 // some "C" functions. The solver used here is generalized sparse 00010 // solver. The user can choose three different ordering schema. 00011 // 00012 // What: "@(#) SymArpackSolver.C, revA" 00013 00014 00015 #include "SymArpackSOE.h" 00016 #include "SymArpackSolver.h" 00017 #include <f2c.h> 00018 #include <math.h> 00019 #include <Channel.h> 00020 #include <FEM_ObjectBroker.h> 00021 #include <AnalysisModel.h> 00022 #include <DOF_GrpIter.h> 00023 #include <DOF_Group.h> 00024 #include <FE_EleIter.h> 00025 #include <FE_Element.h> 00026 #include <Integrator.h> 00027 00028 extern "C" { 00029 #include <FeStructs.h> 00030 #include <globalVars.h> 00031 #include <tim.h> 00032 } 00033 00034 #ifdef _WIN32 00035 extern "C" int DSAUPD(int *ido, char* bmat, 00036 int *n, char *which, 00037 int *nev, 00038 double *tol, double *resid, int *ncv, double *v, 00039 int *ldv, 00040 int *iparam, int *ipntr, double *workd, double *workl, 00041 int *lworkl, int *info); 00042 00043 extern "C" int DSEUPD(bool *rvec, char *howmny, 00044 logical *select, double *d, double *z, 00045 int *ldz, double *sigma, char *bmat, 00046 int *n, char *which, 00047 int *nev, double *tol, double *resid, int *ncv, 00048 double *v, 00049 int *ldv, int *iparam, int *ipntr, double *workd, 00050 double *workl, int *lworkl, int *info); 00051 00052 #else 00053 00054 extern "C" int dsaupd_(int *ido, char* bmat, int *n, char *which, int *nev, 00055 double *tol, double *resid, int *ncv, double *v, int *ldv, 00056 int *iparam, int *ipntr, double *workd, double *workl, 00057 int *lworkl, int *info); 00058 00059 extern "C" int dseupd_(bool *rvec, char *howmny, logical *select, double *d, double *z, 00060 int *ldz, double *sigma, char *bmat, int *n, char *which, 00061 int *nev, double *tol, double *resid, int *ncv, double *v, 00062 int *ldv, int *iparam, int *ipntr, double *workd, 00063 double *workl, int *lworkl, int *info); 00064 00065 #endif 00066 00067 SymArpackSolver::SymArpackSolver(int numE) 00068 :EigenSolver(EigenSOLVER_TAGS_SymArpackSolver), 00069 theSOE(0), factored(false), theNev(numE), value(0), eigenV(0) 00070 { 00071 // nothing to do. 00072 } 00073 00074 00075 SymArpackSolver::~SymArpackSolver() 00076 { 00077 if (value != 0) 00078 delete [] value; 00079 00080 if (vector !=0) 00081 delete [] vector; 00082 00083 if (eigenV !=0) 00084 delete eigenV; 00085 } 00086 00087 00088 extern "C" int pfsfct(int neqns, double *diag, double **penv, int nblks, int *xblk, 00089 OFFDBLK **begblk, OFFDBLK *first, int *rowblks); 00090 00091 extern "C" void pfsslv(int neqns, double *diag, double **penv, int nblks, 00092 int *xblk, double *rhs, OFFDBLK **begblk); 00093 00094 00095 00096 int 00097 SymArpackSolver::solve(void) 00098 { 00099 if (theSOE == 0) { 00100 opserr << "WARNING SymArpackSolver::solve(void)- "; 00101 opserr << " No EigenSOE object has been set\n"; 00102 return -1; 00103 } 00104 00105 int nblks = theSOE->nblks; 00106 int *xblk = theSOE->xblk; 00107 // int *invp = theSOE->invp; 00108 double *diag = theSOE->diag; 00109 double **penv = theSOE->penv; 00110 int *rowblks = theSOE->rowblks; 00111 OFFDBLK **begblk = theSOE->begblk; 00112 OFFDBLK *first = theSOE->first; 00113 00114 int n = theSOE->size; 00115 00116 // check for quick return 00117 if (n == 0) 00118 return 0; 00119 00120 // timer (FACTOR); 00121 if (factored == false) { 00122 00123 //factor the matrix 00124 //call the "C" function to do the numerical factorization. 00125 int factor; 00126 factor = pfsfct(n, diag, penv, nblks, xblk, begblk, first, rowblks); 00127 if (factor>0) { 00128 opserr << "In SymArpackSolver: error in factorization.\n"; 00129 return -1; 00130 } 00131 factored = true; 00132 } 00133 00134 00135 int nev = theNev; 00136 int ncv = getNCV(n, nev); 00137 00138 // set up the space for ARPACK functions. 00139 int ldv = n; 00140 00141 double *v = new double[ldv * ncv]; 00142 double *workl = new double[ncv * (ncv+8)]; 00143 double *workd = new double[3 * n]; 00144 double *d = new double[ncv * 2]; 00145 double *resid = new double[n]; 00146 00147 // char *which = new char[2]; 00148 char bmat = 'G'; 00149 // which = "LM"; 00150 00151 int *iparam = new int[11]; 00152 int *iwork = new int[n]; 00153 00154 int lworkl = ncv*ncv + 8*ncv; 00155 int maxitr, mode, nconv; 00156 00157 double tol = 0.0; 00158 int info = 0; 00159 maxitr = 1000; 00160 mode = 3; 00161 00162 iparam[2] = maxitr; 00163 iparam[6] = mode; 00164 00165 bool rvec = true; 00166 char howmy = 'A'; 00167 logical *select = new logical[ncv]; 00168 00169 iparam[0] = 1; 00170 int ido = 0; 00171 int *ipntr = new int[11]; 00172 int ierr = 0; 00173 unsigned int sizeWhich =2; 00174 unsigned int sizeBmat =1; 00175 unsigned int sizeHowmany =1; 00176 while (1) { 00177 00178 #ifdef _WIN32 00179 DSAUPD(&ido, &bmat, &n, "LM", &nev, &tol, resid, &ncv, v, &ldv, 00180 iparam, ipntr, workd, workl, &lworkl, &info); 00181 #else 00182 dsaupd_(&ido, &bmat, &n, "LM", &nev, &tol, resid, &ncv, v, &ldv, 00183 iparam, ipntr, workd, workl, &lworkl, &info); 00184 #endif 00185 if (ido == -1) { 00186 00187 myMv(n, &workd[ipntr[0]-1], &workd[ipntr[1]-1]); 00188 00189 pfsslv(n, diag, penv, nblks, xblk, &workd[ipntr[1] - 1], begblk); 00190 continue; 00191 } else if (ido == 1) { 00192 double ratio = 1.0; 00193 myCopy(n, &workd[ipntr[2]-1], &workd[ipntr[1]-1]); 00194 00195 pfsslv(n, diag, penv, nblks, xblk, &workd[ipntr[1] - 1], begblk); 00196 00197 continue; 00198 } else if (ido == 2) { 00199 myMv(n, &workd[ipntr[0]-1], &workd[ipntr[1]-1]); 00200 00201 continue; 00202 } 00203 00204 break; 00205 } 00206 00207 if (info < 0) { 00208 opserr << "BandArpackSolver::Error with _saupd info = " << info <<endln; 00209 return info; 00210 } else { 00211 if (info == 1) { 00212 opserr << "BandArpackSolver::Maximum number of iteration reached." << endln; 00213 } else if (info == 3) { 00214 opserr << "BandArpackSolver::No Shifts could be applied during implicit," << endln; 00215 opserr << "Arnoldi update, try increasing NCV." << endln; 00216 } 00217 double sigma = theSOE->shift; 00218 if (iparam[4] > 0) { 00219 #ifdef _WIN32 00220 DSEUPD(&rvec, &howmy, select, d, v, &ldv, &sigma, &bmat, &n, "LM", 00221 &nev, &tol, resid, &ncv, v, &ldv, iparam, ipntr, workd, 00222 workl, &lworkl, &info); 00223 #else 00224 00225 00226 dseupd_(&rvec, &howmy, select, d, v, &ldv, &sigma, &bmat, &n, "LM", 00227 &nev, &tol, resid, &ncv, v, &ldv, iparam, ipntr, workd, 00228 workl, &lworkl, &info); 00229 #endif 00230 if (info != 0) { 00231 opserr << "BandArpackSolver::Error with dseupd_" << info; 00232 return -1; 00233 } 00234 } 00235 } 00236 00237 value = d; 00238 vector = v; 00239 00240 theSOE->factored = true; 00241 00242 delete [] workl; 00243 delete [] workd; 00244 delete [] resid; 00245 delete [] iparam; 00246 delete [] iwork; 00247 #ifdef _WIN32 00248 // cannot invoke the destructor on select .. causes seg in windows!!! 00249 #else 00250 delete [] select; 00251 #endif 00252 00253 delete [] ipntr; 00254 00255 00256 return 0; 00257 } 00258 00259 00260 int 00261 SymArpackSolver::setSize() 00262 { 00263 // nothing to do 00264 // if iPiv not big enough, free it and get one large enough 00265 int size = theSOE->size; 00266 00267 if (eigenV == 0 || eigenV->Size() != size) { 00268 if (eigenV != 0) 00269 delete eigenV; 00270 00271 eigenV = new Vector(size); 00272 if (eigenV == 0 || eigenV->Size() != size) { 00273 opserr << "WARNING BandGenLinLapackSolver::setSize() "; 00274 opserr << " - ran out of memory for eigenVector of size "; 00275 opserr << theSOE->size << endln; 00276 return -2; 00277 } 00278 } 00279 00280 return 0; 00281 } 00282 00283 int 00284 SymArpackSolver::setEigenSOE(SymArpackSOE &theEigenSOE) 00285 { 00286 theSOE = &theEigenSOE; 00287 return 0; 00288 } 00289 00290 00291 int 00292 SymArpackSolver::getNCV(int n, int nev) 00293 { 00294 int result; 00295 if (2*nev > nev+8) { 00296 result = nev+8; 00297 } else { 00298 result = 2*nev; 00299 } 00300 00301 if (result >= n) { 00302 result = n; 00303 } 00304 00305 return result; 00306 } 00307 00308 00309 void 00310 SymArpackSolver::myMv(int n, double *v, double *result) 00311 { 00312 double *tempX = new double[n]; 00313 int *invp = theSOE->invp; 00314 int i; 00315 00316 for (i=0; i<n; i++) { 00317 tempX[i] = v[invp[i]]; 00318 } 00319 for (i=0; i<n; i++) { 00320 v[i] = tempX[i]; 00321 } 00322 00323 Vector x(v, n); 00324 Vector y(result,n); 00325 00326 y.Zero(); 00327 AnalysisModel *theAnalysisModel = theSOE->theModel; 00328 00329 // loop over the FE_Elements 00330 FE_Element *elePtr; 00331 FE_EleIter &theEles = theAnalysisModel->getFEs(); 00332 while((elePtr = theEles()) != 0) { 00333 const Vector &a = elePtr->getM_Force(x,1.0); 00334 y.Assemble(a,elePtr->getID(),1.0); 00335 } 00336 00337 // loop over the DOF_Groups 00338 Integrator *theIntegrator = 0; 00339 DOF_Group *dofPtr; 00340 DOF_GrpIter &theDofs = theAnalysisModel->getDOFs(); 00341 while ((dofPtr = theDofs()) != 0) { 00342 const Vector &a = dofPtr->getM_Force(x,1.0); 00343 y.Assemble(a,dofPtr->getID(),1.0); 00344 } 00345 00346 for (i=0; i<n; i++) { 00347 tempX[invp[i]] = result[i]; 00348 } 00349 for (i=0; i<n; i++) { 00350 result[i] = tempX[i]; 00351 } 00352 delete [] tempX; 00353 } 00354 00355 void 00356 SymArpackSolver::myCopy(int n, double *v, double *result) 00357 { 00358 for (int i=0; i<n; i++) { 00359 result[i] = v[i]; 00360 } 00361 } 00362 00363 const Vector & 00364 SymArpackSolver::getEigenvector(int mode) 00365 { 00366 int *invp = theSOE->invp; 00367 00368 if (mode <= 0 || mode > theNev) { 00369 opserr << "BandArpackSOE::mode is out of range(1 - nev)"; 00370 exit (0); 00371 } 00372 int size = theSOE->size; 00373 00374 int index = (mode - 1) * size; 00375 for (int i=0; i<size; i++) { 00376 (*eigenV)[i] = vector[index + invp[i]]; 00377 } 00378 00379 return *eigenV; 00380 } 00381 00382 00383 double 00384 SymArpackSolver::getEigenvalue(int mode) 00385 { 00386 if (mode <= 0 || mode > theNev) { 00387 opserr << "BandArpackSOE::mode is out of range(1 - nev)"; 00388 exit (0); 00389 } 00390 return value[mode - 1]; 00391 } 00392 00393 00394 int 00395 SymArpackSolver::sendSelf(int cTAg, Channel &theChannel) 00396 { 00397 // nothing to do 00398 return 0; 00399 } 00400 00401 int 00402 SymArpackSolver::recvSelf(int cTag, 00403 Channel &theChannel, FEM_ObjectBroker &theBroker) 00404 { 00405 // nothing to do 00406 return 0; 00407 } 00408 00409 00410
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