My title

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