OpenSees » SRC » system_of_eqn » linearSOE » itpack
ItpackLinSolver.cpp

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.3 $
00022 // $Date: 2003/02/14 23:02:02 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/itpack/ItpackLinSolver.cpp,v $
00024                                                                         
00025 // Written: MHS
00026 // Created: Sept 2001
00027 //
00028 // Description: This file contains the class definition for ItpackLinSolver.
00029 // ItpackLinSolver is an abstract base class and thus no objects of it's type
00030 // can be instantiated. It has pure virtual functions which must be
00031 // implemented in it's derived classes.  Instances of ItpackLinSolver 
00032 // are used to solve a system of equations of type ItpackLinSOE.
00033 
00034 #include <ItpackLinSolver.h>
00035 #include <ItpackLinSOE.h>
00036 #include <ID.h>
00037 
00038 ItpackLinSolver::ItpackLinSolver(int meth, int iter, double om)
00039   :LinearSOESolver(SOLVER_TAGS_Itpack),
00040    theSOE(0), IA(0), JA(0), n(0),
00041    iwksp(0), wksp(0), nwksp(0), maxIter(iter),
00042    method(meth), omega(om)
00043 {
00044 
00045 }    
00046 
00047 ItpackLinSolver::ItpackLinSolver()
00048   :LinearSOESolver(SOLVER_TAGS_Itpack),
00049    theSOE(0), IA(0), JA(0), n(0),
00050    iwksp(0), wksp(0), nwksp(0), maxIter(0),
00051    method(0), omega(0.0)
00052 {
00053 
00054 }    
00055 
00056 ItpackLinSolver::~ItpackLinSolver()    
00057 {
00058   if (IA != 0)
00059     delete [] IA;
00060 
00061   if (JA != 0)
00062     delete [] JA;
00063 
00064   if (iwksp != 0)
00065     delete [] iwksp;
00066 
00067   if (wksp != 0)
00068     delete [] wksp;
00069 }    
00070 
00071 int 
00072 ItpackLinSolver::setLinearSOE(ItpackLinSOE &theItpackSOE)
00073 {
00074   theSOE = &theItpackSOE;
00075 
00076   return 0;
00077 }
00078 
00079 int
00080 ItpackLinSolver::setSize(void)
00081 {
00082   // Get number of equations from SOE
00083   n = theSOE->size;
00084 
00085   if (n > 0) {
00086     if (iwksp != 0)
00087       delete [] iwksp;
00088     iwksp = new int[3*n];
00089   }
00090 
00091   // Should be 2*maxIter for symmetric storage, 4*maxIter for nonsymmetric
00092   int ncg = 4*maxIter;
00093   
00094   // Order of the black subsystem
00095   int nb = theSOE->nnz; // I think this is the most it could be
00096 
00097   switch(method) {
00098   case ItpackJCG:
00099     nwksp = 4*n + ncg;
00100     break;
00101   case ItpackJSI: case ItpackJ:
00102     nwksp = 2*n;
00103     break;
00104   case ItpackSOR: case ItpackGS: case ItpackSORFixed:
00105     nwksp = n;
00106     break;
00107   case ItpackSSORCG:
00108     nwksp = 6*n + ncg;
00109     break;
00110   case ItpackSSORSI: case ItpackSSORFixed:
00111     nwksp = 5*n;
00112     break;
00113   case ItpackRSCG:
00114     nwksp = n + 3*nb + ncg;
00115     break;
00116   case ItpackRSSI: case ItpackRS:
00117     nwksp = n + nb;
00118     break;
00119   default:
00120     nwksp = 6*n + ncg;
00121     break;
00122   }
00123 
00124   if (nwksp > 0) {
00125     if (wksp != 0)
00126       delete [] wksp;
00127     wksp = new double[nwksp];
00128   }
00129 
00130   // Get number of nonzeros from the SOE
00131   int nnz = theSOE->nnz;
00132   
00133   if (nnz > 0) {
00134     if (JA != 0)
00135       delete [] JA;
00136     JA = new int [nnz];
00137   }
00138 
00139   int *jaPtr = theSOE->colA;
00140   int i;
00141   // Add one for FORTRAN indexing
00142   for (i = 0; i < nnz; i++)
00143     JA[i] = jaPtr[i] + 1;
00144 
00145   if (n > 0) {
00146     if (IA != 0)
00147       delete [] IA;
00148     IA = new int [n+1];
00149   }
00150   
00151   int *iaPtr = theSOE->rowStartA;
00152   // Add one for FORTRAN indexing
00153   for (i = 0; i <= n; i++) 
00154     IA[i] = iaPtr[i] + 1;
00155 
00156   return 0;
00157 }
00158 
00159 int
00160 ItpackLinSolver::sendSelf(int commitTag, Channel &theChannel)
00161 {
00162   return -1;
00163 }
00164 
00165 int
00166 ItpackLinSolver::recvSelf(int commitTag, Channel &theChannel, 
00167                           FEM_ObjectBroker &theBroker)
00168 {
00169   return -1;
00170 }
00171 
00172 #ifdef _WIN32
00173 
00174 extern "C" int _stdcall DFAULT(int *IPARM, double *RPARM);
00175 
00176 extern "C" int _stdcall VFILL(int *N, double *U, double *VAL);
00177 
00178 extern "C" int _stdcall JCG(int *N, int *IA, int *JA, double *A, double *RHS,
00179                             double *U, int *IWKSP, int *NW, double *WKSP,
00180                             int *IPARM, double *RPARM, int *IER);
00181 
00182 extern "C" int _stdcall JSI(int *N, int *IA, int *JA, double *A, double *RHS,
00183                             double *U, int *IWKSP, int *NW, double *WKSP,
00184                             int *IPARM, double *RPARM, int *IER);
00185 
00186 extern "C" int _stdcall SOR(int *N, int *IA, int *JA, double *A, double *RHS,
00187                             double *U, int *IWKSP, int *NW, double *WKSP,
00188                             int *IPARM, double *RPARM, int *IER);
00189 
00190 extern "C" int _stdcall SSORCG(int *N, int *IA, int *JA, double *A, double *RHS,
00191                                double *U, int *IWKSP, int *NW, double *WKSP,
00192                                int *IPARM, double *RPARM, int *IER);
00193 
00194 extern "C" int _stdcall SSORSI(int *N, int *IA, int *JA, double *A, double *RHS,
00195                                double *U, int *IWKSP, int *NW, double *WKSP,
00196                                int *IPARM, double *RPARM, int *IER);
00197 
00198 extern "C" int _stdcall RSCG(int *N, int *IA, int *JA, double *A, double *RHS,
00199                              double *U, int *IWKSP, int *NW, double *WKSP,
00200                              int *IPARM, double *RPARM, int *IER);
00201 
00202 extern "C" int _stdcall RSSI(int *N, int *IA, int *JA, double *A, double *RHS,
00203                              double *U, int *IWKSP, int *NW, double *WKSP,
00204                              int *IPARM, double *RPARM, int *IER);
00205 
00206 #define DFAULT dfault_
00207 #define VFILL  vfill_
00208 #define JCG    jcg_
00209 #define JSI    jsi_
00210 #define SOR    sor_
00211 #define SSORCG ssorcg_
00212 #define SSORSI ssorsi_
00213 #define RSCG   rscg_
00214 #define RSSI   rssi_
00215 
00216 #else
00217 
00218 extern "C" int dfault_(int *iparm, double *rparm);
00219 
00220 extern "C" int vfill_(int *n, double *u, double *val);
00221 
00222 extern "C" int jcg_(int *n, int *ia, int *ja, double *a, double *rhs,
00223                     double *u, int *iwksp, int *nw, double *wksp,
00224                     int *iparm, double *rparm, int *ier);
00225 
00226 extern "C" int jsi_(int *n, int *ia, int *ja, double *a, double *rhs,
00227                     double *u, int *iwksp, int *nw, double *wksp,
00228                     int *iparm, double *rparm, int *ier);
00229 
00230 extern "C" int sor_(int *n, int *ia, int *ja, double *a, double *rhs,
00231                     double *u, int *iwksp, int *nw, double *wksp,
00232                     int *iparm, double *rparm, int *ier);
00233 
00234 extern "C" int ssorcg_(int *n, int *ia, int *ja, double *a, double *rhs,
00235                        double *u, int *iwksp, int *nw, double *wksp,
00236                        int *iparm, double *rparm, int *ier);
00237 
00238 extern "C" int ssorsi_(int *n, int *ia, int *ja, double *a, double *rhs,
00239                        double *u, int *iwksp, int *nw, double *wksp,
00240                        int *iparm, double *rparm, int *ier);
00241 
00242 extern "C" int rscg_(int *n, int *ia, int *ja, double *a, double *rhs,
00243                      double *u, int *iwksp, int *nw, double *wksp,
00244                      int *iparm, double *rparm, int *ier);
00245 
00246 extern "C" int rssi_(int *n, int *ia, int *ja, double *a, double *rhs,
00247                      double *u, int *iwksp, int *nw, double *wksp,
00248                      int *iparm, double *rparm, int *ier);
00249 
00250 #endif
00251 
00252 int
00253 ItpackLinSolver::solve(void)
00254 {
00255   // Let ITPACK fill in default parameter values
00256   dfault_(iparm, rparm);
00257 
00258   // Override defaults for "textbook" methods
00259   switch (method) {
00260   case ItpackJ:
00261     iparm[5] = 0; iparm[6] = 2; break;
00262   case ItpackGS:
00263     iparm[5] = 0; break;
00264   case ItpackSORFixed:
00265     iparm[5] = 0; rparm[4] = omega; break;
00266   case ItpackSSORFixed:
00267     iparm[5] = 0; rparm[4] = omega; break;
00268   case ItpackRS:
00269     iparm[5] = 0; break;
00270   default:
00271     break;
00272   }
00273 
00274   // Overwrite default max number of iterations
00275   iparm[0] = maxIter;
00276 
00277   // Sparse matrix storage scheme (0 = symmetric, 1 = nonsymmetric)
00278   iparm[4] = 1;
00279 
00280   double *aPtr = theSOE->A;
00281   double *xPtr = theSOE->X;
00282   double *bPtr = theSOE->B;
00283 
00284   int *iaPtr = IA;
00285   int *jaPtr = JA;
00286 
00287   // Number of non-zero entries in matrix
00288   int nnz = iaPtr[n]-1;
00289 
00290   // Copy original ordering of column indices from the SOE
00291   // because ITPACK will reorder the sparse matrix representation
00292   if (theSOE->Aformed == false) {
00293     int *soeColA = theSOE->colA;
00294     for (int i = 0; i < nnz; i++) {
00295       jaPtr[i] = soeColA[i] + 1;  // Add one for FORTRAN indexing
00296     }
00297   }
00298 
00299   int ier = 0;
00300 
00301   // Fill the x vector with zeros as initial guess to solution of Ax=b
00302   //double val = 0.0;
00303   //vfill_(&n, xPtr, &val);
00304 
00305   switch (method) {
00306   case ItpackJCG:
00307     jcg_(&n, iaPtr, jaPtr, aPtr, bPtr, xPtr,
00308          iwksp, &nwksp, wksp, iparm, rparm, &ier);
00309     break;
00310   case ItpackJSI: case ItpackJ:
00311     jsi_(&n, iaPtr, jaPtr, aPtr, bPtr, xPtr,
00312          iwksp, &nwksp, wksp, iparm, rparm, &ier);
00313     break;
00314   case ItpackSOR: case ItpackGS: case ItpackSORFixed:
00315     sor_(&n, iaPtr, jaPtr, aPtr, bPtr, xPtr,
00316          iwksp, &nwksp, wksp, iparm, rparm, &ier);
00317     break;
00318   case ItpackSSORCG:
00319     ssorcg_(&n, iaPtr, jaPtr, aPtr, bPtr, xPtr,
00320             iwksp, &nwksp, wksp, iparm, rparm, &ier);
00321     break;
00322   case ItpackSSORSI: case ItpackSSORFixed:
00323     ssorsi_(&n, iaPtr, jaPtr, aPtr, bPtr, xPtr,
00324             iwksp, &nwksp, wksp, iparm, rparm, &ier);
00325     break;
00326   case ItpackRSCG:
00327     rscg_(&n, iaPtr, jaPtr, aPtr, bPtr, xPtr,
00328           iwksp, &nwksp, wksp, iparm, rparm, &ier);
00329     break;
00330   case ItpackRSSI: case ItpackRS:
00331     rssi_(&n, iaPtr, jaPtr, aPtr, bPtr, xPtr,
00332           iwksp, &nwksp, wksp, iparm, rparm, &ier);
00333     break;
00334   default:
00335     g3ErrorHandler->fatal("%s -- unknown method type %d",
00336                           "ItpackLinSolver::solve()", method);
00337     break;    
00338   }
00339 
00340   // Tell the SOE that solve() has been called
00341   theSOE->Aformed = true;
00342 
00343   if (ier > 0) {
00344     opserr << "ItpackLinSolver::solve() -- returned ier = " << ier << endln;
00345     return -ier;
00346   }
00347   else
00348     return 0;
00349 }
Generated on Mon Oct 23 15:05:28 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation