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

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