Amesos Package Browser (Single Doxygen Collection) Development
SpoolesOO.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: Direct Sparse Solver Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "SpoolesOO.h"
00030 #ifdef EPETRA_MPI
00031 #include "Epetra_MpiComm.h"
00032 #else
00033 #include "Epetra_Comm.h"
00034 #endif
00035 #include "Epetra_Map.h"
00036 #include "Epetra_Vector.h"
00037 #include "Epetra_RowMatrix.h"
00038 #include "Epetra_Operator.h"
00039 #include "Epetra_Import.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #include <sstream>
00042 #include <vector>
00043 
00044 //
00045 //  SPOOLES include file:
00046 //
00047 extern "C" {
00048 #include "BridgeMPI.h"
00049 }
00050 
00051 //=============================================================================
00052 SpoolesOO::SpoolesOO(Epetra_RowMatrix * A, 
00053      Epetra_MultiVector * X,
00054      Epetra_MultiVector * B) {
00055   //  AllocAzArrays();
00056   SetSpoolesDefaults();
00057 
00058   inConstructor_ = true;  // Shut down complaints about zero pointers for a while
00059   SetUserMatrix(A);
00060   
00061   SetLHS(X);
00062   SetRHS(B);
00063   inConstructor_ = false;
00064 }
00065 
00066 //=============================================================================
00067 SpoolesOO::SpoolesOO() {
00068   //  AllocAzArrays();
00069   SetSpoolesDefaults();
00070 }
00071 
00072 //=============================================================================
00073 SpoolesOO::~SpoolesOO(void) {
00074   //  DeleteMemory();
00075   //  DeleteAzArrays();
00076 }
00077 
00078 //=============================================================================
00079 int SpoolesOO::SetUserMatrix(Epetra_RowMatrix * UserMatrix) {
00080 
00081   if (UserMatrix == 0 && inConstructor_ == true) return(0);
00082   if (UserMatrix == 0) EPETRA_CHK_ERR(-1);
00083 
00084   UserMatrix_ = UserMatrix;
00085 
00086   return(0);
00087 }
00088 
00089 //=============================================================================
00090 int SpoolesOO::SetLHS(Epetra_MultiVector * X) {
00091 
00092   if (X == 0 && inConstructor_ == true) return(0);
00093   if (X == 0) EPETRA_CHK_ERR(-1);
00094   X_ = X;
00095   X_->ExtractView(&x_, &x_LDA_);
00096   return(0);
00097 }
00098 //=============================================================================
00099 int SpoolesOO::SetRHS(Epetra_MultiVector * B) {
00100 
00101   if (B == 0 && inConstructor_ == true) return(0);
00102   if (B == 0) EPETRA_CHK_ERR(-1);
00103   B_ = B;
00104   B_->ExtractView(&b_, &b_LDA_);
00105 
00106   return(0);
00107 }
00108 int SpoolesOO::SetSpoolesDefaults() {
00109 
00110  UserOperator_ = 0;
00111  UserMatrix_ = 0;
00112  // PrecOperator_ = 0;
00113  // PrecMatrix_ = 0;
00114  X_ = 0;
00115  B_ = 0;
00116  
00117  x_LDA_ = 0;
00118  x_ = 0;
00119  b_LDA_ = 0;
00120  b_ = 0;
00121 
00122  return(0);
00123 
00124 }
00125 
00126 //=============================================================================
00127 
00128 int SpoolesOO::Solve() { 
00129   FILE *msgFile = fopen("msgFile", "w" ) ; 
00130   FILE *matFile = 0 ;
00131 
00132   Epetra_RowMatrix   *RowMatrixA = (GetUserMatrix()) ; 
00133   const Epetra_Comm &Comm = RowMatrixA->Comm();
00134   int iam =  Comm.MyPID() ;
00135   bool verbose = false ; // Other option is (iam == 0 ) ; 
00136   verbose = ( iam == 0 ) ; 
00137   if ( verbose ) matFile = fopen("SPO.m", "w" ) ; 
00138 
00139   Epetra_CrsMatrix   *matA = dynamic_cast<Epetra_CrsMatrix*>(RowMatrixA) ; 
00140   assert( matA != NULL ) ; // SuperludistOO shows how to convert a 
00141   // non-CrsMatrix into a CrsMatrix
00142 
00143   Epetra_MultiVector   *vecX = GetLHS() ; 
00144   Epetra_MultiVector   *vecB = GetRHS() ; 
00145 
00146   const Epetra_Map &matAmap = matA->RowMap() ; 
00147   const Epetra_MpiComm & comm1 = dynamic_cast<const Epetra_MpiComm &> (Comm);
00148   MPI_Comm MPIC = comm1.Comm() ;
00149 
00150   int IsLocal = ( matAmap.NumMyElements() == 
00151       matAmap.NumGlobalElements() )?1:0;
00152   Comm.Broadcast( &IsLocal, 1, 0 ) ; 
00153 
00154   EPETRA_CHK_ERR( 1 - IsLocal  ); // SuperludistOO shows hows to 
00155   // deal with distributed matrices.  
00156 
00157   int nArows = matA->NumGlobalRows() ; 
00158   int nAcols = matA->NumGlobalCols() ; 
00159   int nproc =  Comm.NumProc() ;
00160   
00161   assert( vecX->NumVectors() == 1 ) ; 
00162   assert( vecB->NumVectors() == 1 ) ; 
00163 
00164   //  assert( nArows == vecX->MyLength() ) ; 
00165   //  assert( nAcols == vecB->MyLength() ) ;
00166 
00167   //  assert( matAmap.DistributedGlobal() == false ) ; 
00168   if ( iam == 0 ) { 
00169     assert( matAmap.NumMyElements() == matAmap.NumGlobalElements() ) ;
00170   } else { 
00171     assert( matAmap.NumMyElements() == 0 ) ;
00172   } 
00173 
00174   int numentries = matA->NumGlobalNonzeros();
00175   vector <int>rowindices( numentries ) ; 
00176   vector <int>colindices( numentries ) ; 
00177   vector <double>values( numentries ) ; 
00178   int NumEntries;
00179   double *RowValues;
00180 
00181   int *ColIndices;
00182   int numrows = matA->NumGlobalRows();
00183   int entrynum = 0 ; 
00184 
00185   //
00186   //  The following lines allow us time to attach the debugger
00187   //
00188   char hatever;
00189   //  if ( iam == 0 )  cin >> hatever ; 
00190   Comm.Barrier();
00191 
00192 
00193   InpMtx *mtxA; 
00194   if ( iam==0 ) {
00195     //
00196     // Create mtxA on proc 0 
00197     //
00198     for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00199       assert( matA->ExtractMyRowView( MyRow, NumEntries, RowValues, ColIndices ) == 0 ) ;
00200       for ( int j = 0; j < NumEntries; j++ ) { 
00201   rowindices[entrynum] = MyRow ; 
00202   colindices[entrynum] = ColIndices[j] ; 
00203   values[entrynum] = RowValues[j] ; 
00204   entrynum++; 
00205       }
00206     }
00207     mtxA = InpMtx_new() ; 
00208     
00209     InpMtx_init( mtxA, 1, SPOOLES_REAL, 0, 0 ) ;   // I can't figure out 
00210     //  what the right bound for the number of vectors is - Ken 
00211     
00212     if ( GetTrans() ) { 
00213       InpMtx_inputRealTriples( mtxA, numentries, &colindices[0], 
00214              &rowindices[0], &values[0] ) ; 
00215     } else { 
00216       InpMtx_inputRealTriples( mtxA, numentries, &rowindices[0], 
00217              &colindices[0], &values[0] ) ; 
00218     }
00219   } else {
00220     mtxA = 0 ; 
00221   }
00222 
00223   //  return 1;    OK to here
00224   
00225   DenseMtx *mtxX = 0 ; 
00226   DenseMtx *mtxY = 0 ;
00227   double *bValues ;
00228   double *xValues ;
00229   if ( iam == 0 ) { 
00230     //
00231     //  Convert Epetra_Vector x and Epetra_Vector b arrays
00232     //
00233     int bLda, xLda ; 
00234     
00235     assert( vecB->ExtractView( &bValues, &bLda ) == 0 )  ; 
00236     assert( vecX->ExtractView( &xValues, &xLda ) == 0 ) ; 
00237   
00238   
00239     //
00240     //  SPOOLES matrices
00241     //
00242     mtxX = DenseMtx_new() ; 
00243     mtxY = DenseMtx_new() ; 
00244 #ifdef OLDWAY
00245     DenseMtx_init( mtxX, SPOOLES_REAL, -1, -1, numrows, 1, 1, numrows );
00246     DenseMtx_init( mtxY, SPOOLES_REAL, -1, -1, numrows, 1, 1, numrows );
00247 #else
00248     DenseMtx_init( mtxX, SPOOLES_REAL, 0, 0, numrows, 1, 1, numrows );
00249     DenseMtx_init( mtxY, SPOOLES_REAL, 0, 0, numrows, 1, 1, numrows );
00250 #endif
00251     int nrow, nnrhs ; 
00252     DenseMtx_dimensions( mtxY, &nrow, &nnrhs ) ;
00253     assert( nrow = numrows ) ; 
00254     assert( nnrhs == 1 ) ; 
00255     
00256   
00257    if (verbose) InpMtx_writeForMatlab( mtxA, "mtxA", matFile ) ; 
00258     //
00259     //  This is a maximally inefficient way to create the right hand side
00260     //  matrix, but I could not find anything better in the documentation.
00261     //
00262     for ( int i = 0 ; i < numrows; i ++ ) 
00263       {
00264   DenseMtx_setRealEntry( mtxY, i, 0, bValues[i] );
00265       }
00266     if ( verbose ) DenseMtx_writeForMatlab( mtxY, "mtxY", matFile ) ; 
00267     DenseMtx_zero(mtxX) ;
00268   }
00269   
00270   
00271   int type = SPOOLES_REAL ;
00272   int symmetryflag = SPOOLES_NONSYMMETRIC ;
00273   // SPOOLES requires a message level and a message File 
00274   int msglvl = 0;        
00275   int rc;                // return code 
00276   
00277   /*
00278     --------------------------------
00279     create and setup a Bridge object
00280     --------------------------------
00281   */
00282   BridgeMPI     *bridgeMPI = BridgeMPI_new() ;
00283   BridgeMPI_setMPIparams(bridgeMPI, nproc, iam, MPIC ) ;
00284   BridgeMPI_setMatrixParams(bridgeMPI, numrows, type, symmetryflag) ;
00285   double tau = 100.0 ; 
00286   double droptol = 0.0 ; 
00287   int lookahead = 0 ; 
00288   PatchAndGoInfo *PatchAndGo = 0 ; 
00289   BridgeMPI_setFactorParams(bridgeMPI, FRONTMTX_DENSE_FRONTS, 
00290           SPOOLES_PIVOTING, tau, droptol, lookahead, 
00291           PatchAndGo ) ; 
00292   BridgeMPI_setMessageInfo(bridgeMPI, msglvl, msgFile) ;
00293   assert( type == SPOOLES_REAL ) ; 
00294   assert( msglvl >= 0 && msglvl <= 2 ) ; 
00295     
00296   rc = BridgeMPI_setup(bridgeMPI, mtxA) ;
00297   if ( rc != 1 ) {
00298     fprintf(stderr, "\n error return %d from BridgeMPI_setup()", rc) ;
00299     exit(-1) ;
00300   }
00301   Comm.Barrier();
00302   
00303   int nfront, nfind, nfent, nsolveops;
00304   double nfactorops;
00305   rc = BridgeMPI_factorStats(bridgeMPI, type, symmetryflag, &nfront,
00306            &nfind, &nfent, &nsolveops, &nfactorops) ;
00307   if ( rc != 1 ) {
00308     fprintf(stderr,
00309       "\n error return %d from BridgeMPI_factorStats()", rc) ;
00310     exit(-1) ;
00311   }
00312   
00313   /*
00314     --------------------------------
00315     setup the parallel factorization
00316     --------------------------------
00317   */
00318   rc = BridgeMPI_factorSetup(bridgeMPI, 0, 0.0) ;
00319   if (rc != 1 ) {
00320     std::stringstream Message ; 
00321     
00322     Message <<  " SPOOLES factorsetup failed with return code " << 
00323       rc ;
00324     string mess = Message.str() ; 
00325     throw( mess ) ; 
00326   }
00327   
00328   /*
00329     -----------------
00330     factor the matrix
00331     -----------------
00332   */
00333   int permuteflag = 1 ;
00334   int errorcode ;
00335   rc = BridgeMPI_factor(bridgeMPI, mtxA, permuteflag, &errorcode) ;
00336   assert( permuteflag == 1 ) ; 
00337   
00338   if (rc != 1 ) {
00339     std::stringstream Message ; 
00340     
00341     Message <<  " SPOOLES factorization failed with return code " << 
00342       rc << " and error code " << errorcode  ; 
00343     string mess = Message.str() ; 
00344     throw( mess ) ; 
00345   }
00346   /*
00347     ----------------
00348     solve the system
00349     ----------------
00350   */
00351   
00352   rc = BridgeMPI_solveSetup(bridgeMPI) ;
00353   if (rc != 1 ) {
00354     std::stringstream Message ; 
00355     
00356     Message <<  " SPOOLES BridgeMPI_solveSetup failed with return code" << 
00357       rc << endl ; 
00358     string mess = Message.str() ; 
00359     throw( mess ) ; 
00360   }
00361 
00362   assert( permuteflag == 1 ) ; 
00363   rc = BridgeMPI_solve(bridgeMPI, permuteflag, mtxX, mtxY) ;
00364   assert( permuteflag == 1 ) ; 
00365   if (rc != 1 ) {
00366     std::stringstream Message ; 
00367     
00368     Message <<  " SPOOLES BridgeMPI_solve failed with return code" << 
00369       rc << endl ; 
00370     string mess = Message.str() ; 
00371     throw( mess ) ; 
00372   }
00373   
00374   if ( verbose ) DenseMtx_writeForMatlab( mtxX, "mtxXX", matFile ) ; 
00375   if ( verbose ) fclose(matFile);
00376 
00377   //  Result->SolveTime().Time_First( ) ;
00378   //
00379   //  Here we copy the values of B back from mtxX into xValues
00380   //
00381   if ( iam == 0 ) { 
00382     for ( int i = 0 ; i < numrows; i ++ ) 
00383       {
00384   DenseMtx_realEntry( mtxX, i, 0, &xValues[i] );
00385       }
00386     //  DenseMtx_writeForMatlab( mtxX, "mtxX", matFile ) ; 
00387   }
00388   
00389 
00390   if ( iam == 0 ) {
00391     InpMtx_free(mtxA) ;
00392     DenseMtx_free(mtxX) ;
00393     DenseMtx_free(mtxY) ;
00394   }
00395   BridgeMPI_free(bridgeMPI) ;
00396 
00397   Comm.Barrier();
00398 
00399   return(0) ; 
00400 }
00401 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines