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

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