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/ProfileSPDLinDirectSkypackSolver.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/system_of_eqn/linearSOE/ProfileSPD/ProfileSPDLinDirectSkypackSolver.C
00027 //
00028 // Written: fmk 
00029 // Created: 03/98
00030 // Revision: A
00031 //
00032 // Description: This file contains the class definition for 
00033 // ProfileSPDLinDirectSkypackSolver. ProfileSPDLinDirectSkypackSolver 
00034 // is a subclass of LinearSOESOlver. It solves a ProfileSPDLinSOE object using
00035 // the Skypack library developed by Osni Marques, software available at
00036 //   http://www.nersc.gov/research/SCG/Osni/marques_software.html
00037 
00038 // What: "@(#) ProfileSPDLinDirectSkypackSolver.C, revA"
00039 
00040 #include <ProfileSPDLinDirectSkypackSolver.h>
00041 #include <ProfileSPDLinSOE.h>
00042 #include <f2c.h>
00043 #include <math.h>
00044 
00045 #include <math.h>
00046 #include <Channel.h>
00047 #include <FEM_ObjectBroker.h>
00048 
00049 #include <Timer.h>
00050 
00051 ProfileSPDLinDirectSkypackSolver::ProfileSPDLinDirectSkypackSolver()
00052 :ProfileSPDLinSolver(SOLVER_TAGS_ProfileSPDLinDirectSkypackSolver),
00053  mCols(0), mRows(0),rw(0),tw(0), index(0),
00054  size(0), invD(0)
00055 {
00056 
00057 }
00058 
00059 ProfileSPDLinDirectSkypackSolver::ProfileSPDLinDirectSkypackSolver(int Mcols, int Mrows)
00060 :ProfileSPDLinSolver(SOLVER_TAGS_ProfileSPDLinDirectSkypackSolver),
00061  mCols(Mcols), mRows(Mrows),rw(0),tw(0), index(0),
00062  size(0), invD(0)
00063 {
00064     if (mCols != 0 && mRows != 0) {
00065  rw = new double[mRows*mCols];
00066  tw = new double[mRows*mRows];
00067  if (mCols > mRows)
00068      index = new int[mCols];
00069  else
00070      index = new int[mRows];
00071     } else { // make sure mCols and mRows == 0
00072  mCols = 0;
00073  mRows = 0;
00074     }
00075 
00076     // check we got the space requested
00077     if (rw == 0 || tw == 0 || index == 0) {
00078  cerr << "WARNING ProfileSPDLinDirectSkypackSolver::ProfileSPDLinDirectSkypack";
00079  cerr << "Solver() - ran out of memory for work areas, setting mCols and mRows = 0\n";
00080  mCols = 0; mRows = 0;
00081     }
00082 }
00083 
00084     
00085 ProfileSPDLinDirectSkypackSolver::~ProfileSPDLinDirectSkypackSolver()
00086 {
00087  if (invD != 0) delete [] invD;
00088  if (rw != 0) delete [] rw;
00089  if (tw != 0) delete [] tw;
00090  if (index != 0) delete [] index;
00091 }
00092 
00093 int
00094 ProfileSPDLinDirectSkypackSolver::setSize(void)
00095 {
00096     int result = 0;
00097 
00098     if (theSOE == 0) {
00099  cerr << "WARNING ProfileSPDLinDirectSkypackSolver::setSize()";
00100  cerr << " No system has been set\n";
00101  return -1;
00102     }
00103 
00104     // check for quick return 
00105     if (theSOE->size == 0)
00106  return 0;
00107     
00108     size = theSOE->size;
00109     block[0] = 1;
00110     block[1] = size;
00111     block[2] = 1;
00112 
00113     // delete the old work areas
00114     if (invD != 0) 
00115  delete [] invD;
00116 
00117     // create the new work areas
00118     invD = new double[size];
00119     if (invD == 0) {
00120       cerr << "Warning :ProfileSPDLinDirectSkypackSolver::setSize():";
00121       cerr << " ran out of memory for work area invD\n";
00122       result = -2;;
00123     }    
00124 
00125     return result;
00126 }
00127 
00128 extern "C" int skysf2_(double *A, double *D, int *DGPNT, int *JMIN, int *JMAX); 
00129 
00130 extern "C" int skypf2_(int *index, int *JMAX, int *JMIN, int *MCOLS, int *MROWS, int *DGPNT, 
00131          double *A, double *RW, double *TW);
00132 
00133 extern "C" int skyss_(int *LDX, int *N, int *NRHS, 
00134         double *A, double *D, double *X, int *DGPNT,
00135         int *BLOCK, int *NBLOCK,
00136         char *FNAME, int *FUNIT, int *INFO);
00137 
00138 int 
00139 ProfileSPDLinDirectSkypackSolver::solve(void)
00140 {
00141     // check for quick returns
00142     if (theSOE == 0) {
00143  cerr << "ProfileSPDLinDirectSkypackSolver::solve(void): ";
00144  cerr << " - No ProfileSPDSOE has been assigned\n";
00145  return -1;
00146     }
00147     
00148     if (theSOE->size == 0)
00149  return 0;    
00150     
00151 
00152     // check that work area invD has been created
00153     if (invD == 0) {
00154  cerr << "ProfileSPDLinDirectSkypackSolver::solve(void): ";
00155  cerr << " - no space for invD - has setSize() been called?\n";
00156  return -1;
00157     } 
00158 
00159     // set some pointers
00160     double *A = theSOE->A;
00161     double *B = theSOE->B;
00162     double *X = theSOE->X;
00163     int *iDiagLoc = theSOE->iDiagLoc;
00164     int theSize = theSOE->size;
00165     // copy B into X
00166     for (int ii=0; ii<theSize; ii++)
00167  X[ii] = B[ii];
00168 
00169     char *FILE = "INCORE";
00170     
00171     if (theSOE->isAfactored == false)  {
00172       
00173  // FACTOR 
00174  if (mRows == 0 || mCols == 0) { // factor using skysf2_
00175      int JMIN =1;
00176      int JMAX = size;
00177      skysf2_(A, invD, iDiagLoc, &JMIN, &JMAX); 
00178  }
00179  else { // factor using skypf2_
00180      int JMIN =1;
00181      int JMAX = size;
00182      int MCOLS = mCols;
00183      int MROWS = mRows;
00184      double *RW = rw;
00185      double *TW = tw;
00186      int *INDEX = index;
00187      skypf2_(INDEX, &JMAX, &JMIN, &MCOLS, &MROWS, iDiagLoc, A, RW, TW);
00188 
00189      // set up invD
00190      for (int i=0; i<size; i++)
00191   invD[i] = 1.0/A[iDiagLoc[i]-1]; // iDiagLoc has fortran array indexing
00192  }
00193       
00194  // mark the system as having been factored
00195  theSOE->isAfactored = true;
00196  theSOE->numInt = 0;      
00197     }
00198 
00199     /*
00200      * now do the forward and back substitution
00201      */
00202 
00203     int LDX = theSize;
00204     int NRHS = 1;
00205     int numBlock = 1;
00206     int fileFD = 0;
00207     int INFO;
00208     A = theSOE->A;
00209     X = theSOE->X;
00210     int *BLOCK = &block[0];
00211     iDiagLoc = theSOE->iDiagLoc;
00212     
00213     skyss_(&LDX, &theSize, &NRHS, A, invD, X, iDiagLoc, BLOCK, &numBlock, 
00214     FILE,  &fileFD, &INFO);
00215       
00216     // return
00217     if (INFO < 0) {
00218  cerr << "WARNING ProfileSPDLinDirectSkypackSolver::solve()";
00219  cerr << " error value returned from skyss()\n";
00220     }    
00221 
00222     return INFO;
00223 
00224 }
00225 
00226 
00227 int 
00228 ProfileSPDLinDirectSkypackSolver::setProfileSOE(ProfileSPDLinSOE &theNewSOE)
00229 {
00230     theSOE = &theNewSOE;
00231     return 0;
00232 }
00233  
00234 int
00235 ProfileSPDLinDirectSkypackSolver::sendSelf(int cTag,
00236         Channel &theChannel)
00237 {
00238     if (size != 0)
00239  cerr << "ProfileSPDLinDirectSkypackSolver::sendSelf - does not send itself YET\n"; 
00240     return 0;
00241 }
00242 
00243 
00244 int 
00245 ProfileSPDLinDirectSkypackSolver::recvSelf(int cTag,
00246         Channel &theChannel, 
00247         FEM_ObjectBroker &theBroker)
00248 {
00249     return 0;
00250 }
00251 
00252