BandArpackSOE.cpp

Go to the documentation of this file.
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             opserr << "WARNING BandGenLinSOE::BandGenLinSOE :";
00085             opserr << " ran out of memory for A (size,super,sub) (";
00086             opserr << 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         opserr << "WARNING:BandArpackSOE::setSize :";
00105         opserr << " 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         opserr << "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 

Generated on Mon Oct 23 15:05:27 2006 for OpenSees by doxygen 1.5.0