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/SuperLU.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/system_of_eqn/linearSOE/sparseGEN/SuperLU.h
00027 //
00028 // Written: fmk 
00029 // Created: Tue 11/96
00030 // Revision: A
00031 //
00032 // Description: This file contains the implementation of SuperLU
00033 //
00034 // What: "@(#) SuperLU.h, revA"
00035 
00036 #include <SuperLU.h>
00037 #include <SparseGenColLinSOE.h>
00038 #include <f2c.h>
00039 #include <math.h>
00040 #include <Channel.h>
00041 #include <FEM_ObjectBroker.h>
00042 
00043 
00044 SuperLU::SuperLU(int perm, double thrsh, int relx, int panel)
00045 :SparseGenColLinSolver(SOLVER_TAGS_SuperLU),perm_r(0),perm_c(0),
00046  etree(0), sizePerm(0),
00047  relax(relx), permSpec(perm),  panelSize(panel), thresh(thrsh)
00048 {
00049   refact[0] = 'N';    
00050 }
00051 
00052 
00053 SuperLU::~SuperLU()
00054 {
00055   if (perm_r != 0)
00056  delete [] perm_r;
00057   if (perm_c != 0)
00058  delete [] perm_c;
00059   if (etree != 0)
00060  delete [] etree;
00061 }
00062 
00063 extern "C" void  dgstrf (char *refact, SuperMatrix* A, double pivThresh, 
00064     double dropTol, int relax, 
00065     int panelSize, int *etree,
00066     void *work, int lwork, int *perm_r, int *perm_c, 
00067     SuperMatrix *L, SuperMatrix *U, int *info);
00068 
00069 extern "C" void  dgstrs (char *trans, SuperMatrix *L, SuperMatrix *U, 
00070     int *perm_r, int *perm_c, 
00071     SuperMatrix *B, int *info);    
00072 
00073 
00074 extern "C" void    StatInit(int, int);
00075 
00076 extern "C" void    dCreate_CompCol_Matrix(SuperMatrix *, int, int, int, double *,
00077        int *, int *, Stype_t, Dtype_t, Mtype_t);
00078 
00079 
00080 extern "C" void    get_perm_c(int, SuperMatrix *, int *);
00081 
00082 extern "C" void    sp_preorder (char*, SuperMatrix*, int*, int*, SuperMatrix*);
00083 
00084 extern "C" void    dCreate_Dense_Matrix(SuperMatrix *, int, int, double *, int,
00085                               Stype_t, Dtype_t, Mtype_t);
00086 
00087 
00088 int
00089 SuperLU::solve(void)
00090 {
00091     if (theSOE == 0) {
00092  cerr << "WARNING SuperLU::solve(void)- ";
00093  cerr << " No LinearSOE object has been set\n";
00094  return -1;
00095     }
00096     
00097     int n = theSOE->size;
00098     
00099     // check for quick return
00100     if (n == 0)
00101  return 0;
00102 
00103     if (sizePerm == 0) {
00104  cerr << "WARNING SuperLU::solve(void)- ";
00105  cerr << " size for row and col permutations 0 - has setSize() been called?\n";
00106  return -1;
00107     }
00108 
00109     // first copy B into X
00110     double *Xptr = theSOE->X;
00111     double *Bptr = theSOE->B;
00112     for (int i=0; i<n; i++)
00113  *(Xptr++) = *(Bptr++);
00114 
00115     if (theSOE->factored == false) {
00116  // factor the matrix
00117  int info;
00118  dgstrf(refact, &AC, thresh, 0.0, relax, panelSize, etree, NULL, 0,
00119         perm_r, perm_c, &L, &U, &info);
00120 
00121  if (info != 0) { 
00122     cerr << "WARNING SuperLU::solve(void)- ";
00123     cerr << " Error " << info << " returned in factorization dgstrf()\n";
00124     return -info;
00125  }
00126  refact[0] = 'Y';
00127  theSOE->factored = true;
00128     } 
00129 
00130     // do forward and backward substitution
00131     char trans[1]; trans[0] = 'N';
00132     int info;
00133     dgstrs (trans, &L, &U, perm_r, perm_c, &B, &info);    
00134 
00135     if (info != 0) { 
00136        cerr << "WARNING SuperLU::solve(void)- ";
00137        cerr << " Error " << info << " returned in substitution dgstrs()\n";
00138        return -info;
00139     }
00140 
00141     return 0;
00142 }
00143 
00144 
00145 
00146 
00147 int
00148 SuperLU::setSize()
00149 {
00150     int n = theSOE->size;
00151     if (n > 0) {
00152 
00153       // create space for the permutation vectors 
00154       // and the elimination tree
00155       if (sizePerm < n) {
00156 
00157  if (perm_r != 0)
00158    delete [] perm_r;
00159  perm_r = new int[n];  
00160 
00161  if (perm_c != 0)
00162    delete [] perm_c;
00163  perm_c = new int[n];  
00164 
00165  if (etree != 0)
00166    delete [] etree;
00167  etree = new int[n];  
00168 
00169  if (perm_r == 0 || perm_c == 0 || etree == 0) {
00170    cerr << "WARNING SuperLU::setSize()";
00171    cerr << " - ran out of memory\n";
00172    sizePerm = 0;
00173    return -1;
00174  }  
00175  sizePerm = n;
00176       }
00177 
00178       // initialisation
00179       StatInit(panelSize, relax);
00180 
00181       // create the SuperMatrix A 
00182       dCreate_CompCol_Matrix(&A, n, n, theSOE->nnz, theSOE->A, 
00183         theSOE->rowA, theSOE->colStartA, 
00184         NC, _D, GE);
00185 
00186       // obtain and apply column permutation to give SuperMatrix AC
00187       get_perm_c(permSpec, &A, perm_c);
00188 
00189       sp_preorder(refact, &A, perm_c, etree, &AC);
00190 
00191       // create the rhs SuperMatrix B 
00192       dCreate_Dense_Matrix(&B, n, 1, theSOE->X, n, DN, _D, GE);
00193  
00194       // set the refact variable to 'N' after first factorization with new size 
00195       // can set to 'Y'.
00196       refact[0] = 'N';
00197 
00198     } else if (n == 0)
00199  return 0;
00200     else {
00201  cerr << "WARNING SuperLU::setSize()";
00202  cerr << " - order of system <  0\n";
00203  return -1; 
00204     }
00205  
00206     return 0;
00207 }
00208 
00209 int
00210 SuperLU::sendSelf(int cTag, Channel &theChannel)
00211 {
00212     // nothing to do
00213     return 0;
00214 }
00215 
00216 int
00217 SuperLU::recvSelf(int ctag,
00218     Channel &theChannel, 
00219     FEM_ObjectBroker &theBroker)
00220 {
00221     // nothing to do
00222     return 0;
00223 }
00224 
00225 
00226 
00227 
00228 
00229 
00230 
00231 
00232 
00233 
00234 
00235 
00236 
00237 
00238