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