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.2 $
00022 // $Date: 2001/02/17 06:32:38 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/profileSPD/ProfileSPDLinDirectSolver.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/system_of_eqn/linearSOE/ProfileSPD/ProfileSPDLinSOESolver.C
00027 //
00028 // Written: fmk 
00029 // Created: Febuary 1997
00030 // Revision: A
00031 //
00032 // Description: This file contains the implementation for ProfileSPDLinSOESolver
00033 // Description: This file contains the class definition for 
00034 // ProfileSPDLinDirectSolver. ProfileSPDLinDirectSolver is a subclass 
00035 // of LinearSOESOlver. It solves a ProfileSPDLinSOE object using
00036 // the LDL^t factorization.
00037 
00038 // What: "@(#) ProfileSPDLinDirectSolver.C, revA"
00039 
00040 #include <ProfileSPDLinDirectSolver.h>
00041 #include <ProfileSPDLinSOE.h>
00042 #include <math.h>
00043 #include <Channel.h>
00044 #include <FEM_ObjectBroker.h>
00045 
00046 ProfileSPDLinDirectSolver::ProfileSPDLinDirectSolver(double tol)
00047 :ProfileSPDLinSolver(SOLVER_TAGS_ProfileSPDLinDirectSolver),
00048  minDiagTol(tol), size(0), RowTop(0), topRowPtr(0), invD(0)
00049 {
00050 
00051 }
00052 
00053     
00054 ProfileSPDLinDirectSolver::~ProfileSPDLinDirectSolver()
00055 {
00056     if (RowTop != 0) delete [] RowTop;
00057     if (topRowPtr != 0) free((void *)topRowPtr);
00058     if (invD != 0) delete [] invD;
00059 }
00060 
00061 int
00062 ProfileSPDLinDirectSolver::setSize(void)
00063 {
00064 
00065     if (theSOE == 0) {
00066  cerr << "ProfileSPDLinDirectSolver::setSize()";
00067  cerr << " No system has been set\n";
00068  return -1;
00069     }
00070 
00071     // check for quick return 
00072     if (theSOE->size == 0)
00073  return 0;
00074     
00075     size = theSOE->size;
00076     
00077     if (RowTop != 0) delete [] RowTop;
00078     if (topRowPtr != 0) delete [] topRowPtr;
00079     if (invD != 0) delete [] invD;
00080 
00081     RowTop = new int[size];
00082 
00083     // we cannot use topRowPtr = new (double *)[size] with the cxx compiler
00084     topRowPtr = (double **)malloc(size *sizeof(double *));
00085 
00086     invD = new double[size]; 
00087  
00088     if (RowTop == 0 || topRowPtr == 0 || invD == 0) {
00089  cerr << "Warning :ProfileSPDLinDirectSolver::ProfileSPDLinDirectSolver :";
00090  cerr << " ran out of memory for work areas \n";
00091  return -1;
00092     }
00093 
00094 
00095     // set some pointers
00096     double *A = theSOE->A;
00097     int *iDiagLoc = theSOE->iDiagLoc;
00098 
00099     // set RowTop and topRowPtr info
00100 
00101     RowTop[0] = 0;
00102     topRowPtr[0] = A;
00103     for (int j=1; j<size; j++) {
00104  int icolsz = iDiagLoc[j] - iDiagLoc[j-1];
00105  RowTop[j] = j - icolsz +  1;
00106  topRowPtr[j] = &A[iDiagLoc[j-1]]; // FORTRAN array indexing in iDiagLoc
00107     }
00108 
00109     size = theSOE->size;
00110     return 0;
00111 }
00112 
00113 
00114 int 
00115 ProfileSPDLinDirectSolver::solve(void)
00116 {
00117 
00118     // check for quick returns
00119     if (theSOE == 0) {
00120  cerr << "ProfileSPDLinDirectSolver::solve(void): ";
00121  cerr << " - No ProfileSPDSOE has been assigned\n";
00122  return -1;
00123     }
00124     
00125     if (theSOE->size == 0)
00126  return 0;
00127 
00128     // set some pointers
00129     double *B = theSOE->B;
00130     double *X = theSOE->X;
00131     int theSize = theSOE->size;
00132     // copy B into X
00133     for (int ii=0; ii<theSize; ii++)
00134  X[ii] = B[ii];
00135 
00136     /*
00137       for (int iii=0; iii<theSize; iii++) {
00138       int rowiiitop = RowTop[iii];
00139       double *ajiptr = topRowPtr[iii];
00140       cerr << "\n COLUMN " << iii << " TopRow " << rowiiitop << " -> ";
00141       for (int jjj = rowiiitop; jjj <=iii; jjj++)
00142       cerr << *ajiptr++ << " ";
00143       }
00144       cerr << endl;
00145 
00146       for (int iii=0; iii<theSOE->size; iii++) {
00147       cerr << "COLUMN " << iii << " Biii -> " << X[iii] << endl;
00148       }
00149       cerr << endl;
00150       */
00151 
00152     
00153     if (theSOE->isAfactored == false)  {
00154 
00155  // FACTOR & SOLVE
00156  double *ajiPtr, *akjPtr, *akiPtr, *bjPtr;    
00157  
00158  // if the matrix has not been factored already factor it into U^t D U
00159  // storing D^-1 in invD as we go
00160 
00161  double a00 = theSOE->A[0];
00162  if (a00 <= 0.0) {
00163    cerr << "ProfileSPDLinDirectSolver::solve() - ";
00164    cerr << " aii < 0 (i, aii): (0,0)\n"; 
00165    return(-2);
00166  }    
00167  
00168         invD[0] = 1.0/theSOE->A[0]; 
00169  
00170  // for every col across 
00171  for (int i=1; i<theSize; i++) {
00172 
00173      int rowitop = RowTop[i];
00174      ajiPtr = topRowPtr[i];
00175 
00176      for (int j=rowitop; j<i; j++) {
00177   double tmp = *ajiPtr;
00178   int rowjtop = RowTop[j];
00179 
00180   if (rowitop > rowjtop) {
00181 
00182       akjPtr = topRowPtr[j] + (rowitop-rowjtop);
00183       akiPtr = topRowPtr[i];
00184 
00185       for (int k=rowitop; k<j; k++) 
00186    tmp -= *akjPtr++ * *akiPtr++ ;
00187 
00188       *ajiPtr++ = tmp;
00189   }
00190   else {
00191       akjPtr = topRowPtr[j];
00192       akiPtr = topRowPtr[i] + (rowjtop-rowitop);
00193 
00194       for (int k=rowjtop; k<j; k++) 
00195    tmp -= *akjPtr++ * *akiPtr++ ;
00196 
00197       *ajiPtr++ = tmp;
00198   }
00199      }
00200 
00201      /* now form i'th col of [U] and determine [dii] */
00202 
00203      double aii = theSOE->A[theSOE->iDiagLoc[i] -1]; // FORTRAN ARRAY INDEXING
00204      ajiPtr = topRowPtr[i];
00205      double *bjPtr  = &X[rowitop];  
00206      double tmp = 0;     
00207      
00208      for (int jj=rowitop; jj<i; jj++) {
00209   double aji = *ajiPtr;
00210   double lij = aji * invD[jj];
00211   tmp -= lij * *bjPtr++;   
00212   *ajiPtr++ = lij;
00213   aii = aii - lij*aji;
00214      }
00215      
00216      // check that the diag > the tolerance specified
00217      if (aii == 0.0) {
00218   cerr << "ProfileSPDLinDirectSolver::solve() - ";
00219   cerr << " aii < 0 (i, aii): (" << i << ", " << aii << ")\n"; 
00220   return(-2);
00221      }
00222      if (fabs(aii) <= minDiagTol) {
00223   cerr << "ProfileSPDLinDirectSolver::solve() - ";
00224   cerr << " aii < minDiagTol (i, aii): (" << i;
00225   cerr << ", " << aii << ")\n"; 
00226   return(-2);
00227      }  
00228      invD[i] = 1.0/aii; 
00229      X[i] += tmp;     
00230  }
00231 
00232  theSOE->isAfactored = true;
00233  theSOE->numInt = 0;
00234  
00235  
00236  // divide by diag term 
00237  bjPtr = X; 
00238  double *aiiPtr = invD;
00239  for (int j=0; j<theSize; j++) 
00240      *bjPtr++ = *aiiPtr++ * X[j];
00241 
00242 
00243  // now do the back substitution storing result in X
00244  for (int k=(theSize-1); k>0; k--) {
00245 
00246      int rowktop = RowTop[k];
00247      double bk = X[k];
00248      double *ajiPtr = topRowPtr[k];   
00249 
00250      for (int j=rowktop; j<k; j++) 
00251   X[j] -= *ajiPtr++ * bk;
00252  }      
00253     }
00254 
00255     else {
00256 
00257  // JUST DO SOLVE
00258 
00259  // do forward substitution 
00260  for (int i=1; i<theSize; i++) {
00261      
00262      int rowitop = RowTop[i];     
00263      double *ajiPtr = topRowPtr[i];
00264      double *bjPtr  = &X[rowitop];  
00265      double tmp = 0;     
00266      
00267      for (int j=rowitop; j<i; j++) 
00268   tmp -= *ajiPtr++ * *bjPtr++; 
00269      
00270      X[i] += tmp;
00271  }
00272 
00273  // divide by diag term 
00274  double *bjPtr = X; 
00275  double *aiiPtr = invD;
00276  for (int j=0; j<theSize; j++) 
00277      *bjPtr++ = *aiiPtr++ * X[j];
00278 
00279 
00280  // now do the back substitution storing result in X
00281  for (int k=(theSize-1); k>0; k--) {
00282 
00283      int rowktop = RowTop[k];
00284      double bk = X[k];
00285      double *ajiPtr = topRowPtr[k];   
00286 
00287      for (int j=rowktop; j<k; j++) 
00288   X[j] -= *ajiPtr++ * bk;
00289  }     
00290     }    
00291     
00292     /*
00293     cerr << "BBBB " << theSOE->getB();
00294     cerr << "XXXX " << theSOE->getX();
00295     */
00296     
00297     return 0;
00298 }
00299 
00300 double
00301 ProfileSPDLinDirectSolver::getDeterminant(void) 
00302 {
00303    int theSize = theSOE->size;
00304    double determinant = 1.0;
00305    for (int i=0; i<theSize; i++)
00306      determinant *= invD[i];
00307    determinant = 1.0/determinant;
00308    return determinant;
00309 }
00310 
00311 int 
00312 ProfileSPDLinDirectSolver::setProfileSOE(ProfileSPDLinSOE &theNewSOE)
00313 {
00314     if (theSOE != 0) {
00315  cerr << "ProfileSPDLinDirectSolver::setProfileSOE() - ";
00316  cerr << " has already been called \n"; 
00317  return -1;
00318     }
00319     
00320     theSOE = &theNewSOE;
00321     return 0;
00322 }
00323  
00324 
00325 int 
00326 ProfileSPDLinDirectSolver::factor(int n)
00327 {
00328 
00329     // check for quick returns
00330     if (theSOE == 0) {
00331  cerr << "ProfileSPDLinDirectSolver::factor: ";
00332  cerr << " - No ProfileSPDSOE has been assigned\n";
00333  return -1;
00334     }
00335 
00336     int theSize = theSOE->size;    
00337     if (n > theSize) {
00338  cerr << "ProfileSPDLinDirectSolver::factor: ";
00339  cerr << " - n " << n << " greater than size of system" << theSize << endl;
00340  return -1;
00341     }
00342 
00343     if (theSize == 0 || n == 0)
00344  return 0;
00345 
00346 
00347     // set some pointers
00348     if (theSOE->isAfactored == false)  {
00349 
00350  // FACTOR & SOLVE
00351  double *ajiPtr, *akjPtr, *akiPtr;    
00352  
00353  // if the matrix has not been factored already factor it into U^t D U
00354  // storing D^-1 in invD as we go
00355     
00356  invD[0] = 1.0/theSOE->A[0]; 
00357  
00358  // for every col across 
00359  for (int i=1; i<n; i++) {
00360 
00361      int rowitop = RowTop[i];
00362      ajiPtr = topRowPtr[i];
00363 
00364      for (int j=rowitop; j<i; j++) {
00365   double tmp = *ajiPtr;
00366   int rowjtop = RowTop[j];
00367 
00368   if (rowitop > rowjtop) {
00369 
00370       akjPtr = topRowPtr[j] + (rowitop-rowjtop);
00371       akiPtr = topRowPtr[i];
00372 
00373       for (int k=rowitop; k<j; k++) 
00374    tmp -= *akjPtr++ * *akiPtr++ ;
00375 
00376       *ajiPtr++ = tmp;
00377   }
00378   else {
00379       akjPtr = topRowPtr[j];
00380       akiPtr = topRowPtr[i] + (rowjtop-rowitop);
00381 
00382       for (int k=rowjtop; k<j; k++) 
00383    tmp -= *akjPtr++ * *akiPtr++ ;
00384 
00385       *ajiPtr++ = tmp;
00386   }
00387      }
00388 
00389      /* now form i'th col of [U] and determine [dii] */
00390 
00391      double aii = theSOE->A[theSOE->iDiagLoc[i] -1]; // FORTRAN ARRAY INDEXING
00392      ajiPtr = topRowPtr[i];
00393      
00394      for (int jj=rowitop; jj<i; jj++) {
00395   double aji = *ajiPtr;
00396   double lij = aji * invD[jj];
00397   *ajiPtr++ = lij;
00398   aii = aii - lij*aji;
00399      }
00400      
00401      // check that the diag > the tolerance specified
00402      if (aii <= 0.0) {
00403   cerr << "ProfileSPDLinDirectSolver::solve() - ";
00404   cerr << " aii < 0 (i, aii): (" << i << ", " << aii << ")\n"; 
00405   return(-2);
00406      }
00407      if (aii <= minDiagTol) {
00408   cerr << "ProfileSPDLinDirectSolver::solve() - ";
00409   cerr << " aii < minDiagTol (i, aii): (" << i;
00410   cerr << ", " << aii << ")\n"; 
00411   return(-2);
00412      }  
00413      invD[i] = 1.0/aii; 
00414  }
00415 
00416  theSOE->isAfactored = true;
00417  theSOE->numInt = n;
00418  
00419     } 
00420     return 0;
00421 }
00422 
00423 
00424 /*
00425     double *ajiPtr, *akjPtr, *akiPtr, *bjPtr;    
00426  
00427     // if the matrix has not been factored already factor it into U^t D U
00428     // storing D^-1 in invD as we go
00429     
00430     if (theSOE->isAfactored == false) {
00431  
00432  invD[0] = 1.0/theSOE->A[0]; 
00433  
00434  // for every col across 
00435  for (int i=1; i<n; i++) {
00436 
00437      int rowitop = RowTop[i];
00438      ajiPtr = topRowPtr[i];
00439 
00440      for (int j=rowitop; j<i; j++) {
00441   double tmp = *ajiPtr;
00442   int rowjtop = RowTop[j];
00443 
00444   if (rowitop > rowjtop) {
00445 
00446       akjPtr = topRowPtr[j] + (rowitop-rowjtop);
00447       akiPtr = topRowPtr[i];
00448 
00449       for (int k=rowitop; k<j; k++) 
00450    tmp -= *akjPtr++ * *akiPtr++ ;
00451 
00452       *ajiPtr++ = tmp;
00453   }
00454   else {
00455       akjPtr = topRowPtr[j];
00456       akiPtr = topRowPtr[i] + (rowjtop-rowitop);
00457 
00458       for (int k=rowjtop; k<j; k++) 
00459    tmp -= *akjPtr++ * *akiPtr++ ;
00460 
00461       *ajiPtr++ = tmp;
00462   }
00463      }
00464 
00465      // now form i'th col of [U] and determine [dii] 
00466 
00467      double aii = theSOE->A[theSOE->iDiagLoc[i] -1]; // FORTRAN ARRAY INDEXING
00468      ajiPtr = topRowPtr[i];
00469      
00470      for (int jj=rowitop; jj<i; jj++) {
00471   double aji = *ajiPtr;
00472   double lij = aji * invD[jj];
00473   *ajiPtr++ = lij;
00474   aii = aii - lij*aji;
00475      }
00476      
00477      // check that the diag > the tolerance specified
00478      if (aii <= 0.0)     return(-2);
00479      if (aii <= minDiagTol) return(-1);
00480      invD[i] = 1.0/aii; 
00481  }
00482     }
00483     
00484     theSOE->isAfactored = true;
00485     if (n == theSize)
00486  theSOE->numInt = 0;
00487     else
00488  theSOE->numInt = n;
00489     
00490     return 0;
00491 }
00492 */
00493 
00494 int
00495 ProfileSPDLinDirectSolver::sendSelf(int cTag,
00496         Channel &theChannel)
00497 {
00498     return 0;
00499 }
00500 
00501 
00502 int 
00503 ProfileSPDLinDirectSolver::recvSelf(int cTag,
00504         Channel &theChannel, 
00505         FEM_ObjectBroker &theBroker)
00506 {
00507     return 0;
00508 }
00509 
00510