OpenSees » SRC » system_of_eqn » linearSOE » mumps

MumpsSolver.cpp

Go to the documentation of this file.
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 // Description: This file contains the implementation of MumpsSolver
00022 
00023 // $Revision: 1.3 $
00024 // $Date: 2006/04/04 22:59:50 $
00025 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/mumps/MumpsSolver.cpp,v $
00026 
00027 // Written: fmk 
00028 // Created: 02/06
00029 
00030 #include <MumpsSolver.h>
00031 #include <MumpsSOE.h>
00032 #include <Channel.h>
00033 #include <FEM_ObjectBroker.h>
00034 #include <OPS_Globals.h>
00035 
00036 #define ICNTL(I) icntl[(I)-1] /* macro s.t. indices match documentation */
00037 
00038 extern "C" {
00039 #include <mpi.h>
00040 }
00041 
00042 MumpsSolver::MumpsSolver()
00043   :LinearSOESolver(SOLVER_TAGS_MumpsSolver),
00044    theMumpsSOE(0)
00045 {
00046   init = false;
00047   id.job=-1; 
00048   id.par=1; 
00049 
00050   id.comm_fortran=MPI_COMM_WORLD;
00051 }
00052 
00053 
00054 MumpsSolver::MumpsSolver(int ICNTL7)
00055   :LinearSOESolver(SOLVER_TAGS_MumpsSolver),
00056    theMumpsSOE(0)
00057 {
00058   init = false;
00059   id.job=-1; 
00060   id.par=1; 
00061 
00062   id.comm_fortran=MPI_COMM_WORLD;
00063   id.ICNTL(7)=ICNTL7;;
00064 }
00065 
00066 MumpsSolver::~MumpsSolver()
00067 {
00068   id.job=-2; 
00069   dmumps_c(&id); /* Terminate instance */
00070 }
00071 
00072 int
00073 MumpsSolver::solve(void)
00074 {
00075   int nnz = theMumpsSOE->nnz;
00076   int n = theMumpsSOE->size;
00077   int *rowA = theMumpsSOE->rowA;
00078   int *colA = theMumpsSOE->colA;
00079 
00080   double *X = theMumpsSOE->X;
00081   double *B = theMumpsSOE->B;
00082 
00083   // increment row and col A values by 1 for mumps fortran indexing
00084   for (int i=0; i<nnz; i++) {
00085     rowA[i]++;
00086     colA[i]++;
00087     //    opserr << rowA[i] << " " << colA[i] << " " << theMumpsSOE->A[i] << endln;
00088   }
00089   
00090   for (int i=0; i<n; i++)
00091     X[i] = B[i];
00092 
00093   int info = 0;
00094   if (theMumpsSOE->factored == false) {
00095 
00096     // factor the matrix
00097     id.n   = theMumpsSOE->size; 
00098     id.nz  = theMumpsSOE->nnz; 
00099     id.irn = theMumpsSOE->rowA;
00100     id.jcn = theMumpsSOE->colA;
00101     id.a   = theMumpsSOE->A; 
00102     id.rhs = theMumpsSOE->X;
00103 
00104     // No outputs 
00105     id.ICNTL(1)=-1; id.ICNTL(2)=-1; id.ICNTL(3)=-1; id.ICNTL(4)=0;
00106     // Call the MUMPS package to factor & solve the system
00107     id.job = 5;
00108     dmumps_c(&id);
00109 
00110     theMumpsSOE->factored = true;
00111   } else {
00112     // factor the matrix
00113     id.n   = theMumpsSOE->size; 
00114     id.nz  = theMumpsSOE->nnz; 
00115     id.irn = theMumpsSOE->rowA;
00116     id.jcn = theMumpsSOE->colA;
00117     id.a   = theMumpsSOE->A; 
00118     id.rhs = theMumpsSOE->X;
00119 
00120     // No outputs 
00121     id.ICNTL(1)=-1; id.ICNTL(2)=-1; id.ICNTL(3)=-1; id.ICNTL(4)=0;
00122     // Call the MUMPS package to factor & solve the system
00123     id.job = 3;
00124     dmumps_c(&id);
00125   }     
00126 
00127 
00128   
00129   info = id.infog[0];
00130   if (info != 0) {      
00131     opserr << "WARNING MumpsSolver::solve(void)- ";
00132     opserr << " Error " << info << " returned in substitution dmumps()\n";
00133     return info;
00134   }
00135 
00136   // decrement row and col A values by 1 to return to C++ indexing
00137   for (int i=0; i<nnz; i++) {
00138     rowA[i]--;
00139     colA[i]--;
00140   }
00141 
00142   return 0;
00143 }
00144 
00145 
00146 int 
00147 MumpsSOE::setMumpsSolver(MumpsSolver &newSolver)
00148 {
00149   return this->setSolver(newSolver);
00150 }
00151 
00152 
00153 int
00154 MumpsSolver::setSize()
00155 {
00156   if (init == false) {
00157     id.sym=theMumpsSOE->matType; 
00158     dmumps_c(&id);
00159     init = true;
00160   }
00161 
00162   int nnz = theMumpsSOE->nnz;
00163   int *rowA = theMumpsSOE->rowA;
00164   int *colA = theMumpsSOE->colA;
00165 
00166   // increment row and col A values by 1 for mumps fortran indexing
00167   for (int i=0; i<nnz; i++) {
00168     rowA[i]++;
00169     colA[i]++;
00170   }
00171 
00172   // analyze the matrix
00173   id.n   = theMumpsSOE->size; 
00174   id.nz  = theMumpsSOE->nnz; 
00175   id.irn = theMumpsSOE->rowA;
00176   id.jcn = theMumpsSOE->colA;
00177   id.a   = theMumpsSOE->A; 
00178   id.rhs = theMumpsSOE->X;
00179   
00180   // No outputs 
00181   id.ICNTL(1)=-1; id.ICNTL(2)=-1; id.ICNTL(3)=-1; id.ICNTL(4)=0;
00182   // Call the MUMPS package to factor & solve the system
00183   id.job = 1;
00184   dmumps_c(&id);
00185 
00186   int info = id.infog[0];
00187   if (info != 0) {      
00188     opserr << "WARNING MumpsSolver::setSize(void)- ";
00189     opserr << " Error " << info << " returned in substitution dmumps()\n";
00190     return info;
00191   }
00192 
00193   // decrement row and col A values by 1 to return to C++ indexing
00194   for (int i=0; i<nnz; i++) {
00195     rowA[i]--;
00196     colA[i]--;
00197   }
00198   return info;
00199 }
00200 
00201 int
00202 MumpsSolver::sendSelf(int cTag, Channel &theChannel)
00203 {
00204   // nothing to do
00205   return 0;
00206 }
00207 
00208 int
00209 MumpsSolver::recvSelf(int ctag,
00210                       Channel &theChannel, 
00211                       FEM_ObjectBroker &theBroker)
00212 {
00213   // nothing to do
00214   return 0;
00215 }
00216 
00217 int 
00218 MumpsSolver::setLinearSOE(MumpsSOE &theSOE)
00219 {
00220   theMumpsSOE = &theSOE;
00221   return 0;
00222 }
00223 
00224 
00225 
00226 
00227 
00228 
00229 
00230 
00231 
00232 
00233 
00234 
00235

Generated on Mon Oct 23 15:05:28 2006 for OpenSees by doxygen 1.5.0