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

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