IFPACK Development
Ifpack_SILU.cpp
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) 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 "Ifpack_ConfigDefs.h"
00030 #include "Ifpack_SILU.h"
00031 #ifdef HAVE_IFPACK_SUPERLU
00032 
00033 #include "Ifpack_CondestType.h"
00034 #include "Epetra_ConfigDefs.h"
00035 #include "Epetra_Comm.h"
00036 #include "Epetra_Map.h"
00037 #include "Epetra_RowMatrix.h"
00038 #include "Epetra_Vector.h"
00039 #include "Epetra_MultiVector.h"
00040 #include "Epetra_CrsGraph.h"
00041 #include "Epetra_CrsMatrix.h"
00042 #include "Teuchos_ParameterList.hpp"
00043 #include "Teuchos_RefCountPtr.hpp"
00044 
00045 // SuperLU includes
00046 extern "C" {int dfill_diag(int n, NCformat *Astore);}
00047 
00048 using Teuchos::RefCountPtr;
00049 using Teuchos::rcp;
00050 
00051 // Define this macro to see some timers for some of these functions
00052 #define ENABLE_IFPACK_ILU_TEUCHOS_TIMERS
00053 
00054 #ifdef ENABLE_IFPACK_ILU_TEUCHOS_TIMERS
00055 #  include "Teuchos_TimeMonitor.hpp"
00056 #endif
00057 
00058 //==============================================================================
00059 Ifpack_SILU::Ifpack_SILU(Epetra_RowMatrix* Matrix_in) :
00060   A_(rcp(Matrix_in,false)),
00061   Comm_(Matrix_in->Comm()),
00062   UseTranspose_(false),
00063   NumMyDiagonals_(0),
00064   DropTol_(1e-4),
00065   FillTol_(1e-2),
00066   FillFactor_(10.0),
00067   DropRule_(9), 
00068   Condest_(-1.0),
00069   IsInitialized_(false),
00070   IsComputed_(false),
00071   NumInitialize_(0),
00072   NumCompute_(0),
00073   NumApplyInverse_(0),
00074   InitializeTime_(0.0),
00075   ComputeTime_(0.0),
00076   ApplyInverseTime_(0.0),
00077   Time_(Comm()),
00078   etree_(0),
00079   perm_r_(0),
00080   perm_c_(0)
00081 {
00082   Teuchos::ParameterList List;
00083   SetParameters(List);
00084   SY_.Store=0;
00085 }
00086 
00087 
00088 
00089 //==============================================================================
00090 void Ifpack_SILU::Destroy()
00091 {
00092   if(IsInitialized_){
00093     // SuperLU cleanup
00094     StatFree(&stat_);
00095 
00096     Destroy_CompCol_Permuted(&SAc_);
00097     Destroy_SuperNode_Matrix(&SL_);
00098     Destroy_CompCol_Matrix(&SU_);
00099 
00100     // Make sure NOT to clean up Epetra's memory
00101     Destroy_SuperMatrix_Store(&SA_);
00102     if(SY_.Store) Destroy_SuperMatrix_Store(&SY_);
00103     SY_.Store=0;
00104 
00105     // Cleanup stuff I allocated
00106     delete [] etree_;etree_=0;
00107     delete [] perm_r_;perm_r_=0;
00108     delete [] perm_c_;perm_c_=0;  
00109   }
00110 }
00111 
00112 //==========================================================================
00113 int Ifpack_SILU::SetParameters(Teuchos::ParameterList& List)
00114 {
00115   DropTol_=List.get("fact: drop tolerance",DropTol_);
00116   FillTol_=List.get("fact: zero pivot threshold",FillTol_);
00117   FillFactor_=List.get("fact: maximum fill factor",FillFactor_);
00118   DropRule_=List.get("fact: silu drop rule",DropRule_);
00119 
00120   // set label
00121   sprintf(Label_, "IFPACK SILU (drop=%d, zpv=%f, ffact=%f, rthr=%f)",
00122       DropRule(),FillTol(),FillFactor(),DropTol());
00123   return(0);
00124 }
00125 
00126 //==========================================================================
00127 int Ifpack_SILU::Initialize() 
00128 {
00129 
00130 #ifdef ENABLE_IFPACK_ILU_TEUCHOS_TIMERS
00131   TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::Initialize");
00132 #endif
00133 
00134   Time_.ResetStartTime();
00135 
00136   // reset this object
00137   Destroy();
00138 
00139   IsInitialized_ = false;
00140 
00141   Epetra_CrsMatrix* CrsMatrix;
00142   CrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(&*A_);
00143 
00144   if(CrsMatrix && CrsMatrix->RowMap().SameAs(CrsMatrix->ColMap())){
00145     // Case #1: Use original matrix
00146     Aover_=rcp(CrsMatrix,false);
00147   }
00148   else if(CrsMatrix){
00149     // Case #2: Extract using CrsDataPointers
00150     int size = A_->MaxNumEntries();
00151     int N=A_->NumMyRows();
00152     Aover_ = rcp(new Epetra_CrsMatrix(Copy,A_->RowMatrixRowMap(), size));
00153     vector<int> Indices(size);
00154     vector<double> Values(size);
00155 
00156     int i,j,ct,*rowptr,*colind;
00157     double *values;
00158     CrsMatrix->ExtractCrsDataPointers(rowptr,colind,values);
00159  
00160     // Use the fact that EpetraCrsMatrices always number the off-processor columns *LAST*   
00161     for(i=0;i<N;i++){
00162       for(j=rowptr[i],ct=0;j<rowptr[i+1];j++){
00163     if(colind[j]<N){
00164       Indices[ct]=CrsMatrix->GCID(colind[j]);
00165       Values[ct]=values[j];
00166       ct++;
00167     }
00168       }
00169       Aover_->InsertGlobalValues(CrsMatrix->GRID(i),ct,&Values[0],&Indices[0]);
00170     }
00171     IFPACK_CHK_ERR(Aover_->FillComplete(CrsMatrix->RowMap(),CrsMatrix->RowMap()));  
00172   }
00173   else{
00174     // Case #3: Extract using copys
00175     int size = A_->MaxNumEntries();
00176     Aover_ = rcp(new Epetra_CrsMatrix(Copy,A_->RowMatrixRowMap(), size));
00177     if (Aover_.get() == 0) IFPACK_CHK_ERR(-5); // memory allocation error
00178 
00179     vector<int> Indices(size);
00180     vector<double> Values(size);
00181 
00182     // extract each row at-a-time, and insert it into
00183     // the graph, ignore all off-process entries
00184     for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
00185       int NumEntries;
00186       int GlobalRow = A_->RowMatrixRowMap().GID(i);
00187       IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries, 
00188                       &Values[0], &Indices[0]));
00189       // convert to global indices
00190       for (int j = 0 ; j < NumEntries ; ++j) {
00191     Indices[j] = A_->RowMatrixColMap().GID(Indices[j]); 
00192       }
00193       IFPACK_CHK_ERR(Aover_->InsertGlobalValues(GlobalRow,NumEntries,
00194                         &Values[0],&Indices[0]));
00195     }    
00196     IFPACK_CHK_ERR(Aover_->FillComplete(A_->RowMatrixRowMap(),A_->RowMatrixRowMap()));
00197   }
00198 
00199   // Finishing touches
00200   Aover_->OptimizeStorage();
00201   Graph_=rcp(const_cast<Epetra_CrsGraph*>(&Aover_->Graph()),false); 
00202 
00203   IsInitialized_ = true;
00204   NumInitialize_++;
00205   InitializeTime_ += Time_.ElapsedTime();
00206 
00207   return(0);
00208 }
00209 
00210 //==========================================================================
00211 int Ifpack_SILU::Compute() 
00212 {
00213 
00214 #ifdef ENABLE_IFPACK_ILU_TEUCHOS_TIMERS
00215   TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::Compute");
00216 #endif
00217 
00218   if (!IsInitialized()) 
00219     IFPACK_CHK_ERR(Initialize());
00220 
00221   Time_.ResetStartTime();
00222   IsComputed_ = false;
00223 
00224   // Initialize the SuperLU statistics & options variables.
00225   StatInit(&stat_);
00226   ilu_set_default_options(&options_);
00227   options_.ILU_DropTol=DropTol_;
00228   options_.ILU_FillTol=FillTol_;
00229   options_.ILU_DropRule=DropRule_;
00230   options_.ILU_FillFactor=FillFactor_;
00231 
00232   // Grab pointers from Aover_
00233   int *rowptr,*colind;
00234   double *values;
00235   Aover_->ExtractCrsDataPointers(rowptr,colind,values);
00236   int N=Aover_->NumMyRows();
00237 
00238   // Copy the data over to SuperLU land - mark as a transposed CompCol Matrix
00239   dCreate_CompCol_Matrix(&SA_,N,N,Aover_->NumMyNonzeros(),
00240              values,colind,rowptr,SLU_NC,SLU_D,SLU_GE);
00241 
00242   // Fill any zeros on the diagonal
00243   // Commented out for now
00244   dfill_diag(N, (NCformat*)SA_.Store);
00245 
00246   // Allocate SLU memory
00247   etree_=new int [N];
00248   perm_r_=new int[N];
00249   perm_c_=new int[N];
00250 
00251   // Get column permutation
00252   int permc_spec=options_.ColPerm;
00253   if ( permc_spec != MY_PERMC && options_.Fact == DOFACT )
00254     get_perm_c(permc_spec,&SA_,perm_c_);
00255 
00256   // Preorder by column permutation
00257   sp_preorder(&options_, &SA_, perm_c_, etree_, &SAc_);
00258 
00259   // Call the factorization
00260   int panel_size = sp_ienv(1);
00261   int relax      = sp_ienv(2);
00262   int info=0;
00263   dgsitrf(&options_,&SAc_,relax,panel_size,etree_,NULL,0,perm_c_,perm_r_,&SL_,&SU_,&stat_,&info);
00264   if(info<0) IFPACK_CHK_ERR(info);
00265 
00266   IsComputed_ = true;
00267   NumCompute_++;
00268   ComputeTime_ += Time_.ElapsedTime();
00269   return 0;
00270 }
00271 
00272 //=============================================================================
00273 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00274 int Ifpack_SILU::Solve(bool Trans, const Epetra_MultiVector& X, 
00275                       Epetra_MultiVector& Y) const 
00276 {
00277 
00278 #ifdef ENABLE_IFPACK_ILU_TEUCHOS_TIMERS
00279   TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::ApplyInverse - Solve");
00280 #endif
00281 
00282   if (!IsComputed())
00283     IFPACK_CHK_ERR(-3);
00284   int nrhs=X.NumVectors();
00285   int N=A_->NumMyRows();
00286 
00287   // Copy X over to Y
00288   Y=X;
00289 
00290   // Move to SuperLU land
00291   // NTS: Need to do epetra-style realloc-if-nrhs-changes thing
00292   if(SY_.Store) Destroy_SuperMatrix_Store(&SY_);
00293   SY_.Store=0;
00294   dCreate_Dense_Matrix(&SY_,N,nrhs,Y[0],N,SLU_DN,SLU_D,SLU_GE);
00295 
00296   int info;
00297   dgstrs(TRANS,&SL_,&SU_,perm_c_,perm_r_,&SY_,&stat_,&info);
00298   if(!info) IFPACK_CHK_ERR(info);
00299 
00300   return(info);
00301 }
00302 
00303 //=============================================================================
00304 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00305 int Ifpack_SILU::Multiply(bool Trans, const Epetra_MultiVector& X, 
00306                 Epetra_MultiVector& Y) const 
00307 {
00308 
00309   if (!IsComputed())
00310     IFPACK_CHK_ERR(-1);
00311 
00312   return(0);
00313 }
00314 
00315 //=============================================================================
00316 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00317 int Ifpack_SILU::ApplyInverse(const Epetra_MultiVector& X, 
00318                              Epetra_MultiVector& Y) const
00319 {
00320 
00321 #ifdef ENABLE_IFPACK_ILU_TEUCHOS_TIMERS
00322   TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::ApplyInverse");
00323 #endif
00324 
00325   if (!IsComputed())
00326     IFPACK_CHK_ERR(-3);
00327 
00328   if (X.NumVectors() != Y.NumVectors())
00329     IFPACK_CHK_ERR(-2);
00330 
00331   Time_.ResetStartTime();
00332 
00333   // AztecOO gives X and Y pointing to the same memory location,
00334   // need to create an auxiliary vector, Xcopy
00335   Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
00336   if (X.Pointers()[0] == Y.Pointers()[0])
00337     Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
00338   else
00339     Xcopy = Teuchos::rcp( &X, false );
00340 
00341   IFPACK_CHK_ERR(Solve(Ifpack_SILU::UseTranspose(), *Xcopy, Y));
00342 
00343   ++NumApplyInverse_;
00344   ApplyInverseTime_ += Time_.ElapsedTime();
00345 
00346   return(0);
00347 
00348 }
00349 
00350 //=============================================================================
00351 double Ifpack_SILU::Condest(const Ifpack_CondestType CT, 
00352                            const int MaxIters, const double Tol,
00353                               Epetra_RowMatrix* Matrix_in)
00354 {
00355 
00356 #ifdef ENABLE_IFPACK_ILU_TEUCHOS_TIMERS
00357   TEUCHOS_FUNC_TIME_MONITOR("Ifpack_SILU::Condest");
00358 #endif
00359 
00360   if (!IsComputed()) // cannot compute right now
00361     return(-1.0);
00362 
00363   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00364 
00365   return(Condest_);
00366 }
00367 
00368 //=============================================================================
00369 std::ostream&
00370 Ifpack_SILU::Print(std::ostream& os) const
00371 {
00372   if (!Comm().MyPID()) {
00373     os << endl;
00374     os << "================================================================================" << endl;
00375     os << "Ifpack_SILU: " << Label() << endl << endl;
00376     os << "Dropping rule      = "<< DropRule() << endl;
00377     os << "Zero pivot thresh  = "<< FillTol() << endl;
00378     os << "Max fill factor    = "<< FillFactor() << endl;
00379     os << "Drop tolerance     = "<< DropTol() << endl;
00380     os << "Condition number estimate = " << Condest() << endl;
00381     os << "Global number of rows     = " << A_->NumGlobalRows() << endl;
00382     if (IsComputed_) {
00383       // Internal SuperLU info
00384       int fnnz=0;
00385       if(SL_.Store) fnnz+=((SCformat*)SL_.Store)->nnz;
00386       if(SU_.Store) fnnz+=((NCformat*)SU_.Store)->nnz;
00387       int lufill=fnnz - A_->NumMyRows();
00388       os << "No. of nonzeros in L+U    = "<< lufill<<endl;
00389       os << "Fill ratio: nnz(F)/nnz(A) = "<<(double)lufill / (double)A_->NumMyNonzeros()<<endl;
00390     }
00391     os << endl;
00392     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00393     os << "-----           -------   --------------       ------------     --------" << endl;
00394     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
00395        << "  " << std::setw(15) << InitializeTime() 
00396        << "              0.0              0.0" << endl;
00397     os << "Compute()       "   << std::setw(5) << NumCompute() 
00398        << "  " << std::setw(15) << ComputeTime()
00399        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
00400     if (ComputeTime() != 0.0) 
00401       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
00402     else
00403       os << "  " << std::setw(15) << 0.0 << endl;
00404     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00405        << "  " << std::setw(15) << ApplyInverseTime()
00406        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
00407     if (ApplyInverseTime() != 0.0) 
00408       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
00409     else
00410       os << "  " << std::setw(15) << 0.0 << endl;
00411     os << "================================================================================" << endl;
00412     os << endl;
00413   }
00414 
00415   return(os);
00416 }
00417 
00418 #endif
 All Classes Files Functions Variables Enumerations Friends