Amesos Package Browser (Single Doxygen Collection) Development
Amesos_Superludist.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 // TO DO: use Stat structure to print statistics???
00030 // allow users to specify usermap ???
00031 
00032 #include "Amesos_Superludist.h"
00033 #include "Epetra_Map.h"
00034 #include "Epetra_Import.h"
00035 #include "Epetra_CrsMatrix.h"
00036 #include "Epetra_Util.h"
00037 // #include "CrsMatrixTranspose.h"
00038 #include "superlu_ddefs.h"
00039 #include "supermatrix.h"
00040 //  SuperLU defines Reduce to be a macro in util.h, this conflicts with Reduce() in Epetra_MultiVector.h
00041 #undef Reduce
00042 
00043 class Amesos_Superlu_Pimpl {
00044 public:
00045   //   Teuchos::RCP<klu_symbolic> Symbolic_ ;
00046   //   Teuchos::RCP<klu_numeric> Numeric_ ;
00047   fact_t FactOption_; 
00048   //  Here are the structures used by Superlu
00049   SuperMatrix SuperluA_;
00050   ScalePermstruct_t ScalePermstruct_;
00051   LUstruct_t LUstruct_;
00052   SOLVEstruct_t SOLVEstruct_; 
00054   gridinfo_t grid_;
00056   superlu_options_t options_;
00057   
00058 } ;
00059 
00060 
00061 // ====================================================================== 
00062 int Superludist_NumProcRows( int NumProcs ) {
00063 #ifdef TFLOP
00064   //  Else, parameter 6 of DTRSV CTNLU is incorrect 
00065   return 1;
00066 #else
00067   int i;
00068   int NumProcRows ;
00069   for ( i = 1; i*i <= NumProcs; i++ ) 
00070     ;
00071   bool done = false ;
00072   for ( NumProcRows = i-1 ; done == false ; ) {
00073     int NumCols = NumProcs / NumProcRows ; 
00074     if ( NumProcRows * NumCols == NumProcs ) 
00075       done = true; 
00076     else 
00077       NumProcRows-- ; 
00078   }
00079   return NumProcRows;
00080 #endif
00081 }
00082 
00083 // ====================================================================== 
00084 int SetNPRowAndCol(const int MaxProcesses, int& nprow, int& npcol)
00085 {
00086   nprow = Superludist_NumProcRows(MaxProcesses);
00087   if (nprow < 1 ) nprow = 1;
00088   npcol = MaxProcesses / nprow;
00089   
00090   if( nprow <=0 || npcol <= 0 || MaxProcesses <=0 ) {
00091     std::cerr << "Amesos_Superludist: wrong value for MaxProcs ("
00092    << MaxProcesses << "), or nprow (" << nprow 
00093    << ") or npcol (" << npcol << ")" << std::endl;
00094     AMESOS_CHK_ERR(-1);
00095   }
00096   return(0);
00097 }
00098 
00099 //=============================================================================
00100 Amesos_Superludist::Amesos_Superludist(const Epetra_LinearProblem &prob) :
00101   PrivateSuperluData_( rcp( new Amesos_Superlu_Pimpl() ) ),
00102   Problem_(&prob),
00103   RowMatrixA_(0), 
00104   GridCreated_(0), 
00105   FactorizationDone_(0), 
00106   FactorizationOK_(false),    
00107   NumGlobalRows_(0), 
00108   nprow_(0),       
00109   npcol_(0),       
00110   PrintNonzeros_(false),
00111   ColPerm_("NOT SET"),
00112   RowPerm_("NOT SET"),
00113   perm_c_(0),         
00114   perm_r_(0),         
00115   IterRefine_("NOT SET"),
00116   ReplaceTinyPivot_(true),
00117   MtxConvTime_(-1),
00118   MtxRedistTime_(-1),
00119   VecRedistTime_(-1),
00120   NumFactTime_(-1),
00121   SolveTime_(-1),
00122   OverheadTime_(-1)
00123 {
00124   Redistribute_ = true;
00125   AddZeroToDiag_ = false;
00126   PrivateSuperluData_->FactOption_ = SamePattern_SameRowPerm ;
00127   ReuseSymbolic_ = false ; 
00128 
00129   MaxProcesses_ = - 1; 
00130   Equil_ = true;
00131   ColPerm_ = "MMD_AT_PLUS_A";
00132   perm_c_ = 0;
00133   RowPerm_ = "LargeDiag";
00134   perm_r_ = 0;
00135   IterRefine_ = "DOUBLE";
00136   ReplaceTinyPivot_ = true;
00137   
00138   PrintNonzeros_ = false;
00139 
00140   ComputeTrueResidual_ = false;
00141   ComputeVectorNorms_ = false;
00142   PrintStatus_ = false ; 
00143   PrintTiming_ = false ; 
00144   
00145   Teuchos::ParameterList ParamList;
00146   SetParameters(ParamList); 
00147 }
00148 
00149 //=============================================================================
00150 Amesos_Superludist::~Amesos_Superludist(void) 
00151 {
00152   if (PrintTiming_) PrintTiming();
00153   if (PrintStatus_) PrintStatus();
00154 
00155   if ( FactorizationDone_ ) {
00156     SUPERLU_FREE( PrivateSuperluData_->SuperluA_.Store );
00157     ScalePermstructFree(&PrivateSuperluData_->ScalePermstruct_);
00158     Destroy_LU(NumGlobalRows_, &PrivateSuperluData_->grid_, &PrivateSuperluData_->LUstruct_);
00159     LUstructFree(&PrivateSuperluData_->LUstruct_);
00160     if ( PrivateSuperluData_->options_.SolveInitialized ) {
00161       dSolveFinalize(&PrivateSuperluData_->options_, &PrivateSuperluData_->SOLVEstruct_ ) ; 
00162     }
00163   }
00164   if ( GridCreated_ ) {
00165 #ifndef IFPACK_SUBCOMM_CODE
00166     // TODO: Cleanup of MPI comms and types in SuperLU causes an MPI error
00167     // when used in conjuction with Ifpack_SubdomainFilter. This hack WILL
00168     // cause memory leaks when Amesos_Superlu is created and destroyed
00169     // multiple times per execution (only with IFPACK_SUBCOMM_CODE defined)
00170     superlu_gridexit(&PrivateSuperluData_->grid_);
00171 #endif
00172   }
00173 }
00174 
00175 // ====================================================================== 
00176 int Amesos_Superludist::SetParameters( Teuchos::ParameterList &ParameterList ) 
00177 {
00178   // retrive general parameters
00179 
00180   SetStatusParameters( ParameterList );
00181 
00182   SetControlParameters( ParameterList );
00183 
00184   if (ParameterList.isParameter("Redistribute"))
00185     Redistribute_ = ParameterList.get<bool>("Redistribute");  
00186 
00187   // parameters for Superludist only
00188 
00189   if (ParameterList.isSublist("Superludist") ) 
00190   {
00191     const Teuchos::ParameterList& SuperludistParams = 
00192       ParameterList.sublist("Superludist") ;
00193 
00194     if( SuperludistParams.isParameter("ReuseSymbolic") )
00195       ReuseSymbolic_ = SuperludistParams.get<bool>("ReuseSymbolic");
00196     std::string FactOption = "NotSet";
00197 
00198     if( SuperludistParams.isParameter("Fact") )
00199       FactOption = SuperludistParams.get<std::string>("Fact");
00200 
00201     if( FactOption == "SamePattern_SameRowPerm" ) PrivateSuperluData_->FactOption_ = SamePattern_SameRowPerm;
00202     else if( FactOption == "SamePattern" ) PrivateSuperluData_->FactOption_ = SamePattern;
00203     else if ( FactOption != "NotSet" ) 
00204       AMESOS_CHK_ERR(-2); // input not valid
00205 
00206     if (SuperludistParams.isParameter("Equil"))
00207       Equil_ = SuperludistParams.get<bool>("Equil");
00208 
00209     if (SuperludistParams.isParameter("ColPerm"))
00210       ColPerm_ = SuperludistParams.get<std::string>("ColPerm");
00211 
00212     if (ColPerm_ == "MY_PERMC")
00213       if( SuperludistParams.isParameter("perm_c"))
00214         perm_c_ = SuperludistParams.get<int*>("perm_c");
00215 
00216     if (SuperludistParams.isParameter("RowPerm"))
00217       RowPerm_ = SuperludistParams.get<std::string>("RowPerm");
00218     if( RowPerm_ == "MY_PERMR" ) {
00219       if (SuperludistParams.isParameter("perm_r"))
00220         perm_r_ = SuperludistParams.get<int*>("perm_r");
00221     }
00222 
00223     if (SuperludistParams.isParameter("IterRefine"))
00224       IterRefine_ = SuperludistParams.get<std::string>("IterRefine");
00225 
00226     if (SuperludistParams.isParameter("ReplaceTinyPivot"))
00227       ReplaceTinyPivot_ = SuperludistParams.get<bool>("ReplaceTinyPivot");
00228 
00229     if (SuperludistParams.isParameter("PrintNonzeros"))
00230       PrintNonzeros_ = SuperludistParams.get<bool>("PrintNonzeros");
00231   }
00232 
00233   return(0);
00234 }
00235 
00236 // ====================================================================== 
00237 // Tasks of this method:
00238 // 1) To set the required number of processes;
00239 // 2) To set nprow_ and npcol_
00240 // 3) To create a linear distribution (map) with elements only in the
00241 //    active processes
00242 // 4) To redistribute the matrix from the original to the linear map.
00243 // ====================================================================== 
00244 int Amesos_Superludist::RedistributeA()
00245 {
00246   ResetTimer(0);
00247   
00248   if (NumGlobalRows_ != RowMatrixA_->NumGlobalRows())
00249     AMESOS_CHK_ERR(-1); // something has changed
00250 
00251   int iam = Comm().MyPID();
00252 
00253   SetMaxProcesses(MaxProcesses_, *RowMatrixA_);
00254   SetNPRowAndCol(MaxProcesses_, nprow_, npcol_);
00255 
00256   int m_per_p = NumGlobalRows_ / MaxProcesses_ ;
00257   int remainder = NumGlobalRows_ - ( m_per_p * MaxProcesses_ );
00258   int MyFirstElement = iam * m_per_p + EPETRA_MIN( iam, remainder );
00259   int MyFirstNonElement = (iam+1) * m_per_p + EPETRA_MIN( iam+1, remainder );
00260   int NumExpectedElements = MyFirstNonElement - MyFirstElement ; 
00261 
00262   if (iam >= MaxProcesses_)
00263     NumExpectedElements = 0; 
00264 
00265   if (NumGlobalRows_ !=  RowMatrixA_->NumGlobalRows())
00266     AMESOS_CHK_ERR(-1);
00267 
00268   const Epetra_Map &OriginalMap = RowMatrixA_->RowMatrixRowMap();
00269 
00270   UniformMap_       = rcp(new Epetra_Map(NumGlobalRows_, NumExpectedElements, 
00271                                          0, Comm()));
00272   CrsUniformMatrix_ = rcp(new Epetra_CrsMatrix(Copy, *UniformMap_, 0));
00273   UniformMatrix_    = rcp(&CrsUniformMatrix(), false);
00274   Importer_         = rcp(new Epetra_Import(*UniformMap_, OriginalMap));
00275 
00276   CrsUniformMatrix_->Import(*RowMatrixA_, Importer(), Add); 
00277   //  int TracebackMode = Comm().GetTracebackMode();
00278   //  UniformMatrix_->SetTracebackMode(0);
00279   if (AddZeroToDiag_ ) { 
00280      //
00281      //  Add 0.0 to each diagonal entry to avoid empty diagonal entries;
00282      //
00283       int OriginalTracebackMode = CrsUniformMatrix_->GetTracebackMode() ; 
00284       CrsUniformMatrix_->SetTracebackMode( EPETRA_MIN( OriginalTracebackMode, 0) ) ;
00285      double zero = 0.0;
00286      for ( int i = 0 ; i < UniformMap_->NumGlobalElements(); i++ ) 
00287        if ( CrsUniformMatrix_->LRID(i) >= 0 ) 
00288   CrsUniformMatrix_->InsertGlobalValues( i, 1, &zero, &i ) ;
00289      CrsUniformMatrix_->SetTracebackMode( OriginalTracebackMode ) ; 
00290    }
00291   //  UniformMatrix_->SetTracebackMode(TracebackMode);
00292 
00293 
00294   CrsUniformMatrix_->FillComplete(); 
00295 
00296   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
00297 
00298 
00299   
00300   return(0);
00301 }
00302 
00303 
00304 // ====================================================================== 
00305 int Amesos_Superludist::Factor()
00306 {
00307   ResetTimer(1); // for "overhead"
00308 
00309   // FIXME????
00310   //  For now, if you change the shape of a matrix, you need to 
00311   //  create a new Amesos instance.
00312   //  
00313   //
00314   if (NumGlobalRows_ != 0 && NumGlobalRows_ != RowMatrixA_->NumGlobalRows())
00315     AMESOS_CHK_ERR(-5);
00316 
00317   NumGlobalRows_ = RowMatrixA_->NumGlobalRows() ; 
00318 
00319   if (Comm().NumProc() == 1)
00320     Redistribute_ = false;
00321 
00322   // Set the matrix and grid shapes. Time is tracked within
00323   // the RedistributeA() function
00324 
00325   if (Redistribute_)
00326     RedistributeA() ; 
00327   else 
00328   {
00329     if (Comm().NumProc() == 1)
00330     {
00331       nprow_ = 1;
00332       npcol_ = 1;
00333     }
00334     else 
00335     {
00336       if (!(RowMatrixA_->RowMatrixRowMap().LinearMap())) 
00337         AMESOS_CHK_ERR(-2);
00338       SetNPRowAndCol(Comm().NumProc(), nprow_, npcol_);
00339     }
00340 
00341     UniformMatrix_ = rcp(RowMatrixA_, false);
00342   }
00343 
00344   //  Extract Ai_, Ap_ and Aval_ from UniformMatrix_
00345 
00346   ResetTimer(0);
00347 
00348   int MyActualFirstElement = UniformMatrix().RowMatrixRowMap().MinMyGID() ; 
00349   int NumMyElements = UniformMatrix().NumMyRows() ; 
00350   int nnz_loc = UniformMatrix().NumMyNonzeros() ;
00351   Ap_.resize( NumMyElements+1 );
00352   Ai_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ; 
00353   Aval_.resize( EPETRA_MAX( NumMyElements, nnz_loc) ) ; 
00354   
00355   int NzThisRow ;
00356   int Ai_index = 0 ; 
00357   int MyRow;
00358   //int num_my_cols = UniformMatrix().NumMyCols() ; 
00359   double *RowValues;
00360   int *ColIndices;
00361   int MaxNumEntries_ = UniformMatrix().MaxNumEntries();
00362   std::vector<double> RowValuesV_(MaxNumEntries_);
00363   std::vector<int>    ColIndicesV_(MaxNumEntries_);
00364 
00365   Global_Columns_ = UniformMatrix().RowMatrixColMap().MyGlobalElements();
00366 
00367   const Epetra_CrsMatrix *SuperluCrs = dynamic_cast<const Epetra_CrsMatrix *>(&UniformMatrix());
00368 
00369   int ierr;
00370 
00371   for (MyRow = 0; MyRow < NumMyElements ; MyRow++) 
00372   {
00373     if (SuperluCrs != 0) 
00374     {
00375       ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues, 
00376                                       ColIndices);
00377 
00378     }
00379     else {
00380       ierr = UniformMatrix().ExtractMyRowCopy(MyRow, MaxNumEntries_, NzThisRow,
00381                                               &RowValuesV_[0], &ColIndicesV_[0]);
00382       RowValues =  &RowValuesV_[0];
00383       ColIndices = &ColIndicesV_[0];
00384     }
00385 
00386     AMESOS_CHK_ERR(ierr);
00387 
00388     // MS // Added on 15-Mar-05
00389     if (AddToDiag_ != 0.0 || AddZeroToDiag_) {
00390       for (int i = 0 ; i < NzThisRow ; ++i) {
00391         if (ColIndices[i] == MyRow) {
00392           RowValues[i] += AddToDiag_;
00393           break;
00394         }
00395       }
00396     }
00397 
00398     Ap_[MyRow] = Ai_index ; 
00399     for ( int j = 0; j < NzThisRow; j++ ) { 
00400       Ai_[Ai_index] = Global_Columns_[ColIndices[j]] ; 
00401       Aval_[Ai_index] = RowValues[j] ;
00402       Ai_index++;
00403     }
00404   }
00405   assert( NumMyElements == MyRow );
00406   Ap_[ NumMyElements ] = Ai_index ; 
00407 
00408   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00409 
00410   //
00411   //  Setup Superlu's grid 
00412   //
00413   const Epetra_MpiComm & comm1 = dynamic_cast<const Epetra_MpiComm &> (Comm());
00414 
00415   if ( ! GridCreated_ ) {
00416     // NOTE: nprow_ and npcol_ cannot be changed by the user
00417     GridCreated_ = true;
00418     superlu_gridinit(comm1.Comm(), nprow_, npcol_, &PrivateSuperluData_->grid_);
00419   }
00420 
00421   if ( FactorizationDone_ ) {
00422     SUPERLU_FREE( PrivateSuperluData_->SuperluA_.Store );
00423     ScalePermstructFree(&PrivateSuperluData_->ScalePermstruct_);
00424     Destroy_LU(NumGlobalRows_, &PrivateSuperluData_->grid_, &PrivateSuperluData_->LUstruct_);
00425     LUstructFree(&PrivateSuperluData_->LUstruct_);
00426     if ( PrivateSuperluData_->options_.SolveInitialized ) {
00427       dSolveFinalize(&PrivateSuperluData_->options_, &PrivateSuperluData_->SOLVEstruct_ ) ; 
00428     }
00429   }
00430 
00431   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
00432   ResetTimer(0);
00433 
00434   //
00435   //  Only those processes in the grid participate from here on
00436   //
00437   if (Comm().MyPID() < nprow_ * npcol_) {
00438     //
00439     //  Set up Superlu's data structures
00440     //
00441     set_default_options_dist(&PrivateSuperluData_->options_);
00442 
00443     dCreate_CompRowLoc_Matrix_dist( &PrivateSuperluData_->SuperluA_, NumGlobalRows_, NumGlobalRows_, 
00444             nnz_loc, NumMyElements, MyActualFirstElement,
00445             &Aval_[0], &Ai_[0], &Ap_[0], 
00446             SLU_NR_loc, SLU_D, SLU_GE );
00447 
00448     FactorizationDone_ = true;   // i.e. clean up Superlu data structures in the destructor
00449 
00450     ScalePermstructInit(NumGlobalRows_, NumGlobalRows_, &PrivateSuperluData_->ScalePermstruct_);
00451     LUstructInit(NumGlobalRows_, NumGlobalRows_, &PrivateSuperluData_->LUstruct_);
00452 
00453     // stick options from ParameterList to options_ structure
00454     // Here we follow the same order of the SuperLU_dist 2.0 manual (pag 55/56)
00455     
00456     assert( PrivateSuperluData_->options_.Fact == DOFACT );  
00457     PrivateSuperluData_->options_.Fact = DOFACT ;       
00458 
00459     if( Equil_ ) PrivateSuperluData_->options_.Equil = (yes_no_t)YES;
00460     else         PrivateSuperluData_->options_.Equil = NO;
00461 
00462     if( ColPerm_ == "NATURAL" ) PrivateSuperluData_->options_.ColPerm = NATURAL;
00463     else if( ColPerm_ == "MMD_AT_PLUS_A" ) PrivateSuperluData_->options_.ColPerm = MMD_AT_PLUS_A;
00464     else if( ColPerm_ == "MMD_ATA" ) PrivateSuperluData_->options_.ColPerm = MMD_ATA;
00465     //    else if( ColPerm_ == "COLAMD" ) PrivateSuperluData_->options_.ColPerm = COLAMD;     // COLAMD no longer supported in Superludist, as of July 2005
00466     else if( ColPerm_ == "MY_PERMC" ) {
00467       PrivateSuperluData_->options_.ColPerm = MY_PERMC;
00468       PrivateSuperluData_->ScalePermstruct_.perm_c = perm_c_;
00469     }
00470 
00471     if( RowPerm_ == "NATURAL" ) PrivateSuperluData_->options_.RowPerm = (rowperm_t)NATURAL;
00472     if( RowPerm_ == "LargeDiag" ) PrivateSuperluData_->options_.RowPerm = LargeDiag;
00473     else if( ColPerm_ == "MY_PERMR" ) {
00474       PrivateSuperluData_->options_.RowPerm = MY_PERMR;
00475       PrivateSuperluData_->ScalePermstruct_.perm_r = perm_r_;
00476     }
00477 
00478     if( ReplaceTinyPivot_ ) PrivateSuperluData_->options_.ReplaceTinyPivot = (yes_no_t)YES;
00479     else                    PrivateSuperluData_->options_.ReplaceTinyPivot = (yes_no_t)NO;
00480 
00481     if( IterRefine_ == "NO" ) PrivateSuperluData_->options_.IterRefine = (IterRefine_t)NO;
00482     else if( IterRefine_ == "DOUBLE" ) PrivateSuperluData_->options_.IterRefine = DOUBLE;
00483     else if( IterRefine_ == "EXTRA" ) PrivateSuperluData_->options_.IterRefine = EXTRA;
00484 
00485     //  Without the following two lines, SuperLU_DIST cannot be made
00486     //  quiet.
00487 
00488     if (PrintNonzeros_) PrivateSuperluData_->options_.PrintStat = (yes_no_t)YES;
00489     else                PrivateSuperluData_->options_.PrintStat = (yes_no_t)NO;
00490     
00491     SuperLUStat_t stat;
00492     PStatInit(&stat);    /* Initialize the statistics variables. */
00493 
00494     //
00495     //  Factor A using Superludsit (via a call to pdgssvx)
00496     //
00497     int info ;
00498     double berr ;    //  Should be untouched
00499     double xValues;  //  Should be untouched
00500     int nrhs = 0 ;   //  Prevents forward and back solves
00501     int ldx = NumGlobalRows_;     //  Should be untouched
00502 
00503     pdgssvx(&PrivateSuperluData_->options_, &PrivateSuperluData_->SuperluA_, &PrivateSuperluData_->ScalePermstruct_, &xValues, ldx, nrhs, &PrivateSuperluData_->grid_,
00504       &PrivateSuperluData_->LUstruct_, &PrivateSuperluData_->SOLVEstruct_, &berr, &stat, &info);
00505 
00506     if ( PrivateSuperluData_->options_.SolveInitialized ) {
00507       dSolveFinalize(&PrivateSuperluData_->options_, &PrivateSuperluData_->SOLVEstruct_ ) ; 
00508     }
00509     AMESOS_CHK_ERR(info);
00510 
00511     PStatFree(&stat);
00512   }
00513 
00514   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
00515 
00516   return 0;
00517 }
00518 
00519 // ====================================================================== 
00520 //   Refactor - Refactor the matrix 
00521 //
00522 //     Preconditions:
00523 //       The non-zero pattern of the matrix must not have changed since the 
00524 //         previous call to Factor().  Refactor ensures that each process owns 
00525 //         the same number of columns that it did on the previous call to Factor()
00526 //         and returns -4 if a discrepancy is found.  However, that check does not
00527 //         guarantee that no change was made to the non-zero structure of the matrix.
00528 //       No call to SetParameters should be made between the call to Factor()
00529 //         and the call to Refactor().  If the user does not call SetParameters, 
00530 //         as they need never do, they are safe on this.
00531 //
00532 //     Postconditions:
00533 //       The matrix specified by Problem_->Operator() will have been redistributed,
00534 //         converted to the form needed by Superludist and factored.
00535 //       Ai_, Aval_ 
00536 //       SuperluA_
00537 //       SuperLU internal data structures reflecting the LU factorization
00538 //         ScalePermstruct_
00539 //         LUstructInit_
00540 //       
00541 //     Performance notes:
00542 //       Refactor does not allocate or de-allocate memory.
00543 //         
00544 //     Return codes:
00545 //       -4 if we detect a change to the non-zero structure of the matrix.
00546 //
00547 int Amesos_Superludist::ReFactor( ) 
00548 {
00549   ResetTimer(0);
00550   ResetTimer(1);
00551 
00552   //
00553   //  Update Ai_ and Aval_ (while double checking Ap_)
00554   //
00555   if (Redistribute_)  
00556     if(CrsUniformMatrix().Import(*RowMatrixA_, Importer(), Insert)) 
00557       AMESOS_CHK_ERR(-4); 
00558 
00559   MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_, 0);
00560   ResetTimer(0);
00561 
00562   const Epetra_CrsMatrix *SuperluCrs = dynamic_cast<const Epetra_CrsMatrix *>(&UniformMatrix());
00563 
00564   double *RowValues;
00565   int *ColIndices;
00566   int MaxNumEntries_ = UniformMatrix().MaxNumEntries();
00567   int NumMyElements  = UniformMatrix().NumMyRows() ; 
00568   std::vector<int> ColIndicesV_(MaxNumEntries_);
00569   std::vector<double> RowValuesV_(MaxNumEntries_);
00570 
00571   int NzThisRow ;
00572   int Ai_index = 0 ; 
00573   int MyRow;
00574   int ierr;
00575 
00576   for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
00577     if ( SuperluCrs != 0 ) {
00578       ierr = SuperluCrs->ExtractMyRowView(MyRow, NzThisRow, RowValues, 
00579                                           ColIndices);
00580     }
00581     else {
00582       ierr = UniformMatrix().ExtractMyRowCopy( MyRow, MaxNumEntries_,
00583                                               NzThisRow, &RowValuesV_[0], 
00584                                               &ColIndicesV_[0]);
00585       RowValues =  &RowValuesV_[0];
00586       ColIndices = &ColIndicesV_[0];
00587     }
00588 
00589     AMESOS_CHK_ERR(ierr);
00590 
00591     if ( Ap_[MyRow] != Ai_index ) AMESOS_CHK_ERR(-4);
00592     for ( int j = 0; j < NzThisRow; j++ ) { 
00593       //  pdgssvx alters Ai_, so we have to set it again.
00594       Ai_[Ai_index] = Global_Columns_[ColIndices[j]];
00595       Aval_[Ai_index] = RowValues[j] ;  
00596       Ai_index++;
00597     }
00598   }
00599   if( Ap_[ NumMyElements ] != Ai_index ) AMESOS_CHK_ERR(-4); 
00600 
00601   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
00602   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00603   ResetTimer(0);
00604 
00605   if (Comm().MyPID() < nprow_ * npcol_) {
00606 
00607 
00608     //  If we reuse the same options, the code fails on multiprocess runs
00609     set_default_options_dist(&PrivateSuperluData_->options_);   
00610 
00611     if (PrintNonzeros_) PrivateSuperluData_->options_.PrintStat = (yes_no_t)YES;
00612     else                PrivateSuperluData_->options_.PrintStat = (yes_no_t)NO;
00613     
00614 
00615     PrivateSuperluData_->options_.Fact = PrivateSuperluData_->FactOption_;
00616     SuperLUStat_t stat;
00617     PStatInit(&stat);    /* Initialize the statistics variables. */
00618     int info ;
00619     double berr ;    //  Should be untouched
00620     double xValues;  //  Should be untouched
00621     int nrhs = 0 ;   //  Prevents forward and back solves
00622     int ldx = NumGlobalRows_;     //  Should be untouched
00623     pdgssvx(&PrivateSuperluData_->options_, &PrivateSuperluData_->SuperluA_, &PrivateSuperluData_->ScalePermstruct_, &xValues, ldx, nrhs, &PrivateSuperluData_->grid_,
00624             &PrivateSuperluData_->LUstruct_, &PrivateSuperluData_->SOLVEstruct_, &berr, &stat, &info);
00625     PStatFree(&stat);
00626     AMESOS_CHK_ERR( info ) ;
00627   } 
00628 
00629   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
00630 
00631   return 0;
00632 }
00633 
00634 // ====================================================================== 
00635 bool Amesos_Superludist::MatrixShapeOK() const 
00636 { 
00637   if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() != 
00638       GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints()) 
00639     return(false);
00640   else
00641     return(true);
00642 }
00643 
00644 // ====================================================================== 
00645 int Amesos_Superludist::SymbolicFactorization() 
00646 {
00647   FactorizationOK_ = false ; 
00648 
00649   return(0);
00650 }
00651 
00652 // ====================================================================== 
00653 int Amesos_Superludist::NumericFactorization() 
00654 {
00655   RowMatrixA_ = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
00656   if (RowMatrixA_ == 0) 
00657     AMESOS_CHK_ERR(-1); // Linear problem does not contain Epetra_RowMatrix
00658 
00659   // reset factorization
00660   //  FactorizationOK_ = false; 
00661 
00662   if (!MatrixShapeOK())
00663     AMESOS_CHK_ERR(-1); // matrix not square
00664 
00665   CreateTimer(Comm(), 2);
00666 
00667   if (FactorizationOK_ && ReuseSymbolic_)
00668     ReFactor();
00669   else
00670     Factor();
00671 
00672   FactorizationOK_ = true;   
00673 
00674   NumNumericFact_++;
00675   
00676   return(0);
00677 }
00678 
00679 // ====================================================================== 
00680 // We proceed as follows:
00681 // - perform numeric factorization if not done;
00682 // - extract solution and right-hand side;
00683 // - redistribute them if necessary by creating additional
00684 //   working vectors, or use vecX and vecB otherwise;     
00685 // - call SuperLU_DIST's solve routine;                    
00686 // - re-ship data if necessary.                             
00687 // ====================================================================== 
00688 int Amesos_Superludist::Solve() 
00689 {
00690   if (!FactorizationOK_)
00691     AMESOS_CHK_ERR(NumericFactorization());
00692 
00693   ResetTimer(1);
00694 
00695   Epetra_MultiVector* vecX = Problem_->GetLHS(); 
00696   Epetra_MultiVector* vecB = Problem_->GetRHS(); 
00697 
00698   if (vecX == 0 || vecB == 0)
00699     AMESOS_CHK_ERR(-1);
00700 
00701   int nrhs = vecX->NumVectors() ; 
00702   if (vecB->NumVectors() != nrhs)
00703     AMESOS_CHK_ERR(-1);
00704 
00705   double *values;
00706   int ldx;
00707 
00708   RCP<Epetra_MultiVector> vec_uni;
00709 
00710   if (Redistribute_) 
00711   {
00712     vec_uni = Teuchos::rcp(new Epetra_MultiVector(*UniformMap_, nrhs)); 
00713     ResetTimer(0);
00714     vec_uni->Import(*vecB, Importer(), Insert);
00715     VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
00716   } 
00717   else 
00718   {
00719     vecX->Update(1.0, *vecB, 0.0);
00720     vec_uni = Teuchos::rcp(vecX, false); 
00721   }
00722 
00723   //int NumMyElements = vec_uni->MyLength(); 
00724   AMESOS_CHK_ERR(vec_uni->ExtractView(&values, &ldx)); 
00725 
00726   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00727 
00728   ResetTimer(0);
00729   
00730   /* Bail out if I do not belong in the grid. */
00731   if (Comm().MyPID() < nprow_ * npcol_) 
00732   {
00733     int info ;
00734     std::vector<double>berr(nrhs);
00735     SuperLUStat_t stat;
00736     PStatInit(&stat);    /* Initialize the statistics variables. */
00737     
00738     if (!GridCreated_ || !FactorizationDone_)
00739       AMESOS_CHK_ERR(-1); // internal error
00740     PrivateSuperluData_->options_.Fact = FACTORED ;       
00741 
00742     //bool BlockSolve = true ; 
00743     pdgssvx(&PrivateSuperluData_->options_, &PrivateSuperluData_->SuperluA_, &PrivateSuperluData_->ScalePermstruct_, &values[0], ldx, 
00744             nrhs, &PrivateSuperluData_->grid_, &PrivateSuperluData_->LUstruct_, &PrivateSuperluData_->SOLVEstruct_, &berr[0], 
00745             &stat, &info);
00746     AMESOS_CHK_ERR(info);
00747     
00748     PStatFree(&stat);
00749   }
00750 
00751   SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
00752 
00753   ResetTimer(1);
00754   
00755   if (Redistribute_) 
00756   {
00757     ResetTimer(0);
00758     vecX->Export(*vec_uni, Importer(), Insert);
00759     VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
00760   }
00761 
00762   if (ComputeTrueResidual_)
00763     ComputeTrueResidual(*RowMatrixA_, *vecX, *vecB, UseTranspose(), 
00764                         "Amesos_Superludist");
00765 
00766   if (ComputeVectorNorms_)
00767     ComputeVectorNorms(*vecX, *vecB, "Amesos_Superludist");
00768 
00769   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00770   NumSolve_++;
00771   
00772   return(0);
00773 }
00774 
00775 // ====================================================================== 
00776 void Amesos_Superludist::PrintStatus() const
00777 {
00778   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00779     return;
00780   
00781   std::string p = "Amesos_Superludist : ";
00782   int NNZ = RowMatrixA_->NumGlobalNonzeros();
00783 
00784   PrintLine();
00785 
00786   std::cout << p << "Matrix has " << NumGlobalRows_ << " rows"
00787        << " and " << NNZ << " nonzeros" << std::endl;
00788   std::cout << p << "Nonzero elements per row = "
00789        << 1.0 * NNZ / NumGlobalRows_ << std::endl;
00790   std::cout << p << "Percentage of nonzero elements = "
00791        << 100.0 * NNZ /pow(NumGlobalRows_, 2.0) << std::endl;
00792   std::cout << p << "Use transpose = " << UseTranspose() << std::endl;
00793   std::cout << p << "Redistribute = " << Redistribute_ << std::endl;
00794   std::cout << p << "# available processes = " << Comm().NumProc() << std::endl;
00795   std::cout << p << "# processes used in computation = " << nprow_ * npcol_
00796        << " ( = " << nprow_ << "x" << npcol_ << ")" << std::endl;
00797   std::cout << p << "Equil = " << Equil_ << std::endl;
00798   std::cout << p << "ColPerm = " << ColPerm_ << std::endl;
00799   std::cout << p << "RowPerm = " << RowPerm_ << std::endl;
00800   std::cout << p << "IterRefine = " << IterRefine_ << std::endl;
00801   std::cout << p << "ReplaceTinyPivot = " << ReplaceTinyPivot_ << std::endl;
00802   std::cout << p << "AddZeroToDiag = " << AddZeroToDiag_ << std::endl;
00803   std::cout << p << "Redistribute = " << Redistribute_ << std::endl;
00804   
00805   PrintLine();
00806 
00807   return;
00808 }
00809 
00810 // ====================================================================== 
00811 void Amesos_Superludist::PrintTiming() const
00812 {
00813   if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
00814     return;
00815 
00816   double ConTime = GetTime(MtxConvTime_);
00817   double MatTime = GetTime(MtxRedistTime_);
00818   double VecTime = GetTime(VecRedistTime_);
00819   double NumTime = GetTime(NumFactTime_);
00820   double SolTime = GetTime(SolveTime_);
00821   double OveTime = GetTime(OverheadTime_);
00822 
00823   if (NumNumericFact_)
00824     NumTime /= NumNumericFact_;
00825 
00826   if (NumSolve_)
00827     SolTime /= NumSolve_;
00828 
00829   std::string p = "Amesos_Superludist : ";
00830   PrintLine();
00831 
00832   std::cout << p << "Time to convert matrix to Superludist format = "
00833        << ConTime << " (s)" << std::endl;
00834   std::cout << p << "Time to redistribute matrix = "
00835        << MatTime << " (s)" << std::endl;
00836   std::cout << p << "Time to redistribute vectors = "
00837        << VecTime << " (s)" << std::endl;
00838   std::cout << p << "Number of symbolic factorizations = "
00839        << NumSymbolicFact_ << std::endl;
00840   std::cout << p << "Time for sym fact = 0.0 (s), avg = 0.0 (s)" << std::endl;
00841   std::cout << p << "Number of numeric factorizations = "
00842        << NumNumericFact_ << std::endl;
00843   std::cout << p << "Time for num fact = "
00844        << NumTime * NumNumericFact_ << " (s), avg = " << NumTime << " (s)" << std::endl;
00845   std::cout << p << "Number of solve phases = "
00846        << NumSolve_ << std::endl;
00847   std::cout << p << "Time for solve = "
00848        << SolTime * NumSolve_ << " (s), avg = " << SolTime << " (s)" << std::endl;
00849 
00850   double tt = NumTime * NumNumericFact_ + SolTime * NumSolve_;
00851   if (tt != 0)
00852   {
00853     std::cout << p << "Total time spent in Amesos = " << tt << " (s) " << std::endl;
00854     std::cout << p << "Total time spent in the Amesos interface = " << OveTime << " (s)" << std::endl;
00855     std::cout << p << "(the above time does not include SuperLU_DIST time)" << std::endl;
00856     std::cout << p << "Amesos interface time / total time = " << OveTime / tt << std::endl;
00857   }
00858 
00859   PrintLine();
00860 
00861   return;
00862 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines