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