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

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