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