My title

00001 // File: ~/system_of_eqn/linearSOE/symLinSolver/SymSparseLinSolver.C
00002 //
00003 // Written: Jun Peng
00004 // Created: 12/98
00005 // Revision: A
00006 //
00007 // Description: This file contains the class definition for 
00008 // SymSparseinSolver. It solves the SymSparseLinSOE object by calling
00009 // some "C" functions. The solver used here is generalized sparse
00010 // solver. The user can choose three different ordering schema.
00011 //
00012 // What: "@(#) SymSparseLinSolver.C, revA"
00013 
00014 
00015 #include "SymSparseLinSOE.h"
00016 #include "SymSparseLinSolver.h"
00017 #include <f2c.h>
00018 #include <math.h>
00019 #include <Channel.h>
00020 #include <FEM_ObjectBroker.h>
00021 
00022 extern "C" {
00023   #include "FeStructs.h"
00024   #include "globalVars.h"
00025   #include "tim.h"
00026 }
00027 
00028 
00029 SymSparseLinSolver::SymSparseLinSolver()
00030 :LinearSOESolver(SOLVER_TAGS_SymSparseLinSolver),
00031  theSOE(0), factored(false)
00032 {
00033     // nothing to do.
00034 }
00035 
00036 
00037 SymSparseLinSolver::~SymSparseLinSolver()
00038 { 
00039     // nothing to do.
00040 }
00041 
00042 
00043 extern "C" int pfsfct(int neqns, double *diag, double **penv, int nblks, int *xblk,
00044         OFFDBLK **begblk, OFFDBLK *first, int *rowblks);
00045 
00046 extern "C" void pfsslv(int neqns, double *diag, double **penv, int nblks,
00047          int *xblk, double *rhs, OFFDBLK **begblk);
00048 
00049 
00050 int
00051 SymSparseLinSolver::solve(void)
00052 { 
00053     if (theSOE == 0) {
00054  cerr << "WARNING SymSparseLinSolver::solve(void)- ";
00055  cerr << " No LinearSOE object has been set\n";
00056  return -1;
00057     }
00058 
00059     int      nblks = theSOE->nblks;
00060     int      *xblk = theSOE->xblk;
00061     int      *invp = theSOE->invp;
00062     double   *diag = theSOE->diag;
00063     double   **penv = theSOE->penv;
00064     int      *rowblks = theSOE->rowblks;
00065     OFFDBLK  **begblk = theSOE->begblk;
00066     OFFDBLK  *first = theSOE->first;
00067     
00068     
00069     int neq = theSOE->size;
00070 
00071     // check for quick return
00072     if (neq == 0)
00073  return 0;
00074 
00075     // first copy B into X
00076 
00077     for (int i=0; i<neq; i++) {
00078         theSOE->X[i] = theSOE->B[i];
00079     }
00080     double *Xptr = theSOE->X;
00081 
00082     timer (FACTOR);
00083     if (factored == false) {
00084 
00085         //factor the matrix
00086         //call the "C" function to do the numerical factorization.
00087         int factor;
00088  factor = pfsfct(neq, diag, penv, nblks, xblk, begblk, first, rowblks);
00089  if (factor>0) {
00090      cerr << "In SymSparseLinSolver: error in factorization.\n";
00091      return -1;
00092  }
00093  factored = true;
00094     }
00095 
00096     //do forward and backward substitution.
00097     //call the "C" function.
00098     timer (SOLVE);
00099     pfsslv(neq, diag, penv, nblks, xblk, Xptr, begblk);
00100 
00101 
00102     //Since the X we get by solving AX=B is P*X, we need to reordering
00103     //the Xptr to ge the wanted X.
00104 
00105     double *tempX = new double[neq];
00106     if (tempX == 0) {
00107         cerr << "WARNING LawSparseLinSover::SymSparseLinSolver :";
00108  cerr << " ran out of memory for vectors (tempX) ";
00109  return -1;
00110     } 
00111 
00112     for (int m=0; m<neq; m++) {
00113         tempX[m] = Xptr[invp[m]];
00114     }
00115  
00116     for (int k=0; k<neq; k++) {
00117         Xptr[k] = tempX[k];
00118     }
00119  
00120     delete [] tempX;
00121 
00122     timer (END);
00123 
00124     return 0;
00125 }
00126 
00127 
00128 int
00129 SymSparseLinSolver::setSize()
00130 {
00131     // nothing to do
00132     return 0;
00133 }
00134 
00135 
00136 int
00137 SymSparseLinSolver::setLinearSOE(SymSparseLinSOE &theLinearSOE)
00138 {
00139     theSOE = &theLinearSOE;
00140     return 0;
00141 }
00142 
00143 
00144 int
00145 SymSparseLinSolver::sendSelf(int cTAg, Channel &theChannel)
00146 {
00147     // nothing to do
00148     return 0;
00149 }
00150 
00151 
00152 int
00153 SymSparseLinSolver::recvSelf(int cTag,
00154         Channel &theChannel, FEM_ObjectBroker &theBroker)
00155 {
00156     // nothing to do
00157     return 0;
00158 }
00159 
00160 
00161 
00162
SymSparseLinSolver.cpp
