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/profileSPD/ProfileSPDLinDirectBlockSolver.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/system_of_eqn/linearSOE/ProfileSPD/ProfileSPDLinDirectBlockSolver.C
00027 //
00028 // Written: fmk 
00029 // Created: Mar 1998
00030 // Revision: A
00031 //
00032 // Description: This file contains the class definition for 
00033 // ProfileSPDLinDirectBlockSolver. ProfileSPDLinDirectBlockSolver is a subclass 
00034 // of LinearSOESOlver. It solves a ProfileSPDLinSOE object using
00035 // the LDL^t factorization and a block approach
00036 
00037 // What: "@(#) ProfileSPDLinDirectBlockSolver.C, revA"
00038 
00039 #include <ProfileSPDLinDirectBlockSolver.h>
00040 #include <ProfileSPDLinSOE.h>
00041 #include <math.h>
00042 #include <Channel.h>
00043 #include <FEM_ObjectBroker.h>
00044 
00045 #include <Timer.h>
00046 
00047 ProfileSPDLinDirectBlockSolver::ProfileSPDLinDirectBlockSolver(double tol, int blckSize)
00048 :ProfileSPDLinSolver(SOLVER_TAGS_ProfileSPDLinDirectBlockSolver),
00049  minDiagTol(tol), blockSize(blckSize), maxColHeight(0), 
00050  size(0), RowTop(0), topRowPtr(0), invD(0)
00051 {
00052 
00053 }
00054 
00055     
00056 ProfileSPDLinDirectBlockSolver::~ProfileSPDLinDirectBlockSolver()
00057 {
00058     if (RowTop != 0) delete [] RowTop;
00059     if (topRowPtr != 0) free((void *)topRowPtr);
00060     if (invD != 0) delete [] invD;
00061 }
00062 
00063 int
00064 ProfileSPDLinDirectBlockSolver::setSize(void)
00065 {
00066 
00067     if (theSOE == 0) {
00068  cerr << "ProfileSPDLinDirectBlockSolver::setSize()";
00069  cerr << " No system has been set\n";
00070  return -1;
00071     }
00072 
00073     // check for quick return 
00074     if (theSOE->size == 0)
00075  return 0;
00076     if (size != theSOE->size) {    
00077       size = theSOE->size;
00078     
00079       if (RowTop != 0) delete [] RowTop;
00080       if (topRowPtr != 0) delete [] topRowPtr;
00081       if (invD != 0) delete [] invD;
00082 
00083       RowTop = new int[size];
00084 
00085       // we cannot use topRowPtr = new (double *)[size] with the cxx compiler
00086       topRowPtr = (double **)malloc(size *sizeof(double *));
00087 
00088       invD = new double[size]; 
00089  
00090       if (RowTop == 0 || topRowPtr == 0 || invD == 0) {
00091  cerr << "Warning :ProfileSPDLinDirectBlockSolver::ProfileSPDLinDirectBlockSolver :";
00092  cerr << " ran out of memory for work areas \n";
00093  return -1;
00094       }
00095     }
00096 
00097 
00098     // set some pointers
00099     double *A = theSOE->A;
00100     int *iDiagLoc = theSOE->iDiagLoc;
00101 
00102     // set RowTop and topRowPtr info
00103 
00104     maxColHeight = 0;
00105     RowTop[0] = 0;
00106     topRowPtr[0] = A;
00107     for (int j=1; j<size; j++) {
00108  int icolsz = iDiagLoc[j] - iDiagLoc[j-1];
00109         if (icolsz > maxColHeight) maxColHeight = icolsz;
00110  RowTop[j] = j - icolsz +  1;
00111  topRowPtr[j] = &A[iDiagLoc[j-1]]; // FORTRAN array indexing in iDiagLoc
00112     }
00113 
00114     size = theSOE->size;
00115     return 0;
00116 }
00117 
00118 
00119 int 
00120 ProfileSPDLinDirectBlockSolver::solve(void)
00121 {
00122     // check for quick returns
00123     if (theSOE == 0) {
00124  cerr << "ProfileSPDLinDirectBlockSolver::solve(void): ";
00125  cerr << " - No ProfileSPDSOE has been assigned\n";
00126  return -1;
00127     }
00128     
00129     if (theSOE->size == 0)
00130  return 0;
00131 
00132     // set some pointers
00133     double *B = theSOE->B;
00134     double *X = theSOE->X;
00135     int n = theSOE->size;
00136 
00137     // copy B into X
00138     for (int ii=0; ii<n; ii++)
00139  X[ii] = B[ii];
00140     
00141     if (theSOE->isAfactored == false)  {
00142 
00143  // FACTOR 
00144  invD[0] = 1.0/theSOE->A[0]; 
00145         int startRow = 0;
00146  int lastRow = startRow+blockSize-1;
00147         int lastColEffected = lastRow+maxColHeight -1;
00148         int nBlck = n/blockSize;
00149  if ((n % blockSize) != 0)
00150    nBlck++;
00151 
00152  // for every block across      
00153  for (int i=0; i<nBlck; i++) {
00154 
00155    // first factor the diagonal block int Ui,i and Di
00156           int j;
00157    for (j=0; j<blockSize; j++) {
00158             int currentRow = startRow + j;
00159 
00160      if (currentRow < n) { // this is for case when size%blockSize != 0
00161 
00162        int rowjTop = RowTop[currentRow];
00163        double *akjPtr = topRowPtr[currentRow];
00164        int maxRowijTop;
00165        if (rowjTop < startRow) {
00166   akjPtr += startRow-rowjTop; // pointer to start of block row
00167   maxRowijTop = startRow;
00168        } else
00169   maxRowijTop = rowjTop;
00170 
00171        int k;
00172        for (k=maxRowijTop; k<currentRow; k++) {
00173   double tmp = *akjPtr;
00174   int rowkTop = RowTop[k];
00175   int maxRowkjTop;
00176   double *alkPtr, *aljPtr;
00177   if (rowkTop < rowjTop) {
00178     alkPtr = topRowPtr[k] + (rowjTop - rowkTop);
00179     aljPtr = topRowPtr[currentRow];
00180     maxRowkjTop = rowjTop;
00181   } else {
00182     alkPtr = topRowPtr[k];
00183     aljPtr = topRowPtr[currentRow] + (rowkTop - rowjTop);
00184     maxRowkjTop = rowkTop;
00185   }
00186 
00187   for (int l = maxRowkjTop; l<k; l++) 
00188     tmp -= *alkPtr++ * *aljPtr++;
00189   
00190   *akjPtr++ = tmp;
00191        }
00192 
00193        double ajj = *akjPtr;
00194        akjPtr = topRowPtr[currentRow];
00195        double *bjPtr  = &X[rowjTop];  
00196        double tmp = 0;     
00197 
00198        for (k=rowjTop; k<currentRow; k++){
00199   double akj = *akjPtr;
00200   double lkj = akj * invD[k];
00201   tmp -= lkj * *bjPtr++;   
00202   *akjPtr++ = lkj;
00203   ajj = ajj -lkj * akj;
00204        }
00205 
00206        X[currentRow] += tmp;
00207 
00208        // check that the diag > the tolerance specified
00209        if (ajj <= 0.0) {
00210   cerr << "ProfileSPDLinDirectBlockSolver::solve() - ";
00211   cerr << " aii < 0 (i, aii): (" << currentRow << ", " << ajj << ")\n"; 
00212   return(-2);
00213        }
00214        if (ajj <= minDiagTol) {
00215   cerr << "ProfileSPDLinDirectBlockSolver::solve() - ";
00216   cerr << " aii < minDiagTol (i, aii): (" << currentRow;
00217   cerr << ", " << ajj << ")\n"; 
00218   return(-2);
00219        }  
00220        invD[currentRow] = 1.0/ajj; 
00221 
00222      } else 
00223        j = blockSize;
00224 
00225    }
00226 
00227    // now do rest of i'th block row doing a block of columns at a time
00228           // forming Ui,j*Di
00229    int currentCol = startRow + blockSize;
00230    for (j=i+1; j<nBlck; j++) 
00231      for (int k=0; k<blockSize; k++) {
00232 
00233        if (currentCol < n) { // this is for case when size%blockSize != 0
00234 
00235   int rowkTop = RowTop[currentCol];
00236   double *alkPtr = topRowPtr[currentCol];
00237   int maxRowikTop;
00238   if (rowkTop < startRow) {
00239     alkPtr += startRow-rowkTop; // pointer to start of block row
00240     maxRowikTop = startRow;
00241   } else
00242     maxRowikTop = rowkTop;
00243 
00244   for (int l=maxRowikTop; l<=lastRow; l++) {
00245     double tmp = *alkPtr;
00246     int rowlTop = RowTop[l];
00247     int maxRowklTop;
00248     double *amlPtr, *amkPtr;
00249     if (rowlTop < rowkTop) {
00250       amlPtr = topRowPtr[l] + (rowkTop - rowlTop);
00251       amkPtr = topRowPtr[currentCol];
00252       maxRowklTop = rowkTop;
00253     } else {
00254       amlPtr = topRowPtr[l];
00255       amkPtr = topRowPtr[currentCol] + (rowlTop - rowkTop);
00256       maxRowklTop = rowlTop;
00257     }
00258     
00259     for (int m = maxRowklTop; m<l; m++) 
00260       tmp -= *amkPtr++ * *amlPtr++;
00261     
00262     *alkPtr++ = tmp;
00263   }
00264   currentCol++;
00265   if (currentCol > lastColEffected) {
00266     k = blockSize;
00267     j = nBlck;
00268   }
00269 
00270        } else
00271   k = blockSize;
00272 
00273      }
00274 
00275    startRow += blockSize;
00276    lastRow = startRow + blockSize -1;
00277           lastColEffected = lastRow + maxColHeight -1;
00278  }
00279 
00280  theSOE->isAfactored = true;
00281  theSOE->numInt = 0;
00282  
00283  // divide by diag term 
00284  double *bjPtr = X; 
00285  double *aiiPtr = invD;
00286  for (int j=0; j<n; j++) 
00287    *bjPtr++ = *aiiPtr++ * X[j];
00288 
00289     
00290  // now do the back substitution storing result in X
00291  for (int k=(n-1); k>0; k--) {
00292       
00293    int rowktop = RowTop[k];
00294    double bk = X[k];
00295    double *ajiPtr = topRowPtr[k];   
00296 
00297    for (int j=rowktop; j<k; j++) 
00298      X[j] -= *ajiPtr++ * bk;
00299  }     
00300     } else { // just do forward and back substitution
00301       
00302       // do forward substitution 
00303       for (int i=1; i<n; i++) {
00304      
00305  int rowitop = RowTop[i];     
00306  double *ajiPtr = topRowPtr[i];
00307  double *bjPtr  = &X[rowitop];  
00308  double tmp = 0;     
00309      
00310  for (int j=rowitop; j<i; j++) 
00311    tmp -= *ajiPtr++ * *bjPtr++; 
00312      
00313  X[i] += tmp;
00314       }
00315 
00316       // divide by diag term 
00317       double *bjPtr = X; 
00318       double *aiiPtr = invD;
00319       for (int j=0; j<n; j++) 
00320  *bjPtr++ = *aiiPtr++ * X[j];
00321 
00322     
00323       // now do the back substitution storing result in X
00324       for (int k=(n-1); k>0; k--) {
00325       
00326  int rowktop = RowTop[k];
00327  double bk = X[k];
00328  double *ajiPtr = topRowPtr[k];   
00329 
00330  for (int j=rowktop; j<k; j++) 
00331    X[j] -= *ajiPtr++ * bk;
00332       }     
00333     }    
00334     return 0;
00335 }
00336 
00337 int 
00338 ProfileSPDLinDirectBlockSolver::setProfileSOE(ProfileSPDLinSOE &theNewSOE)
00339 {
00340     if (theSOE != 0) {
00341  cerr << "ProfileSPDLinDirectBlockSolver::setProfileSOE() - ";
00342  cerr << " has already been called \n"; 
00343  return -1;
00344     }
00345     
00346     theSOE = &theNewSOE;
00347     return 0;
00348 }
00349  
00350 int
00351 ProfileSPDLinDirectBlockSolver::sendSelf(int cTag, Channel &theChannel)
00352 {
00353     if (size != 0)
00354  cerr << "ProfileSPDLinDirectBlockSolver::sendSelf - does not send itself YET\n"; 
00355     return 0;
00356 }
00357 
00358 
00359 int 
00360 ProfileSPDLinDirectBlockSolver::recvSelf(int cTag,
00361       Channel &theChannel, 
00362       FEM_ObjectBroker &theBroker)
00363 {
00364     return 0;
00365 }
00366 
00367