My title

00001 // File: ~/system_of_eqn/eigenSOE/BandArpackSOE.C
00002 //
00003 // Written: Jun Peng 
00004 // Created: Febuary 1999
00005 // Revision: A
00006 //
00007 // Description: This file contains the class definition for BandArpackSOE
00008 // BandArpackSOE is a subclass of EigenSOE. It uses the LAPACK storage
00009 // scheme to store the components of the K, M matrix, which is a full matrix.
00010 // It uses the ARPACK to do eigenvalue analysis.
00011 
00012 #include <BandArpackSOE.h>
00013 #include <BandArpackSolver.h>
00014 #include <Matrix.h>
00015 #include <Graph.h>
00016 #include <Vertex.h>
00017 #include <VertexIter.h>
00018 #include <f2c.h>
00019 #include <math.h>
00020 #include <Channel.h>
00021 #include <FEM_ObjectBroker.h>
00022 #include <AnalysisModel.h>
00023 
00024 BandArpackSOE::BandArpackSOE(BandArpackSolver &theSolvr,  
00025         AnalysisModel &aModel, double theShift)
00026 :EigenSOE(theSolvr, EigenSOE_TAGS_BandArpackSOE),
00027  size(0), numSubD(0), numSuperD(0), A(0), 
00028  Asize(0), factored(false), shift(theShift), theModel(&aModel)
00029 {
00030     theSolvr.setEigenSOE(*this);
00031 }
00032 
00033 
00034 int
00035 BandArpackSOE::getNumEqn(void) const
00036 {
00037     return size;
00038 }
00039     
00040 BandArpackSOE::~BandArpackSOE()
00041 {
00042     if (A != 0) delete [] A;
00043 }
00044 
00045 int 
00046 BandArpackSOE::setSize(Graph &theGraph)
00047 {
00048     int result = 0;
00049     size = theGraph.getNumVertex();
00050         
00051     // determine the number of superdiagonals and subdiagonals
00052     
00053     numSubD = 0;
00054     numSuperD = 0;
00055 
00056     Vertex *vertexPtr;
00057     VertexIter &theVertices = theGraph.getVertices();
00058     
00059     while ((vertexPtr = theVertices()) != 0) {
00060  int vertexNum = vertexPtr->getTag();
00061  const ID &theAdjacency = vertexPtr->getAdjacency();
00062  for (int i=0; i<theAdjacency.Size(); i++) {
00063      int otherNum = theAdjacency(i);
00064      int diff = vertexNum - otherNum;
00065      if (diff > 0) {
00066   if (diff > numSuperD)
00067       numSuperD = diff;
00068      } else 
00069   if (diff < numSubD)
00070       numSubD = diff;
00071  }
00072     }
00073     numSubD *= -1;
00074 
00075     int newSize = size * (2*numSubD + numSuperD +1);
00076     if (newSize > Asize) { // we have to get another space for A
00077 
00078  if (A != 0) 
00079      delete [] A;
00080 
00081  A = new double[newSize];
00082  
00083         if (A == 0) {
00084             cerr << "WARNING BandGenLinSOE::BandGenLinSOE :";
00085      cerr << " ran out of memory for A (size,super,sub) (";
00086      cerr << size <<", " << numSuperD << ", " << numSubD << ") \n";
00087      Asize = 0; size = 0; numSubD = 0; numSuperD = 0;
00088      result= -1;
00089         }
00090  else  
00091      Asize = newSize;
00092     }
00093 
00094     // zero the matrix
00095     for (int i=0; i<Asize; i++)
00096  A[i] = 0;
00097  
00098     factored = false;
00099 
00100     // invoke setSize() on the Solver
00101     EigenSolver *theSolvr = this->getSolver();
00102     int solverOK = theSolvr->setSize();
00103     if (solverOK < 0) {
00104  cerr << "WARNING:BandArpackSOE::setSize :";
00105  cerr << " solver failed setSize()\n";
00106  return solverOK;
00107     } 
00108    
00109     return result;    
00110 }
00111 
00112 int 
00113 BandArpackSOE::addA(const Matrix &m, const ID &id, double fact)
00114 {
00115     // check for a quick return 
00116     if (fact == 0.0)  return 0;
00117 
00118     // check that m and id are of similar size
00119     int idSize = id.Size();    
00120     if (idSize != m.noRows() && idSize != m.noCols()) {
00121  cerr << "BandArpackSOE::addA() - Matrix and ID not of similar sizes\n";
00122  return -1;
00123     }
00124     
00125     int ldA = 2*numSubD + numSuperD + 1;
00126 
00127     if (fact == 1.0) { // do not need to multiply 
00128  for (int i=0; i<idSize; i++) {
00129      int col = id(i);
00130      if (col < size && col >= 0) {
00131   double *coliiPtr = A + col*ldA + numSubD + numSuperD;
00132   for (int j=0; j<idSize; j++) {
00133       int row = id(j);
00134       if (row <size && row >= 0) {      
00135    int diff = col - row;
00136    if (diff > 0) {
00137        if (diff <= numSuperD) {
00138     double *APtr = coliiPtr - diff;
00139     *APtr += m(j,i);
00140        }   
00141 
00142    } else {
00143        diff *= -1;
00144        if (diff <= numSubD) {
00145     double *APtr = coliiPtr + diff;
00146     *APtr += m(j,i);
00147        }
00148    }
00149       }
00150   }  // for j
00151      } 
00152  }  // for i
00153     } else {
00154  for (int i=0; i<idSize; i++) {
00155      int col = id(i);
00156      if (col < size && col >= 0) {
00157   double *coliiPtr = A + col*ldA + numSubD + numSuperD;
00158   for (int j=0; j<idSize; j++) {
00159       int row = id(j);
00160       if (row <size && row >= 0) {      
00161    int diff = col - row;
00162    if (diff > 0) {
00163        if (diff <= numSuperD) {
00164     double *APtr = coliiPtr - diff;
00165     *APtr += m(j,i) *fact;
00166        }
00167    } else {
00168        diff *= -1;
00169        if (diff <= numSubD) {
00170     double *APtr = coliiPtr + diff;
00171     *APtr += m(j,i) *fact;
00172        }
00173    }
00174       }
00175   }  // for j
00176      } 
00177  }  // for i
00178     }    
00179     return 0;
00180 }
00181 
00182 
00183 void 
00184 BandArpackSOE::zeroA(void)
00185 {
00186     double *Aptr = A;
00187     int theSize = size*(2*numSubD+numSuperD+1);
00188     for (int i=0; i<theSize; i++)
00189  *Aptr++ = 0;
00190     
00191     factored = false;
00192 }
00193 
00194 int 
00195 BandArpackSOE::addM(const Matrix &m, const ID &id, double fact)
00196 {
00197   return this->addA(m, id, -shift);
00198 }   
00199  
00200 void 
00201 BandArpackSOE::zeroM(void)
00202 {
00203   
00204 }
00205 
00206 
00207 double 
00208 BandArpackSOE::getShift(void)
00209 {
00210     return shift;
00211 }
00212 
00213 
00214 int 
00215 BandArpackSOE::sendSelf(int commitTag, Channel &theChannel)
00216 {
00217     return 0;
00218 }
00219     
00220 int 
00221 BandArpackSOE::recvSelf(int commitTag, Channel &theChannel, 
00222    FEM_ObjectBroker &theBroker)
00223 {
00224     return 0;
00225 }
00226
BandArpackSOE.cpp
