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   // only thing I can do, use MPI_COMM_WORLD. This will work in serial as well
00476   MDS.comm_fortran = -987654;
00477 #endif
00478   
00479   MDS.job = -1  ;     //  Initialization
00480   MDS.par = 1 ;       //  Host IS involved in computations
00481 //  MDS.sym = MatrixProperty_;
00482   MDS.sym =  0;       //  MatrixProperty_ is unititalized.  Furthermore MUMPS 
00483                       //  expects only half of the matrix to be provided for
00484                       //  symmetric matrices.  Hence setting MDS.sym to be non-zero
00485                       //  indicating that the matrix is symmetric will only work
00486                       //  if we change ConvertToTriplet to pass only half of the 
00487                       //  matrix.  Bug #2331 and Bug #2332 - low priority
00488 
00489 
00490   RedistrMatrix(true);
00491 
00492   if (Comm().MyPID() < MaxProcs_) 
00493   {
00494     dmumps_c(&(MDS));   //  Initialize MUMPS
00495     static_cast<void>( CheckError( ) );  
00496   }
00497 
00498   MDS.n = Matrix().NumGlobalRows();
00499 
00500   // fix pointers for nonzero pattern of A. Numerical values
00501   // will be entered in PerformNumericalFactorization()
00502   if (Comm().NumProc() != 1) 
00503   {
00504     MDS.nz_loc = RedistrMatrix().NumMyNonzeros();
00505 
00506     if (Comm().MyPID() < MaxProcs_) 
00507     {
00508       MDS.irn_loc = &Row[0]; 
00509       MDS.jcn_loc = &Col[0];
00510     }
00511   } 
00512   else 
00513   {
00514     if (Comm().MyPID() == 0) 
00515     {
00516       MDS.nz = Matrix().NumMyNonzeros();
00517       MDS.irn = &Row[0]; 
00518       MDS.jcn = &Col[0]; 
00519     }
00520   }
00521 
00522   // scaling if provided by the user
00523   if (RowSca_ != 0) 
00524   {
00525     MDS.rowsca = RowSca_;
00526     MDS.colsca = ColSca_;
00527   }
00528 
00529   // given ordering if provided by the user
00530   if (PermIn_ != 0) 
00531   {
00532     MDS.perm_in = PermIn_;
00533   }
00534 
00535   MDS.job = 1;     // Request symbolic factorization
00536 
00537   SetICNTLandCNTL();
00538 
00539   // Perform symbolic factorization
00540 
00541   ResetTimer();
00542 
00543   if (Comm().MyPID() < MaxProcs_) 
00544     dmumps_c(&(MDS));
00545 
00546   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
00547 
00548   int IntWrong = CheckError()?1:0 ; 
00549   int AnyWrong;
00550   Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ; 
00551   bool Wrong = AnyWrong > 0 ; 
00552 
00553 
00554   if ( Wrong ) {
00555       AMESOS_CHK_ERR( StructurallySingularMatrixError ) ; 
00556   }
00557 
00558   IsSymbolicFactorizationOK_ = true ;
00559   NumSymbolicFact_++;  
00560 
00561   return 0;
00562 }
00563 
00564 //=============================================================================
00565 
00566 int Amesos_Mumps::NumericFactorization()
00567 {
00568   IsNumericFactorizationOK_ = false;
00569   
00570   if (IsSymbolicFactorizationOK_ == false)
00571     AMESOS_CHK_ERR(SymbolicFactorization());
00572 
00573   RedistrMatrix(true);
00574   AMESOS_CHK_ERR(ConvertToTriplet(true));
00575 
00576   if (Comm().NumProc() != 1) 
00577   {
00578     if (Comm().MyPID() < MaxProcs_) 
00579       MDS.a_loc = &Val[0];
00580   } 
00581   else 
00582     MDS.a = &Val[0];
00583 
00584   // Request numeric factorization 
00585   MDS.job = 2;
00586   // Perform numeric factorization
00587   ResetTimer();
00588 
00589   if (Comm().MyPID() < MaxProcs_) {
00590     dmumps_c(&(MDS));
00591   }
00592 
00593   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
00594   
00595   int IntWrong = CheckError()?1:0 ; 
00596   int AnyWrong;
00597   Comm().SumAll( &IntWrong, &AnyWrong, 1 ) ; 
00598   bool Wrong = AnyWrong > 0 ; 
00599 
00600 
00601   if ( Wrong ) {
00602       AMESOS_CHK_ERR( NumericallySingularMatrixError ) ; 
00603   }
00604 
00605   IsNumericFactorizationOK_ = true;
00606   NumNumericFact_++;  
00607   return(0);
00608 }
00609 
00610 //=============================================================================
00611 
00612 int Amesos_Mumps::Solve()
00613 { 
00614   if (IsNumericFactorizationOK_ == false)
00615     AMESOS_CHK_ERR(NumericFactorization());
00616 
00617   Epetra_MultiVector* vecX = Problem_->GetLHS(); 
00618   Epetra_MultiVector* vecB = Problem_->GetRHS();
00619   int NumVectors = vecX->NumVectors();
00620 
00621   if ((vecX == 0) || (vecB == 0))
00622     AMESOS_CHK_ERR(-1);
00623 
00624   if (NumVectors != vecB->NumVectors())
00625     AMESOS_CHK_ERR(-1);
00626 
00627   if (Comm().NumProc() == 1) 
00628   {
00629     // do not import any data
00630     for (int j = 0 ; j < NumVectors; j++) 
00631     {
00632       ResetTimer();
00633 
00634       MDS.job = 3;     // Request solve
00635 
00636       for (int i = 0 ; i < Matrix().NumMyRows() ; ++i) 
00637   (*vecX)[j][i] = (*vecB)[j][i];
00638       MDS.rhs = (*vecX)[j];
00639 
00640       dmumps_c(&(MDS)) ;  // Perform solve
00641       static_cast<void>( CheckError( ) );   // Can hang 
00642       SolveTime_ = AddTime("Total solve time", SolveTime_);
00643     }
00644   } 
00645   else 
00646   {
00647     Epetra_MultiVector SerialVector(SerialMap(),NumVectors);
00648 
00649     ResetTimer();
00650     AMESOS_CHK_ERR(SerialVector.Import(*vecB,SerialImporter(),Insert));
00651     VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00652     
00653     for (int j = 0 ; j < NumVectors; j++) 
00654     {
00655       if (Comm().MyPID() == 0)
00656   MDS.rhs = SerialVector[j];
00657 
00658       // solve the linear system and take time
00659       MDS.job = 3;     
00660       ResetTimer();
00661       if (Comm().MyPID() < MaxProcs_) 
00662   dmumps_c(&(MDS)) ;  // Perform solve
00663       static_cast<void>( CheckError( ) );   // Can hang 
00664 
00665       SolveTime_ = AddTime("Total solve time", SolveTime_);
00666     }
00667 
00668     // ship solution back and take timing
00669     ResetTimer();
00670     AMESOS_CHK_ERR(vecX->Export(SerialVector,SerialImporter(),Insert));
00671     VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
00672   }
00673 
00674   if (ComputeTrueResidual_)
00675     ComputeTrueResidual(Matrix(), *vecX, *vecB, UseTranspose(), "Amesos_Mumps");
00676 
00677   if (ComputeVectorNorms_)
00678     ComputeVectorNorms(*vecX, *vecB, "Amesos_Mumps");
00679 
00680   NumSolve_++;  
00681   
00682   return(0) ; 
00683 }
00684 
00685 #if 0
00686 //=============================================================================
00687 Epetra_CrsMatrix * Amesos_Mumps::GetCrsSchurComplement() 
00688 {
00689 
00690   if( IsComputeSchurComplementOK_ ) {
00691     
00692     if( Comm().MyPID() != 0 ) NumSchurComplementRows_ = 0;
00693     
00694     Epetra_Map SCMap(-1,NumSchurComplementRows_, 0, Comm());
00695 
00696     CrsSchurComplement_ = new Epetra_CrsMatrix(Copy,SCMap,NumSchurComplementRows_);
00697 
00698     if( Comm().MyPID() == 0 )
00699       for( int i=0 ; i<NumSchurComplementRows_ ; ++i ) {
00700   for( int j=0 ; j<NumSchurComplementRows_ ; ++j ) {
00701     int pos = i+ j *NumSchurComplementRows_;
00702     CrsSchurComplement_->InsertGlobalValues(i,1,&(MDS.schur[pos]),&j);
00703   }
00704       }
00705     
00706     CrsSchurComplement_->FillComplete();
00707 
00708     return CrsSchurComplement_;
00709   }
00710 
00711   return 0;
00712   
00713 }
00714 
00715 //=============================================================================
00716 Epetra_SerialDenseMatrix * Amesos_Mumps::GetDenseSchurComplement() 
00717 {
00718   
00719   if( IsComputeSchurComplementOK_ ) {
00720     
00721     if( Comm().MyPID() != 0 ) return 0;
00722     
00723     DenseSchurComplement_ = new Epetra_SerialDenseMatrix(NumSchurComplementRows_,
00724               NumSchurComplementRows_);
00725     
00726     for( int i=0 ; i<NumSchurComplementRows_ ; ++i ) {
00727       for( int j=0 ; j<NumSchurComplementRows_ ; ++j ) {
00728   int pos = i+ j *NumSchurComplementRows_;
00729   (*DenseSchurComplement_)(i,j) = MDS.schur[pos];
00730       }
00731     }
00732     
00733     return DenseSchurComplement_;
00734     
00735   }
00736   
00737   return(0);
00738 }
00739 
00740 //=============================================================================
00741 int Amesos_Mumps::ComputeSchurComplement(bool flag, int NumSchurComplementRows,
00742            int * SchurComplementRows)
00743 {
00744   NumSchurComplementRows_ = NumSchurComplementRows;
00745   
00746   SchurComplementRows_ = SchurComplementRows;
00747 
00748   // modify because MUMPS is fortran-driven
00749   if( Comm().MyPID() == 0 )
00750     for( int i=0 ; i<NumSchurComplementRows ; ++i ) SchurComplementRows_[i]++;
00751   
00752   IsComputeSchurComplementOK_ = flag;
00753 
00754   return 0;
00755 }
00756 #endif
00757 
00758 //=============================================================================
00759 void Amesos_Mumps::PrintStatus() const
00760 {
00761 
00762   if (Comm().MyPID() != 0 ) return;
00763 
00764   //  The following lines are commented out to deal with bug #1887 - kss
00765 #ifndef IRIX64
00766   PrintLine();
00767   std::cout << "Amesos_Mumps : Matrix has " << Matrix().NumGlobalRows() << " rows"
00768        << " and " << Matrix().NumGlobalNonzeros() << " nonzeros" << std::endl;
00769   std::cout << "Amesos_Mumps : Nonzero elements per row = "
00770        << 1.0*Matrix().NumGlobalNonzeros()/Matrix().NumGlobalRows() << std::endl;
00771   std::cout << "Amesos_Mumps : Percentage of nonzero elements = "
00772        << 100.0*Matrix().NumGlobalNonzeros()/(pow(Matrix().NumGlobalRows(),2.0)) << std::endl;
00773   std::cout << "Amesos_Mumps : Use transpose = " << UseTranspose_ << std::endl;
00774 //  MatrixProperty_ is unused - see bug #2331 and bug #2332 in this file and bugzilla
00775   if (MatrixProperty_ == 0) std::cout << "Amesos_Mumps : Matrix is general unsymmetric" << std::endl;
00776   if (MatrixProperty_ == 2) std::cout << "Amesos_Mumps : Matrix is general symmetric" << std::endl;
00777   if (MatrixProperty_ == 1) std::cout << "Amesos_Mumps : Matrix is SPD" << std::endl;
00778   std::cout << "Amesos_Mumps : Available process(es) = " << Comm().NumProc() << std::endl;
00779   std::cout << "Amesos_Mumps : Using " << MaxProcs_ << " process(es)" << std::endl;
00780   
00781   std::cout << "Amesos_Mumps : Estimated FLOPS for elimination = "
00782        << MDS.RINFOG(1) << std::endl;
00783   std::cout << "Amesos_Mumps : Total FLOPS for assembly = "
00784        << MDS.RINFOG(2) << std::endl;
00785   std::cout << "Amesos_Mumps : Total FLOPS for elimination = "
00786        << MDS.RINFOG(3) << std::endl;
00787   
00788   std::cout << "Amesos_Mumps : Total real space to store the LU factors = "
00789        << MDS.INFOG(9) << std::endl;
00790   std::cout << "Amesos_Mumps : Total integer space to store the LU factors = "
00791        << MDS.INFOG(10) << std::endl;
00792   std::cout << "Amesos_Mumps : Total number of iterative steps refinement = "
00793        << MDS.INFOG(15) << std::endl;
00794   std::cout << "Amesos_Mumps : Estimated size of MUMPS internal data\n"
00795        << "Amesos_Mumps : for running factorization = "
00796        << MDS.INFOG(16) << " Mbytes" << std::endl;
00797   std::cout << "Amesos_Mumps : for running factorization = "
00798        << MDS.INFOG(17) << " Mbytes" << std::endl;
00799   std::cout << "Amesos_Mumps : Allocated during factorization = "
00800        << MDS.INFOG(19) << " Mbytes" << std::endl;
00801   PrintLine();
00802 #endif
00803 }
00804 
00805 //=============================================================================
00806 int Amesos_Mumps::CheckError() 
00807 {
00808   bool Wrong = ((MDS.INFOG(1) != 0) || (MDS.INFO(1) != 0))
00809                && (Comm().MyPID() < MaxProcs_);
00810   
00811   // an error occurred in MUMPS. Print out information and quit.
00812 
00813   if (Comm().MyPID() == 0 && Wrong) 
00814   {
00815     std::cerr << "Amesos_Mumps : ERROR" << std::endl;
00816     std::cerr << "Amesos_Mumps : INFOG(1) = " << MDS.INFOG(1) << std::endl;
00817     std::cerr << "Amesos_Mumps : INFOG(2) = " << MDS.INFOG(2) << std::endl;
00818   }
00819   
00820   if (MDS.INFO(1) != 0 && Wrong) 
00821   {
00822     std::cerr << "Amesos_Mumps : On process " << Comm().MyPID()
00823    << ", INFO(1) = " << MDS.INFO(1) << std::endl;
00824     std::cerr << "Amesos_Mumps : On process " << Comm().MyPID()
00825    << ", INFO(2) = " << MDS.INFO(2) << std::endl;
00826   }
00827 
00828   if (Wrong) 
00829       return 1 ;
00830   else
00831       return 0 ;
00832 }
00833 
00834 // ======================================================================
00835 void Amesos_Mumps::PrintTiming() const
00836 {
00837   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00838     return;
00839 
00840   double ConTime = GetTime(MtxConvTime_);
00841   double MatTime = GetTime(MtxRedistTime_);
00842   double VecTime = GetTime(VecRedistTime_);
00843   double SymTime = GetTime(SymFactTime_);
00844   double NumTime = GetTime(SymFactTime_);
00845   double SolTime = GetTime(SolveTime_);
00846 
00847   if (NumSymbolicFact_) SymTime /= NumSymbolicFact_;
00848   if (NumNumericFact_)  NumTime /= NumNumericFact_;
00849   if (NumSolve_)        SolTime /= NumSolve_;
00850 
00851   std::string p = "Amesos_Mumps : ";
00852   PrintLine();
00853 
00854   std::cout << p << "Time to convert matrix to MUMPS format = "
00855        << ConTime << " (s)" << std::endl;
00856   std::cout << p << "Time to redistribute matrix = "
00857        << MatTime << " (s)" << std::endl;
00858   std::cout << p << "Time to redistribute vectors = "
00859        << VecTime << " (s)" << std::endl;
00860   std::cout << p << "Number of symbolic factorizations = "
00861        << NumSymbolicFact_ << std::endl;
00862   std::cout << p << "Time for sym fact = "
00863        << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
00864   std::cout << p << "Number of numeric factorizations = "
00865        << NumNumericFact_ << std::endl;
00866   std::cout << p << "Time for num fact = "
00867        << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
00868   std::cout << p << "Number of solve phases = "
00869        << NumSolve_ << std::endl;
00870   std::cout << p << "Time for solve = "
00871        << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
00872 
00873   PrintLine();
00874 }
00875 
00876 // ===========================================================================
00877 Epetra_RowMatrix& Amesos_Mumps::Matrix() 
00878 {
00879   Epetra_RowMatrix* Matrix = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
00880   assert (Matrix != 0);
00881   return(*Matrix);
00882 }
00883 
00884 // ===========================================================================
00885 const Epetra_RowMatrix& Amesos_Mumps::Matrix() const
00886 {
00887   Epetra_RowMatrix* Matrix = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
00888   assert (Matrix != 0);
00889   return(*Matrix);
00890 }
00891 
00892 // ===========================================================================
00893 Epetra_Map& Amesos_Mumps::RedistrMap() 
00894 {
00895   assert (Comm().NumProc() != 1);
00896 
00897   if (RedistrMap_ == Teuchos::null) {
00898     int i = Matrix().NumGlobalRows() / MaxProcs_;
00899     if (Comm().MyPID() == 0)
00900       i += Matrix().NumGlobalRows() % MaxProcs_;
00901     else if (Comm().MyPID() >= MaxProcs_)
00902       i = 0;
00903 
00904     RedistrMap_ = rcp(new Epetra_Map(Matrix().NumGlobalRows(),i,0,Comm()));
00905     assert (RedistrMap_.get() != 0);
00906   }
00907   return(*RedistrMap_);
00908 }
00909 
00910 // ===========================================================================
00911 Epetra_Import& Amesos_Mumps::RedistrImporter()
00912 {
00913   assert (Comm().NumProc() != 1);
00914 
00915   if (RedistrImporter_ == null) {
00916     RedistrImporter_ = rcp(new Epetra_Import(RedistrMap(),Matrix().RowMatrixRowMap()));
00917     assert (RedistrImporter_.get() != 0);
00918   }
00919   return(*RedistrImporter_);
00920 }
00921 
00922 // ===========================================================================
00923 Epetra_RowMatrix& Amesos_Mumps::RedistrMatrix(const bool ImportMatrix)
00924 {
00925   if (Comm().NumProc() == 1) return(Matrix());
00926 
00927   if (ImportMatrix) RedistrMatrix_ = null;
00928 
00929   if (RedistrMatrix_ == null) 
00930   {
00931     RedistrMatrix_ = rcp(new Epetra_CrsMatrix(Copy,RedistrMap(),0));
00932     if (ImportMatrix) {
00933       int ierr = RedistrMatrix_->Import(Matrix(),RedistrImporter(),Insert);
00934       assert (ierr == 0);
00935       ierr = RedistrMatrix_->FillComplete();
00936       assert (ierr == 0);
00937     }
00938   }
00939 
00940   return(*RedistrMatrix_);
00941 }
00942 
00943 // ===========================================================================
00944 Epetra_Map& Amesos_Mumps::SerialMap()
00945 {
00946   if (SerialMap_ == null) 
00947   {
00948     int i = Matrix().NumGlobalRows();
00949     if (Comm().MyPID()) i = 0;
00950     SerialMap_ = rcp(new Epetra_Map(-1,i,0,Comm()));
00951     assert (SerialMap_ != null);
00952   }
00953   return(*SerialMap_);
00954 }
00955 
00956 // ===========================================================================
00957 Epetra_Import& Amesos_Mumps::SerialImporter()
00958 { 
00959   if (SerialImporter_ == null) {
00960     SerialImporter_ = rcp(new Epetra_Import(SerialMap(),Matrix().OperatorDomainMap()));
00961     assert (SerialImporter_ != null);
00962   }
00963   return(*SerialImporter_);
00964 }
00965 
00966 // ===========================================================================
00967 double* Amesos_Mumps::GetRINFO() 
00968 {
00969   return ( MDS.rinfo);
00970 }
00971 
00972 // ===========================================================================
00973 int* Amesos_Mumps::GetINFO() 
00974 {
00975   return (MDS.info);
00976 }
00977 
00978 // ===========================================================================
00979 double* Amesos_Mumps::GetRINFOG()
00980 {
00981   return (MDS.rinfog);
00982 }
00983 
00984 // ===========================================================================
00985 int* Amesos_Mumps::GetINFOG()
00986 {
00987   return (MDS.infog);
00988 }
00989 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines