Main Page   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members  

SymArpackSOE.cpp

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