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