Amesos Package Browser (Single Doxygen Collection) Development
Amesos_Mumps.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 "Amesos_Mumps.h"
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Import.h"
00032 #include "Epetra_RowMatrix.h"
00033 #include "Epetra_CrsMatrix.h"
00034 #include "Epetra_Vector.h"
00035 #include "Epetra_Util.h"
00036 #include "Epetra_Time.h"
00037 #include "Epetra_IntSerialDenseVector.h"
00038 #include "Epetra_SerialDenseVector.h"
00039 #include "Epetra_IntVector.h"
00040 #include "Epetra_Vector.h"
00041 #include "Epetra_SerialDenseMatrix.h"
00042 #include "Epetra_Util.h"
00043 
00044 #define ICNTL(I) icntl[(I)-1]
00045 #define CNTL(I)  cntl[(I)-1]
00046 #define INFOG(I) infog[(I)-1]
00047 #define INFO(I) info[(I)-1]
00048 #define RINFOG(I) rinfog[(I)-1]
00049 
00050 //=============================================================================
00051 Amesos_Mumps::Amesos_Mumps(const Epetra_LinearProblem &prob ) :
00052   IsComputeSchurComplementOK_(false),
00053   NoDestroy_(false),
00054   MaxProcs_(-1),
00055   UseTranspose_(false),
00056   MtxConvTime_(-1), 
00057   MtxRedistTime_(-1), 
00058   VecRedistTime_(-1),
00059   SymFactTime_(-1), 
00060   NumFactTime_(-1), 
00061   SolveTime_(-1),
00062   RowSca_(0),
00063   ColSca_(0),
00064   PermIn_(0),
00065   NumSchurComplementRows_(-1),
00066 #ifdef HAVE_MPI  
00067   MUMPSComm_(0),
00068 #endif
00069   Problem_(&prob)
00070 {
00071   // -777 is for me. It means: I never called MUMPS, so
00072   // SymbolicFactorization should not call Destroy() and ask MUMPS to
00073   // free its space.
00074   MDS.job = -777;
00075   
00076   // load up my default parameters
00077   ICNTL[1]  = -1;  // Turn off error messages  6=on, -1 =off
00078   ICNTL[2]  = -1;  // Turn off diagnostic printing  6=on, -1=off
00079   ICNTL[3]  = -1;  // Turn off global information messages   6=on, -1=off
00080   ICNTL[4]  = -1;  // 3 = most msgs; -1= none  
00081 
00082 #ifdef MUMPS_4_9
00083 
00084   ICNTL[5]  = 0;   // Matrix is given in assembled (i.e. triplet) from
00085   ICNTL[6]  = 7;   // Choose column permutation automatically
00086   ICNTL[7]  = 7;   // Choose ordering method automatically
00087   ICNTL[8]  = 77;  // Choose scaling automatically
00088   ICNTL[9]  = 1;   // Compute A x = b 
00089   ICNTL[10] = 0;   // Maximum steps of iterative refinement
00090   ICNTL[11] = 0;   // Do not collect statistics
00091   ICNTL[12] = 0;   // Automatic choice of ordering strategy
00092   ICNTL[13] = 0;   // Use ScaLAPACK for root node 
00093   ICNTL[14] = 20;  // Increase memory allocation 20% at a time 
00094 
00095   // 15, 16 and 17 are not used
00096   // 18 is set after we know NumProc
00097   // 19 is set after we know Schur status
00098 
00099   ICNTL[20] = 0;   // RHS is given in dense form
00100   ICNTL[21] = 0;   // Solution is assembled at end, not left distributed
00101   ICNTL[22] = 0;   // Do all computations in-core
00102   ICNTL[23] = 0;   // We don't supply maximum MB of working memory
00103   ICNTL[24] = 0;   // Do not perform null pivot detection
00104   ICNTL[25] = 0;   // No null space basis computation, return 1 possible solution
00105   ICNTL[26] = 0;   // Do not condense/reduce Schur RHS
00106   ICNTL[27] = -8;  // Blocking factor for multiple RHSs
00107   ICNTL[28] = 0;   // Automatic choice of sequential/parallel analysis phase
00108   ICNTL[29] = 0;   // Parallel analysis uses PT-SCOTCH or ParMetis? (none)
00109 
00110   // 30 - 40 are not used
00111 
00112 #else
00113   ICNTL[5]  = 0;   // Matrix is given in elemental (i.e. triplet) from
00114   ICNTL[6]  = 7;   // Choose column permutation automatically
00115   ICNTL[7]  = 7;   // Choose ordering method automatically
00116   ICNTL[8]  = 7;   // Choose scaling automatically
00117   ICNTL[8]  = 7;   // Choose scaling automatically
00118   ICNTL[9]  = 1;   // Compute A x = b
00119   ICNTL[10] = 0;   // Maximum steps of iterative refinement
00120   ICNTL[11] = 0;   // Do not collect statistics
00121   ICNTL[12] = 0;   // Use Node level parallelism
00122   ICNTL[13] = 0;   // Use ScaLAPACK for root node
00123   ICNTL[14] = 20;  // Increase memory allocation 20% at a time
00124   ICNTL[15] = 0;   // Minimize memory use (not flop count)
00125   ICNTL[16] = 0;   // Do not perform null space detection
00126   ICNTL[17] = 0;   // Unused (null space dimension)
00127 #endif
00128 
00129   Teuchos::ParameterList ParamList;
00130   SetParameters(ParamList);
00131 }
00132 
00133 //=============================================================================
00134 void Amesos_Mumps::Destroy()
00135 {
00136   if (!NoDestroy_) 
00137   { 
00138     // destroy instance of the package
00139     MDS.job = -2;
00140 
00141     if (Comm().MyPID() < MaxProcs_) dmumps_c(&(MDS));
00142     
00143 #if 0
00144     if (IsComputeSchurComplementOK_ && (Comm().MyPID() == 0)
00145   && MDS.schur) {
00146       delete [] MDS.schur;
00147       MDS.schur = 0;
00148     }
00149 #endif
00150     
00151 #ifdef HAVE_MPI
00152     if (MUMPSComm_) 
00153     {
00154       MPI_Comm_free( &MUMPSComm_ );
00155       MUMPSComm_ = 0;
00156     }
00157 #endif
00158     
00159     if( (verbose_ && PrintTiming_) || verbose_ == 2) PrintTiming();
00160     if( (verbose_ && PrintStatus_) || verbose_ == 2) PrintStatus();
00161   }
00162 }
00163 
00164 //=============================================================================
00165 Amesos_Mumps::~Amesos_Mumps(void)
00166 {
00167   Destroy();
00168 }
00169 
00170 //=============================================================================
00171 int Amesos_Mumps::ConvertToTriplet(const bool OnlyValues)
00172 {
00173 
00174   Epetra_RowMatrix* ptr;
00175   if (Comm().NumProc() == 1)
00176     ptr = &Matrix();
00177   else {
00178     ptr = &RedistrMatrix(true);
00179   }
00180 
00181   ResetTimer();
00182   
00183 #ifdef EXTRA_DEBUG_INFO
00184   Epetra_CrsMatrix* Eptr = dynamic_cast<Epetra_CrsMatrix*>( ptr );
00185   if ( ptr->NumGlobalNonzeros() < 300 ) SetICNTL(4,3 );  // Enable more debug info for small matrices
00186   if ( ptr->NumGlobalNonzeros() < 42 && Eptr ) { 
00187       std::cout << " Matrix = " << std::endl ; 
00188       Eptr->Print( std::cout ) ; 
00189   } else {
00190       assert( Eptr );
00191   }
00192 #endif
00193 
00194   Row.resize(ptr->NumMyNonzeros());
00195   Col.resize(ptr->NumMyNonzeros());
00196   Val.resize(ptr->NumMyNonzeros());
00197 
00198   int MaxNumEntries = ptr->MaxNumEntries();
00199   std::vector<int> Indices;
00200   std::vector<double> Values;
00201   Indices.resize(MaxNumEntries);
00202   Values.resize(MaxNumEntries);
00203 
00204   int count = 0;
00205 
00206   for (int i = 0; i < ptr->NumMyRows() ; ++i) {
00207 
00208     int GlobalRow = ptr->RowMatrixRowMap().GID(i);
00209 
00210     int NumEntries = 0;
00211     int ierr;
00212     ierr = ptr->ExtractMyRowCopy(i, MaxNumEntries,
00213            NumEntries, &Values[0],
00214            &Indices[0]);
00215     AMESOS_CHK_ERR(ierr);
00216 
00217     for (int j = 0 ; j < NumEntries ; ++j) {
00218       if (OnlyValues == false) {
00219   Row[count] = GlobalRow + 1;
00220   Col[count] = ptr->RowMatrixColMap().GID(Indices[j]) + 1;
00221       }
00222       
00223       // MS // Added on 15-Mar-05.
00224       if (AddToDiag_ && Indices[j] == i)
00225         Values[j] += AddToDiag_;
00226 
00227       Val[count] = Values[j];
00228       count++;
00229     }
00230   }
00231 
00232   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
00233   
00234   assert (count <= ptr->NumMyNonzeros());
00235 
00236   return(0);
00237 }
00238 
00239 //=============================================================================
00240 // NOTE: suppose first position is 1 (as in FORTRAN)
00241 void Amesos_Mumps::SetICNTL(int pos, int value)
00242 {
00243   ICNTL[pos] = value;
00244 }
00245 
00246 //=============================================================================
00247 // NOTE: suppose first position is 1 (as in FORTRAN)
00248 void Amesos_Mumps::SetCNTL(int pos, double value)
00249 {
00250   CNTL[pos] = value;
00251 }
00252 
00253 //=============================================================================
00254 void Amesos_Mumps::SetICNTLandCNTL()
00255 {
00256   std::map<int,int>::iterator i_iter;
00257   for (i_iter = ICNTL.begin() ; i_iter != ICNTL.end() ; ++i_iter)
00258   {
00259     int pos = i_iter->first;
00260     int val = i_iter->second;
00261     if (pos < 1 || pos > 40) continue;
00262     MDS.ICNTL(pos) = val;
00263   }
00264 
00265   std::map<int,double>::iterator d_iter;
00266   for (d_iter = CNTL.begin() ; d_iter != CNTL.end() ; ++d_iter)
00267   {
00268     int pos = d_iter->first;
00269     double val = d_iter->second;
00270     if (pos < 1 || pos > 5) continue;
00271     MDS.CNTL(pos) = val;
00272   }
00273 
00274   // fix some options
00275   
00276   if (Comm().NumProc() != 1) MDS.ICNTL(18)= 3;
00277   else                       MDS.ICNTL(18)= 0;
00278   
00279   if (UseTranspose())  MDS.ICNTL(9) = 0; 
00280   else                 MDS.ICNTL(9) = 1;
00281   
00282   MDS.ICNTL(5) = 0;
00283 
00284 #if 0
00285   if (IsComputeSchurComplementOK_) MDS.ICNTL(19) = 1;
00286   else                             MDS.ICNTL(19) = 0;
00287 
00288   // something to do if the Schur complement is required.
00289   if (IsComputeSchurComplementOK_ && Comm().MyPID() == 0) 
00290   {
00291     MDS.size_schur = NumSchurComplementRows_;
00292     MDS.listvar_schur = SchurComplementRows_;
00293     MDS.schur = new double[NumSchurComplementRows_*NumSchurComplementRows_];
00294   }
00295 #endif
00296 }
00297 
00298 //=============================================================================
00299 int Amesos_Mumps::SetParameters( Teuchos::ParameterList & ParameterList)
00300 {
00301   // ========================================= //
00302   // retrive MUMPS' parameters from list.      //
00303   // default values defined in the constructor //
00304   // ========================================= //
00305   
00306   // retrive general parameters
00307 
00308   SetStatusParameters(ParameterList);
00309 
00310   SetControlParameters(ParameterList);
00311 
00312   if (ParameterList.isParameter("NoDestroy"))
00313     NoDestroy_ = ParameterList.get<bool>("NoDestroy");
00314   
00315   // can be use using mumps sublist, via "ICNTL(2)"
00316   //  if (debug_) 
00317   //    MDS.ICNTL(2) = 6; // Turn on Mumps verbose output 
00318 
00319   // retrive MUMPS' specific parameters
00320   
00321   if (ParameterList.isSublist("mumps")) 
00322   {
00323     Teuchos::ParameterList MumpsParams = ParameterList.sublist("mumps") ;
00324 
00325     // ICNTL
00326     for (int i = 1 ; i <= 40 ; ++i)
00327     {
00328       char what[80];
00329       sprintf(what, "ICNTL(%d)", i);
00330       if (MumpsParams.isParameter(what)) 
00331         SetICNTL(i, MumpsParams.get<int>(what));
00332     }
00333 
00334     // CNTL
00335     for (int i = 1 ; i <= 5 ; ++i)
00336     {
00337       char what[80];
00338       sprintf(what, "CNTL(%d)", i);
00339       if (MumpsParams.isParameter(what)) 
00340         SetCNTL(i, MumpsParams.get<double>(what));
00341     }
00342 
00343     // ordering
00344      if (MumpsParams.isParameter("PermIn")) 
00345      {
00346       int* PermIn = 0;
00347       PermIn = MumpsParams.get<int*>("PermIn");
00348       if (PermIn) SetOrdering(PermIn);
00349      }
00350      // Col scaling
00351      if (MumpsParams.isParameter("ColScaling")) 
00352      {
00353        double * ColSca = 0;
00354        ColSca = MumpsParams.get<double *>("ColScaling");
00355        if (ColSca) SetColScaling(ColSca);
00356      }
00357      // Row scaling
00358      if (MumpsParams.isParameter("RowScaling")) 
00359      {
00360        double * RowSca = 0;
00361        RowSca = MumpsParams.get<double *>("RowScaling");
00362        if (RowSca) SetRowScaling(RowSca);
00363      }
00364      // that's all folks
00365   }  
00366 
00367   return(0);
00368 }
00369 
00370 //=============================================================================
00371 
00372 void Amesos_Mumps::CheckParameters() 
00373 {
00374 #ifndef HAVE_AMESOS_MPI_C2F
00375   MaxProcs_ = -3;
00376 #endif
00377 
00378   // check parameters and fix values of MaxProcs_
00379 
00380   int NumGlobalNonzeros, NumRows;
00381   
00382   NumGlobalNonzeros = Matrix().NumGlobalNonzeros(); 
00383   NumRows = Matrix().NumGlobalRows(); 
00384 
00385   // optimal value for MaxProcs == -1
00386   
00387   int OptNumProcs1 = 1 + EPETRA_MAX(NumRows/10000, NumGlobalNonzeros/100000);
00388   OptNumProcs1 = EPETRA_MIN(Comm().NumProc(),OptNumProcs1);
00389 
00390   // optimal value for MaxProcs == -2
00391 
00392   int OptNumProcs2 = (int)sqrt(1.0 *  Comm().NumProc());
00393   if (OptNumProcs2 < 1) OptNumProcs2 = 1;
00394 
00395   // fix the value of MaxProcs
00396 
00397   switch (MaxProcs_) {
00398   case -1:
00399     MaxProcs_ = OptNumProcs1;
00400     break;
00401   case -2:
00402     MaxProcs_ = OptNumProcs2;
00403     break;
00404   case -3:
00405     MaxProcs_ = Comm().NumProc();
00406     break;
00407   }
00408 
00409   // few checks
00410   if (MaxProcs_ > Comm().NumProc()) MaxProcs_ = Comm().NumProc();
00411 //  if ( MaxProcs_ > 1 ) MaxProcs_ = Comm().NumProc();     // Bug - bogus kludge here  - didn't work anyway
00412 }
00413 
00414 //=============================================================================
00415 int Amesos_Mumps::SymbolicFactorization()
00416 {
00417 
00418   // erase data if present. 
00419   if (IsSymbolicFactorizationOK_ && MDS.job != -777)
00420    Destroy();
00421 
00422   IsSymbolicFactorizationOK_ = false;
00423   IsNumericFactorizationOK_ = false;
00424 
00425   CreateTimer(Comm());
00426   
00427   CheckParameters();
00428   AMESOS_CHK_ERR(ConvertToTriplet(false));
00429 
00430 #if defined(HAVE_MPI) && defined(HAVE_AMESOS_MPI_C2F)
00431   if (MaxProcs_ != Comm().NumProc()) 
00432   {
00433     if(MUMPSComm_) 
00434       MPI_Comm_free(&MUMPSComm_);
00435 
00436     std::vector<int> ProcsInGroup(MaxProcs_);
00437     for (int i = 0 ; i < MaxProcs_ ; ++i) 
00438       ProcsInGroup[i] = i;
00439 
00440     MPI_Group OrigGroup, MumpsGroup;
00441     MPI_Comm_group(MPI_COMM_WORLD, &OrigGroup);
00442     MPI_Group_incl(OrigGroup, MaxProcs_, &ProcsInGroup[0], &MumpsGroup);
00443     MPI_Comm_create(MPI_COMM_WORLD, MumpsGroup, &MUMPSComm_);
00444 
00445 #ifdef MUMPS_4_9
00446     MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f( MUMPSComm_);
00447 #else
00448 
00449 #ifndef HAVE_AMESOS_OLD_MUMPS
00450     MDS.comm_fortran = (DMUMPS_INT) MPI_Comm_c2f( MUMPSComm_);
00451 #else
00452     MDS.comm_fortran = (F_INT) MPI_Comm_c2f( MUMPSComm_);
00453 #endif
00454 
00455 #endif
00456 
00457   } 
00458   else 
00459   {
00460     const Epetra_MpiComm* MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
00461     assert (MpiComm != 0);
00462 #ifdef MUMPS_4_9
00463     MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
00464 #else
00465 
00466 #ifndef HAVE_AMESOS_OLD_MUMPS
00467     MDS.comm_fortran = (DMUMPS_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
00468 #else
00469     MDS.comm_fortran = (F_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
00470 #endif
00471 
00472 #endif
00473   }
00474 #else
00475 #ifdef IFPACK_SUBCOMM_CODE
00476   // This next three lines of code were required to make Amesos_Mumps work
00477   // with Ifpack_SubdomainFilter. They may actually be usefull in all cases
00478   // when using MUMPS on a subdomain.
00479   const Epetra_MpiComm* MpiComm = dynamic_cast<const Epetra_MpiComm*>(&Comm());
00480   assert (MpiComm != 0);
00481   MDS.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(MpiComm->GetMpiComm());
00482 #else
00483   // only thing I can do, use MPI_COMM_WORLD. This will work in serial as well
00484   MDS.comm_fortran = -987654;
00485 #endif
00486 #endif
00487   
00488   MDS.job = -1  ;     //  Initialization
00489   MDS.par = 1 ;       //  Host IS involved in computations
00490 //  MDS.sym = MatrixProperty_;
00491   MDS.sym =  0;       //  MatrixProperty_ is unititalized.  Furthermore MUMPS 
00492                       //  expects only half of the matrix to be provided for
00493                       //  symmetric matrices.  Hence setting MDS.sym to be non-zero
00494                       //  indicating that the matrix is symmetric will only work
00495                       //  if we change ConvertToTriplet to pass only half of the 
00496                       //  matrix.  Bug #2331 and Bug #2332 - low priority
00497 
00498 
00499   RedistrMatrix(true);
00500 
00501   if (Comm().MyPID() < MaxProcs_) 
00502   {
00503     dmumps_c(&(MDS));   //  Initialize MUMPS
00504     static_cast<void>( CheckError( ) );  
00505   }
00506 
00507   MDS.n = Matrix().NumGlobalRows();
00508 
00509   // fix pointers for nonzero pattern of A. Numerical values
00510   // will be entered in PerformNumericalFactorization()
00511   if (Comm().NumProc() != 1) 
00512   {
00513     MDS.nz_loc = RedistrMatrix().NumMyNonzeros();
00514 
00515     if (Comm().MyPID() < MaxProcs_) 
00516     {
00517       MDS.irn_loc = &Row[0]; 
00518       MDS.jcn_loc = &Col[0];
00519     }
00520   } 
00521   else 
00522   {
00523     if (Comm().MyPID() == 0) 
00524     {
00525       MDS.nz = Matrix().NumMyNonzeros();
00526       MDS.irn = &Row[0]; 
00527       MDS.jcn = &Col[0]; 
00528     }
00529   }
00530 
00531   // scaling if provided by the user
00532   if (RowSca_ != 0) 
00533   {
00534     MDS.rowsca = RowSca_;
00535     MDS.colsca = ColSca_;
00536   }
00537 
00538   // given ordering if provided by the user
00539   if (PermIn_ != 0) 
00540   {
00541     MDS.perm_in = PermIn_;
00542   }
00543 
00544   MDS.job = 1;     // Request symbolic factorization
00545 
00546   SetICNTLandCNTL();
00547 
00548   // Perform symbolic factorization
00549 
00550   ResetTimer();
00551 
00552   if (Comm().MyPID() < MaxProcs_) 
00553     dmumps_c(&(MDS));
00554 
00555   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
00556 
00557   int IntWrong = CheckError()?1:0 ; 
00558   int AnyWrong;
00559   Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ; 
00560   bool Wrong = AnyWrong > 0 ; 
00561 
00562 
00563   if ( Wrong ) {
00564       AMESOS_CHK_ERR( StructurallySingularMatrixError ) ; 
00565   }
00566 
00567   IsSymbolicFactorizationOK_ = true ;
00568   NumSymbolicFact_++;  
00569 
00570   return 0;
00571 }
00572 
00573 //=============================================================================
00574 
00575 int Amesos_Mumps::NumericFactorization()
00576 {
00577   IsNumericFactorizationOK_ = false;
00578   
00579   if (IsSymbolicFactorizationOK_ == false)
00580     AMESOS_CHK_ERR(SymbolicFactorization());
00581 
00582   RedistrMatrix(true);
00583   AMESOS_CHK_ERR(ConvertToTriplet(true));
00584 
00585   if (Comm().NumProc() != 1) 
00586   {
00587     if (Comm().MyPID() < MaxProcs_) 
00588       MDS.a_loc = &Val[0];
00589   } 
00590   else 
00591     MDS.a = &Val[0];
00592 
00593   // Request numeric factorization 
00594   MDS.job = 2;
00595   // Perform numeric factorization
00596   ResetTimer();
00597 
00598   if (Comm().MyPID() < MaxProcs_) {
00599     dmumps_c(&(MDS));
00600   }
00601 
00602   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
00603   
00604   int IntWrong = CheckError()?1:0 ; 
00605   int AnyWrong;
00606   Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ; 
00607   bool Wrong = AnyWrong > 0 ; 
00608 
00609 
00610   if ( Wrong ) {
00611       AMESOS_CHK_ERR( NumericallySingularMatrixError ) ; 
00612   }
00613 
00614   IsNumericFactorizationOK_ = true;
00615   NumNumericFact_++;  
00616   return(0);
00617 }
00618 
00619 //=============================================================================
00620 
00621 int Amesos_Mumps::Solve()
00622 { 
00623   if (IsNumericFactorizationOK_ == false)
00624     AMESOS_CHK_ERR(NumericFactorization());
00625 
00626   Epetra_MultiVector* vecX = Problem_->GetLHS(); 
00627   Epetra_MultiVector* vecB = Problem_->GetRHS();
00628   int NumVectors = vecX->NumVectors();
00629 
00630   if ((vecX == 0) || (vecB == 0))
00631     AMESOS_CHK_ERR(-1);
00632 
00633   if (NumVectors != vecB->NumVectors())
00634     AMESOS_CHK_ERR(-1);
00635 
00636   if (Comm().NumProc() == 1) 
00637   {
00638     // do not import any data
00639     for (int j = 0 ; j < NumVectors; j++) 
00640     {
00641       ResetTimer();
00642 
00643       MDS.job = 3;     // Request solve
00644 
00645       for (int i = 0 ; i < Matrix().NumMyRows() ; ++i) 
00646   (*vecX)[j][i] = (*vecB)[j][i];
00647       MDS.rhs = (*vecX)[j];
00648 
00649       dmumps_c(&(MDS)) ;  // Perform solve
00650       static_cast<void>( CheckError( ) );   // Can hang 
00651       SolveTime_ = AddTime("Total solve time", SolveTime_);
00652     }
00653   } 
00654   else 
00655   {
00656     Epetra_MultiVector SerialVector(SerialMap(),NumVectors);
00657 
00658     ResetTimer();
00659     AMESOS_CHK_ERR(SerialVector.Import(*vecB,SerialImporter(),Insert));
00660     VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00661     
00662     for (int j = 0 ; j < NumVectors; j++) 
00663     {
00664       if (Comm().MyPID() == 0)
00665   MDS.rhs = SerialVector[j];
00666 
00667       // solve the linear system and take time
00668       MDS.job = 3;     
00669       ResetTimer();
00670       if (Comm().MyPID() < MaxProcs_) 
00671   dmumps_c(&(MDS)) ;  // Perform solve
00672       static_cast<void>( CheckError( ) );   // Can hang 
00673 
00674       SolveTime_ = AddTime("Total solve time", SolveTime_);
00675     }
00676 
00677     // ship solution back and take timing
00678     ResetTimer();
00679     AMESOS_CHK_ERR(vecX->Export(SerialVector,SerialImporter(),Insert));
00680     VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00681   }
00682 
00683   if (ComputeTrueResidual_)
00684     ComputeTrueResidual(Matrix(), *vecX, *vecB, UseTranspose(), "Amesos_Mumps");
00685 
00686   if (ComputeVectorNorms_)
00687     ComputeVectorNorms(*vecX, *vecB, "Amesos_Mumps");
00688 
00689   NumSolve_++;  
00690   
00691   return(0) ; 
00692 }
00693 
00694 #if 0
00695 //=============================================================================
00696 Epetra_CrsMatrix * Amesos_Mumps::GetCrsSchurComplement() 
00697 {
00698 
00699   if( IsComputeSchurComplementOK_ ) {
00700     
00701     if( Comm().MyPID() != 0 ) NumSchurComplementRows_ = 0;
00702     
00703     Epetra_Map SCMap(-1,NumSchurComplementRows_, 0, Comm());
00704 
00705     CrsSchurComplement_ = new Epetra_CrsMatrix(Copy,SCMap,NumSchurComplementRows_);
00706 
00707     if( Comm().MyPID() == 0 )
00708       for( int i=0 ; i<NumSchurComplementRows_ ; ++i ) {
00709   for( int j=0 ; j<NumSchurComplementRows_ ; ++j ) {
00710     int pos = i+ j *NumSchurComplementRows_;
00711     CrsSchurComplement_->InsertGlobalValues(i,1,&(MDS.schur[pos]),&j);
00712   }
00713       }
00714     
00715     CrsSchurComplement_->FillComplete();
00716 
00717     return CrsSchurComplement_;
00718   }
00719 
00720   return 0;
00721   
00722 }
00723 
00724 //=============================================================================
00725 Epetra_SerialDenseMatrix * Amesos_Mumps::GetDenseSchurComplement() 
00726 {
00727   
00728   if( IsComputeSchurComplementOK_ ) {
00729     
00730     if( Comm().MyPID() != 0 ) return 0;
00731     
00732     DenseSchurComplement_ = new Epetra_SerialDenseMatrix(NumSchurComplementRows_,
00733               NumSchurComplementRows_);
00734     
00735     for( int i=0 ; i<NumSchurComplementRows_ ; ++i ) {
00736       for( int j=0 ; j<NumSchurComplementRows_ ; ++j ) {
00737   int pos = i+ j *NumSchurComplementRows_;
00738   (*DenseSchurComplement_)(i,j) = MDS.schur[pos];
00739       }
00740     }
00741     
00742     return DenseSchurComplement_;
00743     
00744   }
00745   
00746   return(0);
00747 }
00748 
00749 //=============================================================================
00750 int Amesos_Mumps::ComputeSchurComplement(bool flag, int NumSchurComplementRows,
00751            int * SchurComplementRows)
00752 {
00753   NumSchurComplementRows_ = NumSchurComplementRows;
00754   
00755   SchurComplementRows_ = SchurComplementRows;
00756 
00757   // modify because MUMPS is fortran-driven
00758   if( Comm().MyPID() == 0 )
00759     for( int i=0 ; i<NumSchurComplementRows ; ++i ) SchurComplementRows_[i]++;
00760   
00761   IsComputeSchurComplementOK_ = flag;
00762 
00763   return 0;
00764 }
00765 #endif
00766 
00767 //=============================================================================
00768 void Amesos_Mumps::PrintStatus() const
00769 {
00770 
00771   if (Comm().MyPID() != 0 ) return;
00772 
00773   //  The following lines are commented out to deal with bug #1887 - kss
00774 #ifndef IRIX64
00775   PrintLine();
00776   std::cout << "Amesos_Mumps : Matrix has " << Matrix().NumGlobalRows() << " rows"
00777        << " and " << Matrix().NumGlobalNonzeros() << " nonzeros" << std::endl;
00778   std::cout << "Amesos_Mumps : Nonzero elements per row = "
00779        << 1.0*Matrix().NumGlobalNonzeros()/Matrix().NumGlobalRows() << std::endl;
00780   std::cout << "Amesos_Mumps : Percentage of nonzero elements = "
00781        << 100.0*Matrix().NumGlobalNonzeros()/(pow(Matrix().NumGlobalRows(),2.0)) << std::endl;
00782   std::cout << "Amesos_Mumps : Use transpose = " << UseTranspose_ << std::endl;
00783 //  MatrixProperty_ is unused - see bug #2331 and bug #2332 in this file and bugzilla
00784   if (MatrixProperty_ == 0) std::cout << "Amesos_Mumps : Matrix is general unsymmetric" << std::endl;
00785   if (MatrixProperty_ == 2) std::cout << "Amesos_Mumps : Matrix is general symmetric" << std::endl;
00786   if (MatrixProperty_ == 1) std::cout << "Amesos_Mumps : Matrix is SPD" << std::endl;
00787   std::cout << "Amesos_Mumps : Available process(es) = " << Comm().NumProc() << std::endl;
00788   std::cout << "Amesos_Mumps : Using " << MaxProcs_ << " process(es)" << std::endl;
00789   
00790   std::cout << "Amesos_Mumps : Estimated FLOPS for elimination = "
00791        << MDS.RINFOG(1) << std::endl;
00792   std::cout << "Amesos_Mumps : Total FLOPS for assembly = "
00793        << MDS.RINFOG(2) << std::endl;
00794   std::cout << "Amesos_Mumps : Total FLOPS for elimination = "
00795        << MDS.RINFOG(3) << std::endl;
00796   
00797   std::cout << "Amesos_Mumps : Total real space to store the LU factors = "
00798        << MDS.INFOG(9) << std::endl;
00799   std::cout << "Amesos_Mumps : Total integer space to store the LU factors = "
00800        << MDS.INFOG(10) << std::endl;
00801   std::cout << "Amesos_Mumps : Total number of iterative steps refinement = "
00802        << MDS.INFOG(15) << std::endl;
00803   std::cout << "Amesos_Mumps : Estimated size of MUMPS internal data\n"
00804        << "Amesos_Mumps : for running factorization = "
00805        << MDS.INFOG(16) << " Mbytes" << std::endl;
00806   std::cout << "Amesos_Mumps : for running factorization = "
00807        << MDS.INFOG(17) << " Mbytes" << std::endl;
00808   std::cout << "Amesos_Mumps : Allocated during factorization = "
00809        << MDS.INFOG(19) << " Mbytes" << std::endl;
00810   PrintLine();
00811 #endif
00812 }
00813 
00814 //=============================================================================
00815 int Amesos_Mumps::CheckError() 
00816 {
00817   bool Wrong = ((MDS.INFOG(1) != 0) || (MDS.INFO(1) != 0))
00818                && (Comm().MyPID() < MaxProcs_);
00819   
00820   // an error occurred in MUMPS. Print out information and quit.
00821 
00822   if (Comm().MyPID() == 0 && Wrong) 
00823   {
00824     std::cerr << "Amesos_Mumps : ERROR" << std::endl;
00825     std::cerr << "Amesos_Mumps : INFOG(1) = " << MDS.INFOG(1) << std::endl;
00826     std::cerr << "Amesos_Mumps : INFOG(2) = " << MDS.INFOG(2) << std::endl;
00827   }
00828   
00829   if (MDS.INFO(1) != 0 && Wrong) 
00830   {
00831     std::cerr << "Amesos_Mumps : On process " << Comm().MyPID()
00832    << ", INFO(1) = " << MDS.INFO(1) << std::endl;
00833     std::cerr << "Amesos_Mumps : On process " << Comm().MyPID()
00834    << ", INFO(2) = " << MDS.INFO(2) << std::endl;
00835   }
00836 
00837   if (Wrong) 
00838       return 1 ;
00839   else
00840       return 0 ;
00841 }
00842 
00843 // ======================================================================
00844 void Amesos_Mumps::PrintTiming() const
00845 {
00846   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00847     return;
00848 
00849   double ConTime = GetTime(MtxConvTime_);
00850   double MatTime = GetTime(MtxRedistTime_);
00851   double VecTime = GetTime(VecRedistTime_);
00852   double SymTime = GetTime(SymFactTime_);
00853   double NumTime = GetTime(SymFactTime_);
00854   double SolTime = GetTime(SolveTime_);
00855 
00856   if (NumSymbolicFact_) SymTime /= NumSymbolicFact_;
00857   if (NumNumericFact_)  NumTime /= NumNumericFact_;
00858   if (NumSolve_)        SolTime /= NumSolve_;
00859 
00860   std::string p = "Amesos_Mumps : ";
00861   PrintLine();
00862 
00863   std::cout << p << "Time to convert matrix to MUMPS format = "
00864        << ConTime << " (s)" << std::endl;
00865   std::cout << p << "Time to redistribute matrix = "
00866        << MatTime << " (s)" << std::endl;
00867   std::cout << p << "Time to redistribute vectors = "
00868        << VecTime << " (s)" << std::endl;
00869   std::cout << p << "Number of symbolic factorizations = "
00870        << NumSymbolicFact_ << std::endl;
00871   std::cout << p << "Time for sym fact = "
00872        << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
00873   std::cout << p << "Number of numeric factorizations = "
00874        << NumNumericFact_ << std::endl;
00875   std::cout << p << "Time for num fact = "
00876        << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
00877   std::cout << p << "Number of solve phases = "
00878        << NumSolve_ << std::endl;
00879   std::cout << p << "Time for solve = "
00880        << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
00881 
00882   PrintLine();
00883 }
00884 
00885 // ===========================================================================
00886 Epetra_RowMatrix& Amesos_Mumps::Matrix() 
00887 {
00888   Epetra_RowMatrix* Matrix = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
00889   assert (Matrix != 0);
00890   return(*Matrix);
00891 }
00892 
00893 // ===========================================================================
00894 const Epetra_RowMatrix& Amesos_Mumps::Matrix() const
00895 {
00896   Epetra_RowMatrix* Matrix = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
00897   assert (Matrix != 0);
00898   return(*Matrix);
00899 }
00900 
00901 // ===========================================================================
00902 Epetra_Map& Amesos_Mumps::RedistrMap() 
00903 {
00904   assert (Comm().NumProc() != 1);
00905 
00906   if (RedistrMap_ == Teuchos::null) {
00907     int i = Matrix().NumGlobalRows() / MaxProcs_;
00908     if (Comm().MyPID() == 0)
00909       i += Matrix().NumGlobalRows() % MaxProcs_;
00910     else if (Comm().MyPID() >= MaxProcs_)
00911       i = 0;
00912 
00913     RedistrMap_ = rcp(new Epetra_Map(Matrix().NumGlobalRows(),i,0,Comm()));
00914     assert (RedistrMap_.get() != 0);
00915   }
00916   return(*RedistrMap_);
00917 }
00918 
00919 // ===========================================================================
00920 Epetra_Import& Amesos_Mumps::RedistrImporter()
00921 {
00922   assert (Comm().NumProc() != 1);
00923 
00924   if (RedistrImporter_ == null) {
00925     RedistrImporter_ = rcp(new Epetra_Import(RedistrMap(),Matrix().RowMatrixRowMap()));
00926     assert (RedistrImporter_.get() != 0);
00927   }
00928   return(*RedistrImporter_);
00929 }
00930 
00931 // ===========================================================================
00932 Epetra_RowMatrix& Amesos_Mumps::RedistrMatrix(const bool ImportMatrix)
00933 {
00934   if (Comm().NumProc() == 1) return(Matrix());
00935 
00936   if (ImportMatrix) RedistrMatrix_ = null;
00937 
00938   if (RedistrMatrix_ == null) 
00939   {
00940     RedistrMatrix_ = rcp(new Epetra_CrsMatrix(Copy,RedistrMap(),0));
00941     if (ImportMatrix) {
00942       int ierr = RedistrMatrix_->Import(Matrix(),RedistrImporter(),Insert);
00943       assert (ierr == 0);
00944       ierr = RedistrMatrix_->FillComplete();
00945       assert (ierr == 0);
00946     }
00947   }
00948 
00949   return(*RedistrMatrix_);
00950 }
00951 
00952 // ===========================================================================
00953 Epetra_Map& Amesos_Mumps::SerialMap()
00954 {
00955   if (SerialMap_ == null) 
00956   {
00957     int i = Matrix().NumGlobalRows();
00958     if (Comm().MyPID()) i = 0;
00959     SerialMap_ = rcp(new Epetra_Map(-1,i,0,Comm()));
00960     assert (SerialMap_ != null);
00961   }
00962   return(*SerialMap_);
00963 }
00964 
00965 // ===========================================================================
00966 Epetra_Import& Amesos_Mumps::SerialImporter()
00967 { 
00968   if (SerialImporter_ == null) {
00969     SerialImporter_ = rcp(new Epetra_Import(SerialMap(),Matrix().OperatorDomainMap()));
00970     assert (SerialImporter_ != null);
00971   }
00972   return(*SerialImporter_);
00973 }
00974 
00975 // ===========================================================================
00976 double* Amesos_Mumps::GetRINFO() 
00977 {
00978   return ( MDS.rinfo);
00979 }
00980 
00981 // ===========================================================================
00982 int* Amesos_Mumps::GetINFO() 
00983 {
00984   return (MDS.info);
00985 }
00986 
00987 // ===========================================================================
00988 double* Amesos_Mumps::GetRINFOG()
00989 {
00990   return (MDS.rinfog);
00991 }
00992 
00993 // ===========================================================================
00994 int* Amesos_Mumps::GetINFOG()
00995 {
00996   return (MDS.infog);
00997 }
00998 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines