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