OpenSees » SRC » system_of_eqn » linearSOE » sparseGEN
SuperLU.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.6 $
00022 // $Date: 2005/04/01 19:56:12 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/sparseGEN/SuperLU.cpp,v $
00024                                                                         
00025                                                                         
00026 //
00027 // Written: fmk 
00028 // Created: Tue 11/96
00029 // Revision: A
00030 //
00031 // Description: This file contains the implementation of SuperLU
00032 //
00033 // What: "@(#) SuperLU.h, revA"
00034 
00035 #include <SuperLU.h>
00036 #include <SparseGenColLinSOE.h>
00037 #include <f2c.h>
00038 #include <math.h>
00039 #include <Channel.h>
00040 #include <FEM_ObjectBroker.h>
00041 
00042 
00043 SuperLU::SuperLU(int perm, 
00044                  double drop_tolerance, 
00045                  int panel, 
00046                  int relx, 
00047                  char symm)
00048 :SparseGenColLinSolver(SOLVER_TAGS_SuperLU),
00049  perm_r(0),perm_c(0), etree(0), sizePerm(0),
00050  relax(relx), permSpec(perm), panelSize(panel), 
00051  drop_tol(drop_tolerance), symmetric(symm)
00052 {
00053   // set_default_options(&options);
00054   options.Fact = DOFACT;
00055   options.Equil = YES;
00056   options.ColPerm = COLAMD;
00057   options.DiagPivotThresh = 1.0;
00058   options.Trans = NOTRANS;
00059   options.IterRefine = NOREFINE;
00060   options.SymmetricMode = NO;
00061   options.PivotGrowth = NO;
00062   options.ConditionNumber = NO;
00063   options.PrintStat = YES;
00064 
00065   if (symmetric == 'Y')
00066     options.SymmetricMode = YES;
00067   L.ncol = 0;
00068   U.ncol = 0;
00069   A.ncol = 0;
00070   B.ncol = 0;
00071 }
00072 
00073 
00074 SuperLU::~SuperLU()
00075 {
00076 
00077   if (perm_r != 0)
00078     delete [] perm_r;
00079   if (perm_c != 0)
00080     delete [] perm_c;
00081   if (etree != 0) {
00082     delete [] etree;
00083     StatFree(&stat);
00084   }
00085 
00086   if (L.ncol != 0)
00087     Destroy_SuperNode_Matrix(&L);
00088   if (U.ncol != 0)
00089     Destroy_CompCol_Matrix(&U);
00090   if (AC.ncol != 0) {
00091     NCPformat *ACstore = (NCPformat *)AC.Store;
00092     SUPERLU_FREE(ACstore->colbeg);
00093     SUPERLU_FREE(ACstore->colend);
00094     SUPERLU_FREE(ACstore);
00095   }
00096   if (A.ncol != 0) {
00097     SUPERLU_FREE(A.Store);
00098   }
00099   if (B.ncol != 0) {
00100     SUPERLU_FREE(B.Store);
00101   }
00102 }
00103 
00104 /*
00105 extern "C" void  dgstrf (char *refact, SuperMatrix* A, double pivThresh, 
00106                          double dropTol, int relax, 
00107                          int panelSize, int *etree,
00108                          void *work, int lwork, int *perm_r, int *perm_c, 
00109                          SuperMatrix *L, SuperMatrix *U, int *info);
00110 
00111 extern "C" void  dgstrs (char *trans, SuperMatrix *L, SuperMatrix *U, 
00112                          int *perm_r, int *perm_c, 
00113                          SuperMatrix *B, int *info);    
00114 
00115 extern "C" void    StatInit(int, int);
00116 
00117 extern "C" void    dCreate_CompCol_Matrix(SuperMatrix *, int, int, int, double *,
00118                                           int *, int *, Stype_t, Dtype_t, Mtype_t);
00119 
00120 
00121 extern "C" void    get_perm_c(int, SuperMatrix *, int *);
00122 
00123 extern "C" void    sp_preorder (char*, SuperMatrix*, int*, int*, SuperMatrix*);
00124 
00125 extern "C" void    dCreate_Dense_Matrix(SuperMatrix *, int, int, double *, int,
00126                                      Stype_t, Dtype_t, Mtype_t);
00127 
00128 */
00129 
00130 int
00131 SuperLU::solve(void)
00132 {
00133     if (theSOE == 0) {
00134         opserr << "WARNING SuperLU::solve(void)- ";
00135         opserr << " No LinearSOE object has been set\n";
00136         return -1;
00137     }
00138     
00139     int n = theSOE->size;
00140     
00141     // check for quick return
00142     if (n == 0)
00143         return 0;
00144 
00145     if (sizePerm == 0) {
00146         opserr << "WARNING SuperLU::solve(void)- ";
00147         opserr << " size for row and col permutations 0 - has setSize() been called?\n";
00148         return -1;
00149     }
00150 
00151     // first copy B into X
00152     double *Xptr = theSOE->X;
00153     double *Bptr = theSOE->B;
00154     for (int i=0; i<n; i++)
00155         *(Xptr++) = *(Bptr++);
00156 
00157     if (theSOE->factored == false) {
00158         // factor the matrix
00159         int info;
00160 
00161         if (L.ncol != 0 && symmetric == 'N') {
00162           Destroy_SuperNode_Matrix(&L);
00163           Destroy_CompCol_Matrix(&U);     
00164         }
00165 
00166         dgstrf(&options, &AC, drop_tol, relax, panelSize,
00167                etree, NULL, 0, perm_c, perm_r, &L, &U, &stat, &info);
00168 
00169 
00170         if (info != 0) {        
00171           opserr << "WARNING SuperLU::solve(void)- ";
00172           opserr << " Error " << info << " returned in factorization dgstrf()\n";
00173           return -info;
00174         }
00175 
00176         if (symmetric == 'Y')
00177           options.Fact= SamePattern_SameRowPerm;
00178         else
00179           options.Fact = SamePattern;
00180         
00181         theSOE->factored = true;
00182     }   
00183 
00184     // do forward and backward substitution
00185     trans_t trans = NOTRANS;
00186     int info;
00187     dgstrs (trans, &L, &U, perm_c, perm_r, &B, &stat, &info);    
00188 
00189     if (info != 0) {    
00190        opserr << "WARNING SuperLU::solve(void)- ";
00191        opserr << " Error " << info << " returned in substitution dgstrs()\n";
00192        return -info;
00193     }
00194 
00195     return 0;
00196 }
00197 
00198 
00199 
00200 
00201 int
00202 SuperLU::setSize()
00203 {
00204     int n = theSOE->size;
00205     if (n > 0) {
00206 
00207       // create space for the permutation vectors 
00208       // and the elimination tree
00209       if (sizePerm < n) {
00210 
00211         if (perm_r != 0)
00212           delete [] perm_r;
00213         perm_r = new int[n];            
00214 
00215         if (perm_c != 0)
00216           delete [] perm_c;
00217         perm_c = new int[n];            
00218 
00219         if (etree != 0)
00220           delete [] etree;
00221         etree = new int[n];             
00222 
00223         if (perm_r == 0 || perm_c == 0 || etree == 0) {
00224           opserr << "WARNING SuperLU::setSize()";
00225           opserr << " - ran out of memory\n";
00226           sizePerm = 0;
00227           return -1;
00228         }               
00229         sizePerm = n;
00230       }
00231 
00232       // initialisation
00233       StatInit(&stat);
00234 
00235       // create the SuperMatrix A       
00236       dCreate_CompCol_Matrix(&A, n, n, theSOE->nnz, theSOE->A, 
00237                              theSOE->rowA, theSOE->colStartA, 
00238                              SLU_NC, SLU_D, SLU_GE);
00239 
00240       // obtain and apply column permutation to give SuperMatrix AC
00241       get_perm_c(permSpec, &A, perm_c);
00242 
00243       sp_preorder(&options, &A, perm_c, etree, &AC);
00244 
00245       // create the rhs SuperMatrix B 
00246       dCreate_Dense_Matrix(&B, n, 1, theSOE->X, n, SLU_DN, SLU_D, SLU_GE);
00247         
00248       // set the refact variable to 'N' after first factorization with new size 
00249       // can set to 'Y'.
00250       options.Fact = DOFACT;
00251 
00252       if (symmetric == 'Y')
00253         options.SymmetricMode=YES;
00254 
00255     } else if (n == 0)
00256         return 0;
00257     else {
00258         opserr << "WARNING SuperLU::setSize()";
00259         opserr << " - order of system <  0\n";
00260         return -1;      
00261     }
00262         
00263     return 0;
00264 }
00265 
00266 int
00267 SuperLU::sendSelf(int cTag, Channel &theChannel)
00268 {
00269     // nothing to do
00270     return 0;
00271 }
00272 
00273 int
00274 SuperLU::recvSelf(int ctag,
00275                   Channel &theChannel, 
00276                   FEM_ObjectBroker &theBroker)
00277 {
00278     // nothing to do
00279     return 0;
00280 }
00281 
00282 
00283 
00284 
00285 
00286 
00287 
00288 
00289 
00290 
00291 
00292 
00293 
00294 
00295
Generated on Mon Oct 23 15:05:29 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation