SparseGenColLinSOE.cpp

Go to the documentation of this file.
00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.7 $
00022 // $Date: 2006/10/02 20:23:22 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/sparseGEN/SparseGenColLinSOE.cpp,v $
00024                                                                         
00025                                                                         
00026 // Written: fmk 
00027 // Created: 04/98
00028 // Revision: A
00029 //
00030 // Description: This file contains the implementation for SparseGenColLinSOE
00031 
00032 #include <SparseGenColLinSOE.h>
00033 #include <SparseGenColLinSolver.h>
00034 #include <Matrix.h>
00035 #include <Graph.h>
00036 #include <Vertex.h>
00037 #include <VertexIter.h>
00038 #include <math.h>
00039 
00040 #include <Channel.h>
00041 #include <FEM_ObjectBroker.h>
00042 
00043 SparseGenColLinSOE::SparseGenColLinSOE(SparseGenColLinSolver &the_Solver)
00044 :LinearSOE(the_Solver, LinSOE_TAGS_SparseGenColLinSOE),
00045  size(0), nnz(0), A(0), B(0), X(0), rowA(0), colStartA(0),
00046  vectX(0), vectB(0),
00047  Asize(0), Bsize(0),
00048  factored(false)
00049 {
00050     the_Solver.setLinearSOE(*this);
00051 }
00052 
00053 
00054 SparseGenColLinSOE::SparseGenColLinSOE(SparseGenColLinSolver &the_Solver, int classTag)
00055   :LinearSOE(the_Solver, classTag),
00056    size(0), nnz(0), A(0), B(0), X(0), rowA(0), colStartA(0),
00057    vectX(0), vectB(0),
00058    Asize(0), Bsize(0),
00059    factored(false)
00060 {
00061   //    the_Solver.setLinearSOE(*this);
00062 }
00063 
00064 
00065 SparseGenColLinSOE::SparseGenColLinSOE(int N, int NNZ, int *ColStartA, int *RowA,
00066                                  SparseGenColLinSolver &the_Solver)
00067 :LinearSOE(the_Solver, LinSOE_TAGS_SparseGenColLinSOE),
00068  size(N), nnz(NNZ), A(0), B(0), X(0), 
00069  rowA(RowA), colStartA(ColStartA), 
00070  vectX(0), vectB(0),
00071  Asize(0), Bsize(0),
00072  factored(false)
00073 {
00074 
00075     A = new double[NNZ];
00076         
00077     if (A == 0) {
00078         opserr << "WARNING :SparseGenColLinSOE::SparseGenColLinSOE :";
00079         opserr << " ran out of memory for A (NNZ) (";
00080         opserr << NNZ << ") \n";
00081         size = 0; nnz = 0;
00082     } else {
00083         // zero the matrix
00084         Asize = NNZ;
00085         for (int i=0; i<Asize; i++)
00086             A[i] = 0;
00087     
00088         B = new double[size];
00089         X = new double[size];
00090         
00091         if (B == 0 || X == 0) {
00092             opserr << "WARNING :SparseGenColLinSOE::SparseGenColLinSOE :";
00093             opserr << " ran out of memory for vectors (size) (";
00094             opserr << size << ") \n";
00095             size = 0; Bsize = 0; NNZ = 0;
00096         } else {
00097             Bsize = size;
00098             // zero the vectors
00099             for (int j=0; j<size; j++) {
00100                 B[j] = 0;
00101                 X[j] = 0;
00102             }
00103         }
00104     }
00105 
00106     vectX = new Vector(X,size);
00107     vectB = new Vector(B,size);    
00108     
00109     the_Solver.setLinearSOE(*this);
00110     
00111     // invoke setSize() on the Solver        
00112     if (the_Solver.setSize() < 0) {
00113         opserr << "WARNING :SparseGenColLinSOE::SparseGenColLinSOE :";
00114         opserr << " solver failed setSize() in constructor\n";
00115     }    
00116 }
00117 
00118     
00119 SparseGenColLinSOE::~SparseGenColLinSOE()
00120 {
00121     if (A != 0) delete [] A;
00122     if (B != 0) delete [] B;
00123     if (X != 0) delete [] X;
00124     if (colStartA != 0) delete [] colStartA;
00125     if (rowA != 0) delete []rowA;
00126     if (vectX != 0) delete vectX;    
00127     if (vectB != 0) delete vectB;        
00128 }
00129 
00130 
00131 int
00132 SparseGenColLinSOE::getNumEqn(void) const
00133 {
00134     return size;
00135 }
00136 
00137 int 
00138 SparseGenColLinSOE::setSize(Graph &theGraph)
00139 {
00140 
00141     int result = 0;
00142     int oldSize = size;
00143     size = theGraph.getNumVertex();
00144 
00145     // fist itearte through the vertices of the graph to get nnz
00146     Vertex *theVertex;
00147     int newNNZ = 0;
00148     VertexIter &theVertices = theGraph.getVertices();
00149     while ((theVertex = theVertices()) != 0) {
00150         const ID &theAdjacency = theVertex->getAdjacency();
00151         newNNZ += theAdjacency.Size() +1; // the +1 is for the diag entry
00152     }
00153     nnz = newNNZ;
00154 
00155     if (newNNZ > Asize) { // we have to get more space for A and rowA
00156         if (A != 0) 
00157             delete [] A;
00158         if (rowA != 0)
00159             delete [] rowA;
00160         
00161         A = new double[newNNZ];
00162         rowA = new int[newNNZ];
00163         
00164         if (A == 0 || rowA == 0) {
00165             opserr << "WARNING SparseGenColLinSOE::SparseGenColLinSOE :";
00166             opserr << " ran out of memory for A and rowA with nnz = ";
00167             opserr << newNNZ << " \n";
00168             size = 0; Asize = 0; nnz = 0;
00169             result =  -1;
00170         } 
00171         
00172         Asize = newNNZ;
00173     }
00174 
00175     // zero the matrix
00176     for (int i=0; i<Asize; i++)
00177         A[i] = 0;
00178         
00179     factored = false;
00180     
00181     if (size > Bsize) { // we have to get space for the vectors
00182         
00183         // delete the old       
00184         if (B != 0) delete [] B;
00185         if (X != 0) delete [] X;
00186         if (colStartA != 0) delete [] colStartA;
00187 
00188         // create the new
00189         B = new double[size];
00190         X = new double[size];
00191         colStartA = new int[size+1]; 
00192         
00193         if (B == 0 || X == 0 || colStartA == 0) {
00194             opserr << "WARNING SparseGenColLinSOE::SparseGenColLinSOE :";
00195             opserr << " ran out of memory for vectors (size) (";
00196             opserr << size << ") \n";
00197             size = 0; Bsize = 0;
00198             result =  -1;
00199         }
00200         else
00201             Bsize = size;
00202     }
00203 
00204     // zero the vectors
00205     for (int j=0; j<size; j++) {
00206         B[j] = 0;
00207         X[j] = 0;
00208     }
00209     
00210     // create new Vectors objects
00211     if (size != oldSize) {
00212         if (vectX != 0)
00213             delete vectX;
00214 
00215         if (vectB != 0)
00216             delete vectB;
00217         
00218         vectX = new Vector(X,size);
00219         vectB = new Vector(B,size);     
00220     }
00221 
00222     // fill in colStartA and rowA
00223     if (size != 0) {
00224       colStartA[0] = 0;
00225       int startLoc = 0;
00226       int lastLoc = 0;
00227       for (int a=0; a<size; a++) {
00228 
00229         theVertex = theGraph.getVertexPtr(a);
00230         if (theVertex == 0) {
00231           opserr << "WARNING:SparseGenColLinSOE::setSize :";
00232           opserr << " vertex " << a << " not in graph! - size set to 0\n";
00233           size = 0;
00234           return -1;
00235         }
00236 
00237         rowA[lastLoc++] = theVertex->getTag(); // place diag in first
00238         const ID &theAdjacency = theVertex->getAdjacency();
00239         int idSize = theAdjacency.Size();
00240         
00241         // now we have to place the entries in the ID into order in rowA
00242         for (int i=0; i<idSize; i++) {
00243 
00244           int row = theAdjacency(i);
00245           bool foundPlace = false;
00246           // find a place in rowA for current col
00247           for (int j=startLoc; j<lastLoc; j++)
00248             if (rowA[j] > row) { 
00249               // move the entries already there one further on
00250               // and place col in current location
00251               for (int k=lastLoc; k>j; k--)
00252                 
00253                 rowA[k] = rowA[k-1];
00254               rowA[j] = row;
00255               foundPlace = true;
00256               j = lastLoc;
00257             }
00258           if (foundPlace == false) // put in at the end
00259             rowA[lastLoc] = row;
00260 
00261           lastLoc++;
00262         }
00263         colStartA[a+1] = lastLoc;;          
00264         startLoc = lastLoc;
00265       }
00266     }
00267 
00268     
00269     // invoke setSize() on the Solver    
00270     LinearSOESolver *the_Solver = this->getSolver();
00271     int solverOK = the_Solver->setSize();
00272     if (solverOK < 0) {
00273         opserr << "WARNING:SparseGenColLinSOE::setSize :";
00274         opserr << " solver failed setSize()\n";
00275         return solverOK;
00276     }    
00277 
00278     return result;
00279 }
00280 
00281 int 
00282 SparseGenColLinSOE::addA(const Matrix &m, const ID &id, double fact)
00283 {
00284     // check for a quick return 
00285     if (fact == 0.0)  
00286         return 0;
00287 
00288     int idSize = id.Size();
00289     
00290     // check that m and id are of similar size
00291     if (idSize != m.noRows() && idSize != m.noCols()) {
00292         opserr << "SparseGenColLinSOE::addA() ";
00293         opserr << " - Matrix and ID not of similar sizes\n";
00294         return -1;
00295     }
00296     
00297     if (fact == 1.0) { // do not need to multiply 
00298       for (int i=0; i<idSize; i++) {
00299         int col = id(i);
00300         if (col < size && col >= 0) {
00301           int startColLoc = colStartA[col];
00302           int endColLoc = colStartA[col+1];
00303           for (int j=0; j<idSize; j++) {
00304             int row = id(j);
00305             if (row <size && row >= 0) {
00306               // find place in A using rowA
00307               for (int k=startColLoc; k<endColLoc; k++)
00308                 if (rowA[k] == row) {
00309                   A[k] += m(j,i);
00310                   k = endColLoc;
00311                 }
00312             }
00313           }  // for j           
00314         } 
00315       }  // for i
00316     } else {
00317       for (int i=0; i<idSize; i++) {
00318         int col = id(i);
00319         if (col < size && col >= 0) {
00320           int startColLoc = colStartA[col];
00321           int endColLoc = colStartA[col+1];
00322           for (int j=0; j<idSize; j++) {
00323             int row = id(j);
00324             if (row <size && row >= 0) {
00325               // find place in A using rowA
00326               for (int k=startColLoc; k<endColLoc; k++)
00327                 if (rowA[k] == row) {
00328                   A[k] += fact * m(j,i);
00329                   k = endColLoc;
00330                 }
00331             }
00332           }  // for j           
00333         } 
00334       }  // for i
00335     }
00336     return 0;
00337 }
00338 
00339     
00340 int 
00341 SparseGenColLinSOE::addB(const Vector &v, const ID &id, double fact)
00342 {
00343     // check for a quick return 
00344     if (fact == 0.0)  return 0;
00345 
00346     int idSize = id.Size();    
00347     // check that m and id are of similar size
00348     if (idSize != v.Size() ) {
00349         opserr << "SparseGenColLinSOE::addB() ";
00350         opserr << " - Vector and ID not of similar sizes\n";
00351         return -1;
00352     }    
00353 
00354     if (fact == 1.0) { // do not need to multiply if fact == 1.0
00355         for (int i=0; i<idSize; i++) {
00356             int pos = id(i);
00357             if (pos <size && pos >= 0)
00358                 B[pos] += v(i);
00359         }
00360     } else if (fact == -1.0) { // do not need to multiply if fact == -1.0
00361         for (int i=0; i<idSize; i++) {
00362             int pos = id(i);
00363             if (pos <size && pos >= 0)
00364                 B[pos] -= v(i);
00365         }
00366     } else {
00367         for (int i=0; i<idSize; i++) {
00368             int pos = id(i);
00369             if (pos <size && pos >= 0)
00370                 B[pos] += v(i) * fact;
00371         }
00372     }   
00373 
00374     return 0;
00375 }
00376 
00377 
00378 int
00379 SparseGenColLinSOE::setB(const Vector &v, double fact)
00380 {
00381     // check for a quick return 
00382     if (fact == 0.0)  return 0;
00383 
00384 
00385     if (v.Size() != size) {
00386         opserr << "WARNING BandGenLinSOE::setB() -";
00387         opserr << " incomptable sizes " << size << " and " << v.Size() << endln;
00388         return -1;
00389     }
00390     
00391     if (fact == 1.0) { // do not need to multiply if fact == 1.0
00392         for (int i=0; i<size; i++) {
00393             B[i] = v(i);
00394         }
00395     } else if (fact == -1.0) {
00396         for (int i=0; i<size; i++) {
00397             B[i] = -v(i);
00398         }
00399     } else {
00400         for (int i=0; i<size; i++) {
00401             B[i] = v(i) * fact;
00402         }
00403     }   
00404     return 0;
00405 }
00406 
00407 void 
00408 SparseGenColLinSOE::zeroA(void)
00409 {
00410     double *Aptr = A;
00411     for (int i=0; i<Asize; i++)
00412         *Aptr++ = 0;
00413 
00414     factored = false;
00415 }
00416         
00417 void 
00418 SparseGenColLinSOE::zeroB(void)
00419 {
00420     double *Bptr = B;
00421     for (int i=0; i<size; i++)
00422         *Bptr++ = 0;
00423 }
00424 
00425 void 
00426 SparseGenColLinSOE::setX(int loc, double value)
00427 {
00428     if (loc < size && loc >=0)
00429         X[loc] = value;
00430 }
00431 
00432 void 
00433 SparseGenColLinSOE::setX(const Vector &x)
00434 {
00435   if (x.Size() == size && vectX != 0)
00436     *vectX = x;
00437 }
00438 
00439 const Vector &
00440 SparseGenColLinSOE::getX(void)
00441 {
00442     if (vectX == 0) {
00443         opserr << "FATAL SparseGenColLinSOE::getX - vectX == 0";
00444         exit(-1);
00445     }
00446     return *vectX;
00447 }
00448 
00449 const Vector &
00450 SparseGenColLinSOE::getB(void)
00451 {
00452     if (vectB == 0) {
00453         opserr << "FATAL SparseGenColLinSOE::getB - vectB == 0";
00454         exit(-1);
00455     }        
00456     return *vectB;
00457 }
00458 
00459 double 
00460 SparseGenColLinSOE::normRHS(void)
00461 {
00462     double norm =0.0;
00463     for (int i=0; i<size; i++) {
00464         double Yi = B[i];
00465         norm += Yi*Yi;
00466     }
00467     return sqrt(norm);
00468     
00469 }    
00470 
00471 
00472 int
00473 SparseGenColLinSOE::setSparseGenColSolver(SparseGenColLinSolver &newSolver)
00474 {
00475     newSolver.setLinearSOE(*this);
00476     
00477     if (size != 0) {
00478         int solverOK = newSolver.setSize();
00479         if (solverOK < 0) {
00480             opserr << "WARNING:SparseGenColLinSOE::setSolver :";
00481             opserr << "the new solver could not setSeize() - staying with old\n";
00482             return -1;
00483         }
00484     }
00485     
00486     return this->LinearSOE::setSolver(newSolver);
00487 }
00488 
00489 
00490 int 
00491 SparseGenColLinSOE::sendSelf(int cTag, Channel &theChannel)
00492 {
00493     return 0;
00494 }
00495 
00496 int 
00497 SparseGenColLinSOE::recvSelf(int cTag, Channel &theChannel, 
00498                              FEM_ObjectBroker &theBroker)  
00499 {
00500     return 0;
00501 }
00502 

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