My title

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 _stdcall DSAUPD(int *ido, char* bmat, unsigned int *, 
00036           int *n, char *which, unsigned int *,
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 _stdcall DSEUPD(bool *rvec, char *howmny, unsigned int *,
00044           logical *select, double *d, double *z,
00045           int *ldz, double *sigma, char *bmat, unsigned int *,
00046           int  *n, char *which, unsigned int *,
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), 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 
00084 
00085 extern "C" int pfsfct(int neqns, double *diag, double **penv, int nblks, int *xblk,
00086         OFFDBLK **begblk, OFFDBLK *first, int *rowblks);
00087 
00088 extern "C" void pfsslv(int neqns, double *diag, double **penv, int nblks, 
00089          int *xblk, double *rhs, OFFDBLK **begblk);
00090 
00091 
00092 
00093 int
00094 SymArpackSolver::solve(void)
00095 { 
00096     if (theSOE == 0) {
00097  cerr << "WARNING SymArpackSolver::solve(void)- ";
00098  cerr << " No EigenSOE object has been set\n";
00099  return -1;
00100     }
00101 
00102     int      nblks = theSOE->nblks;
00103     int      *xblk = theSOE->xblk;
00104     //    int      *invp = theSOE->invp;
00105     double   *diag = theSOE->diag;
00106     double   **penv = theSOE->penv;
00107     int      *rowblks = theSOE->rowblks;
00108     OFFDBLK  **begblk = theSOE->begblk;
00109     OFFDBLK  *first = theSOE->first;
00110     
00111     int n = theSOE->size;
00112 
00113     // check for quick return
00114     if (n == 0)
00115  return 0;
00116 
00117  timer (FACTOR);
00118  if (factored == false) {
00119 
00120     //factor the matrix
00121     //call the "C" function to do the numerical factorization.
00122     int factor;
00123     factor = pfsfct(n, diag, penv, nblks, xblk, begblk, first, rowblks);
00124     if (factor>0) {
00125     cerr << "In SymArpackSolver: error in factorization.\n";
00126     return -1;
00127     }
00128     factored = true;
00129  }
00130 
00131 
00132     int nev = theNev;
00133     int ncv = getNCV(n, nev);
00134 
00135     // set up the space for ARPACK functions.
00136     int ldv = n;
00137     
00138     double *v = new double[ldv * ncv];
00139     double *workl = new double[ncv * (ncv+8)];
00140     double *workd = new double[3 * n];
00141     double *d = new double[ncv * 2];
00142     double *resid = new double[n];
00143  
00144 //    char *which = new char[2];
00145     char bmat = 'G';
00146 //    which = "LM";
00147 
00148     int *iparam = new int[11];
00149     int *iwork = new int[n];
00150 
00151     int lworkl = ncv*ncv + 8*ncv;
00152     int maxitr, mode, nconv;
00153 
00154     double tol = 0.0;
00155     int info = 0;
00156     maxitr = 1000;
00157     mode = 3;
00158 
00159     iparam[2] = maxitr;
00160     iparam[6] = mode; 
00161 
00162     bool rvec = true;
00163     char howmy = 'A';
00164     logical *select = new logical[ncv];
00165 
00166     iparam[0] = 1;
00167     int ido = 0;
00168     int *ipntr = new int[11];
00169     int ierr = 0;
00170     unsigned int sizeWhich =2;
00171     unsigned int sizeBmat =1;
00172     unsigned int sizeHowmany =1;
00173     while (1) { 
00174  
00175 #ifdef _WIN32
00176  DSAUPD(&ido, &bmat, &sizeBmat, &n, "LM", &sizeWhich, &nev, &tol, resid, 
00177         &ncv, v, &ldv,
00178         iparam, ipntr, workd, workl, &lworkl, &info);
00179 #else
00180  dsaupd_(&ido, &bmat, &n, "LM", &nev, &tol, resid, &ncv, v, &ldv,
00181   iparam, ipntr, workd, workl, &lworkl, &info);
00182 #endif
00183       if (ido == -1) {
00184 
00185    myMv(n, &workd[ipntr[0]-1], &workd[ipntr[1]-1]); 
00186 
00187    pfsslv(n, diag, penv, nblks, xblk, &workd[ipntr[1] - 1], begblk);
00188    continue;
00189       } else if (ido == 1) {
00190           double ratio = 1.0;
00191    myCopy(n, &workd[ipntr[2]-1], &workd[ipntr[1]-1]);
00192 
00193    pfsslv(n, diag, penv, nblks, xblk, &workd[ipntr[1] - 1], begblk);
00194 
00195    continue;
00196       } else if (ido == 2) {     
00197    myMv(n, &workd[ipntr[0]-1], &workd[ipntr[1]-1]);
00198 
00199    continue;
00200       }
00201 
00202       break;
00203     }
00204 
00205     if (info < 0) {
00206         cerr << "BandArpackSolver::Error with _saupd info = " << info <<endl;
00207  return info;
00208     } else {
00209         if (info == 1) {
00210      cerr << "BandArpackSolver::Maximum number of iteration reached." << endl;
00211  } else if (info == 3) {
00212      cerr << "BandArpackSolver::No Shifts could be applied during implicit," << endl;
00213      cerr << "Arnoldi update, try increasing NCV." << endl;
00214  }
00215  double sigma = theSOE->shift;
00216  if (iparam[4] > 0) {
00217 #ifdef _WIN32
00218      DSEUPD(&rvec, &howmy, &sizeHowmany, select, d, v, &ldv, &sigma, &bmat,
00219      &sizeBmat, &n, "LM", &sizeWhich,
00220      &nev, &tol, resid, &ncv, v, &ldv, iparam, ipntr, workd,
00221      workl, &lworkl, &info);
00222 #else
00223 
00224 
00225      dseupd_(&rvec, &howmy, select, d, v, &ldv, &sigma, &bmat, &n, "LM",
00226       &nev, &tol, resid, &ncv, v, &ldv, iparam, ipntr, workd,
00227       workl, &lworkl, &info);
00228 #endif
00229      if (info != 0) {
00230          cerr << "BandArpackSolver::Error with dseupd_" << info;
00231   return -1;
00232      }
00233  }
00234     }
00235 
00236     value = d;
00237     vector = v;
00238 
00239     theSOE->factored = true;
00240 
00241     delete [] workl;
00242     delete [] workd;
00243     delete [] resid;
00244     delete [] iparam;
00245     delete [] iwork;
00246 #ifdef _WIN32
00247     // cannot invoke the destructor on select .. causes seg in windows!!!
00248 #else
00249     delete [] select;
00250 #endif    
00251 
00252     delete [] ipntr;
00253 
00254     
00255     return 0;
00256 }
00257 
00258 
00259 int
00260 SymArpackSolver::setSize()
00261 {
00262     // nothing to do
00263     // if iPiv not big enough, free it and get one large enough
00264     int size = theSOE->size;    
00265     
00266     if (eigenV == 0 || eigenV->Size() != size) {
00267  if (eigenV != 0)
00268      delete eigenV;
00269 
00270  eigenV = new Vector(size);
00271  if (eigenV == 0 || eigenV->Size() != size) {
00272      cerr << "WARNING BandGenLinLapackSolver::setSize() ";
00273      cerr << " - ran out of memory for eigenVector of size ";
00274      cerr << theSOE->size << endl;
00275      return -2;     
00276  }
00277     }
00278      
00279     return 0;
00280 }
00281 
00282 int
00283 SymArpackSolver::setEigenSOE(SymArpackSOE &theEigenSOE)
00284 {
00285   theSOE = &theEigenSOE;
00286   return 0;
00287 }
00288 
00289 
00290 int 
00291 SymArpackSolver::getNCV(int n, int nev)
00292 {
00293     int result;
00294     if (2*nev > nev+8) {
00295         result = nev+8;
00296     } else {
00297         result = 2*nev;
00298     }
00299 
00300     if (result >= n) {
00301         result = n;
00302     }
00303 
00304     return result;
00305 }
00306 
00307 
00308 void
00309 SymArpackSolver::myMv(int n, double *v, double *result)
00310 {
00311     double *tempX = new double[n];
00312     int      *invp = theSOE->invp;    
00313     int i;
00314     
00315     for (i=0; i<n; i++) {
00316         tempX[i] = v[invp[i]];
00317     }
00318     for (i=0; i<n; i++) {
00319         v[i] = tempX[i];
00320     }
00321 
00322     Vector x(v, n);
00323     Vector y(result,n);
00324 
00325     y.Zero();
00326     AnalysisModel *theAnalysisModel = theSOE->theModel;
00327 
00328     // loop over the FE_Elements
00329     FE_Element *elePtr;
00330     FE_EleIter &theEles = theAnalysisModel->getFEs();    
00331     while((elePtr = theEles()) != 0) {
00332         elePtr->zeroResidual();
00333  elePtr->addM_Force(x,1.0);
00334  const Vector &a = elePtr->getResidual(0);
00335  y.Assemble(a,elePtr->getID(),1.0);
00336  // const Vector &a = elePtr->getM_Force(x,1.0);
00337  // y.Assemble(a,elePtr->getID(),1.0);
00338     }
00339 
00340     // loop over the DOF_Groups
00341     Integrator *theIntegrator = 0;
00342     DOF_Group *dofPtr;
00343     DOF_GrpIter &theDofs = theAnalysisModel->getDOFs();
00344     while ((dofPtr = theDofs()) != 0) {
00345         dofPtr->zeroUnbalance();
00346  dofPtr->addM_Force(x,1.0);
00347  const Vector &a = dofPtr->getUnbalance(theIntegrator);
00348  y.Assemble(a,dofPtr->getID(),1.0);
00349     }
00350 
00351     for (i=0; i<n; i++) {
00352         tempX[invp[i]] = result[i];
00353     }
00354     for (i=0; i<n; i++) {
00355         result[i] = tempX[i];
00356     }
00357     delete [] tempX;
00358 }
00359 
00360 void
00361 SymArpackSolver::myCopy(int n, double *v, double *result)
00362 {
00363     for (int i=0; i<n; i++) {
00364         result[i] = v[i];
00365     }
00366 }
00367 
00368 const Vector &
00369 SymArpackSolver::getEigenvector(int mode)
00370 {
00371     int      *invp = theSOE->invp;
00372     
00373     if (mode <= 0 || mode > theNev) {
00374         cerr << "BandArpackSOE::mode is out of range(1 - nev)";
00375  exit (0);
00376     }
00377     int size = theSOE->size;
00378 
00379     int index = (mode - 1) * size;
00380     for (int i=0; i<size; i++) {
00381  (*eigenV)[i] = vector[index + invp[i]];
00382     } 
00383 
00384     return *eigenV;      
00385 }
00386 
00387 
00388 double
00389 SymArpackSolver::getEigenvalue(int mode)
00390 {
00391     if (mode <= 0 || mode > theNev) {
00392         cerr << "BandArpackSOE::mode is out of range(1 - nev)";
00393  exit (0);
00394     }
00395     return value[mode - 1];
00396 }
00397 
00398 
00399 int
00400 SymArpackSolver::sendSelf(int cTAg, Channel &theChannel)
00401 {
00402     // nothing to do
00403     return 0;
00404 }
00405 
00406 int
00407 SymArpackSolver::recvSelf(int cTag,
00408         Channel &theChannel, FEM_ObjectBroker &theBroker)
00409 {
00410     // nothing to do
00411     return 0;
00412 }
00413 
00414 
00415
SymArpackSolver.cpp
