Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_IKLU.cpp
Go to the documentation of this file.
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 
00030 #include "Ifpack_ConfigDefs.h"
00031 #include "Ifpack_Preconditioner.h"
00032 #include "Ifpack_IKLU.h"
00033 #include "Ifpack_Condest.h"
00034 #include "Ifpack_Utils.h"
00035 #include "Ifpack_HashTable.h"
00036 #include "Epetra_SerialComm.h"
00037 #include "Epetra_Comm.h"
00038 #include "Epetra_Map.h"
00039 #include "Epetra_RowMatrix.h"
00040 #include "Epetra_CrsMatrix.h"
00041 #include "Epetra_Vector.h"
00042 #include "Epetra_MultiVector.h"
00043 #include "Epetra_Util.h"
00044 #include "Teuchos_ParameterList.hpp"
00045 #include "Teuchos_RefCountPtr.hpp"
00046 #include <functional>
00047 
00048 using namespace Teuchos;
00049 
00050 //==============================================================================
00051 // FIXME: allocate Comm_ and Time_ the first Initialize() call
00052 Ifpack_IKLU::Ifpack_IKLU(const Epetra_RowMatrix* A) :
00053   A_(*A),
00054   Comm_(A->Comm()),
00055   Condest_(-1.0),
00056   Relax_(0.),
00057   Athresh_(0.0),
00058   Rthresh_(1.0),
00059   LevelOfFill_(1.0),
00060   DropTolerance_(1e-12),
00061   IsInitialized_(false),
00062   IsComputed_(false),
00063   UseTranspose_(false),
00064   NumMyRows_(-1),
00065   NumInitialize_(0),
00066   NumCompute_(0),
00067   NumApplyInverse_(0),
00068   InitializeTime_(0.0),
00069   ComputeTime_(0.0),
00070   ApplyInverseTime_(0.0),
00071   ComputeFlops_(0.0),
00072   ApplyInverseFlops_(0.0),
00073   Time_(Comm()),
00074   GlobalNonzeros_(0),
00075   csrA_(0),
00076   cssS_(0),
00077   csrnN_(0)
00078 {
00079   // do nothing here..
00080 }
00081 
00082 //==============================================================================
00083 Ifpack_IKLU::~Ifpack_IKLU()
00084 {
00085   Destroy();
00086 }
00087 
00088 //==============================================================================
00089 void Ifpack_IKLU::Destroy()
00090 {
00091   IsInitialized_ = false;
00092   IsComputed_ = false;
00093   if (csrA_)
00094     csr_spfree( csrA_ );
00095   if (cssS_)
00096     csr_sfree( cssS_ );
00097   if (csrnN_)
00098     csr_nfree( csrnN_ );
00099 }
00100 
00101 //==========================================================================
00102 int Ifpack_IKLU::SetParameters(Teuchos::ParameterList& List)
00103 {
00104   try 
00105   {
00106     LevelOfFill_ = List.get<double>("fact: ilut level-of-fill", LevelOfFill());
00107     if (LevelOfFill_ <= 0.0)
00108       IFPACK_CHK_ERR(-2); // must be greater than 0.0
00109 
00110     Athresh_ = List.get<double>("fact: absolute threshold", Athresh_);
00111     Rthresh_ = List.get<double>("fact: relative threshold", Rthresh_);
00112     Relax_ = List.get<double>("fact: relax value", Relax_);
00113     DropTolerance_ = List.get<double>("fact: drop tolerance", DropTolerance_);
00114 
00115     Label_ = "IFPACK IKLU (fill=" + Ifpack_toString(LevelOfFill())
00116       + ", relax=" + Ifpack_toString(RelaxValue())
00117       + ", athr=" + Ifpack_toString(AbsoluteThreshold())
00118       + ", rthr=" + Ifpack_toString(RelativeThreshold())
00119       + ", droptol=" + Ifpack_toString(DropTolerance())
00120       + ")";
00121     return(0);
00122   }
00123   catch (...)
00124   {
00125     cerr << "Caught an exception while parsing the parameter list" << endl;
00126     cerr << "This typically means that a parameter was set with the" << endl;
00127     cerr << "wrong type (for example, int instead of double). " << endl;
00128     cerr << "please check the documentation for the type required by each parameer." << endl;
00129     IFPACK_CHK_ERR(-1);
00130   }
00131 
00132   return(0);
00133 }
00134 
00135 //==========================================================================
00136 int Ifpack_IKLU::Initialize()
00137 {
00138   // delete previously allocated factorization
00139   Destroy();
00140 
00141   Time_.ResetStartTime();
00142 
00143   if (A_.Comm().NumProc() != 1) {
00144     cout << " There are too many processors !!! " << endl;
00145     cerr << "Ifpack_IKLU can be used with Comm().NumProc() == 1" << endl;
00146     cerr << "only. This class is a subdomain solver for Ifpack_AdditiveSchwarz," << endl;
00147     cerr << "and it is currently not meant to be used otherwise." << endl;
00148     exit(EXIT_FAILURE);
00149   }
00150   
00151   // check dimensions of input matrix only in serial
00152   if (Comm().NumProc() == 1 && Matrix().NumMyRows() != Matrix().NumMyCols())
00153     IFPACK_CHK_ERR(-2);
00154     
00155   NumMyRows_ = Matrix().NumMyRows();
00156   NumMyNonzeros_ = Matrix().NumMyNonzeros();
00157 
00158   int RowNnz, Length = Matrix().MaxNumEntries();
00159   vector<int>    RowIndices(Length);
00160   vector<double> RowValues(Length);
00161 
00162   //cout << "Processor " << Comm().MyPID() << " owns " << NumMyRows_ << " rows and has " << NumMyNonzeros_ << " nonzeros " << endl;
00163   // get general symbolic structure of the matrix
00164   csrA_ = csr_spalloc( NumMyRows_, NumMyRows_, NumMyNonzeros_, 1, 0 );
00165 
00166   // copy the symbolic structure into csrA_
00167   int count = 0;
00168   csrA_->p[0] = 0;
00169   for (int i = 0; i < NumMyRows_; ++i ) {
00170 
00171     IFPACK_CHK_ERR(A_.ExtractMyRowCopy(i,Length,RowNnz,
00172                &RowValues[0],&RowIndices[0]));
00173     for (int j = 0 ; j < RowNnz ; ++j) {
00174       csrA_->j[count++] = RowIndices[j];
00175       //cout << "Row = " << i << ", Column = " << RowIndices[j] << ", Value = " << RowValues[j] << endl;
00176     }
00177     csrA_->p[i+1] = csrA_->p[i] + RowNnz;
00178   }
00179 
00180   // Perform symbolic analysis on the current matrix structure
00181   int order = 1;
00182   cssS_ = csr_sqr( order, csrA_ );
00183   for (int i = 0; i < NumMyRows_; ++i ) {
00184      cout << "AMD Perm (from inside KLU) [" << i << "] = " << cssS_->q[i] << endl;
00185   }
00186 
00187   // nothing else to do here
00188   IsInitialized_ = true;
00189   ++NumInitialize_;
00190   InitializeTime_ += Time_.ElapsedTime();
00191 
00192   return(0);
00193 }
00194 
00195 //==========================================================================
00196 class Ifpack_AbsComp 
00197 {
00198  public:
00199   inline bool operator()(const double& x, const double& y) 
00200   {
00201     return(IFPACK_ABS(x) > IFPACK_ABS(y));
00202   }
00203 };
00204 
00205 //==========================================================================
00206 
00207 int Ifpack_IKLU::Compute() 
00208 {
00209   if (!IsInitialized()) 
00210     IFPACK_CHK_ERR(Initialize());
00211 
00212   Time_.ResetStartTime();
00213   IsComputed_ = false;
00214 
00215   NumMyRows_ = A_.NumMyRows();
00216   int Length = A_.MaxNumEntries();
00217 
00218   bool distributed = (Comm().NumProc() > 1)?true:false;
00219   if (distributed)
00220   {
00221     SerialComm_ = rcp(new Epetra_SerialComm);
00222     SerialMap_ = rcp(new Epetra_Map(NumMyRows_, 0, *SerialComm_));
00223     assert (SerialComm_.get() != 0);
00224     assert (SerialMap_.get() != 0);
00225   }
00226   else
00227     SerialMap_ = rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()), false);
00228 
00229   int RowNnz;
00230   vector<int>    RowIndices(Length);
00231   vector<double> RowValues(Length);
00232 
00233   // copy the values from A_ into csrA_
00234   int count = 0;
00235   for (int i = 0; i < NumMyRows_; ++i ) {
00236 
00237     IFPACK_CHK_ERR(A_.ExtractMyRowCopy(i,Length,RowNnz,
00238                &RowValues[0],&RowIndices[0]));
00239     // make sure each row has the same number of nonzeros
00240     if (RowNnz != (csrA_->p[i+1]-csrA_->p[i])) {
00241       cout << "The number of nonzeros for this row does not math the expected number of nonzeros!!!" << endl;
00242     }
00243     for (int j = 0 ; j < RowNnz ; ++j) {
00244       
00245       csrA_->x[count++] = RowValues[j];
00246       //cout << "Row = " << i << ", Column = " << RowIndices[j] << ", Value = " << RowValues[j] << endl;
00247     }
00248   }
00249   
00250   // compute the lu factors
00251   double tol = 0.1;
00252   csrnN_ = csr_lu( &*csrA_, &*cssS_, tol );
00253 
00254   // Create L and U as a view of the information stored in csrnN_->L and csrnN_->U
00255   csr* L_tmp = csrnN_->L;
00256   csr* U_tmp = csrnN_->U;
00257   vector<int> numEntriesL( NumMyRows_ ), numEntriesU( NumMyRows_ );
00258   for (int i=0; i < NumMyRows_; ++i) {
00259     numEntriesL[i] = ( L_tmp->p[i+1] - L_tmp->p[i] );
00260     numEntriesU[i] = ( U_tmp->p[i+1] - U_tmp->p[i] );
00261   }
00262   L_ = rcp(new Epetra_CrsMatrix(View, *SerialMap_, &numEntriesL[0]));
00263   U_ = rcp(new Epetra_CrsMatrix(View, *SerialMap_, &numEntriesU[0]));
00264 
00265   // Insert the values into L and U
00266   for (int i=0; i < NumMyRows_; ++i) {
00267     L_->InsertGlobalValues( i, numEntriesL[i], &(L_tmp->x[L_tmp->p[i]]), &(L_tmp->j[L_tmp->p[i]]) );
00268     U_->InsertGlobalValues( i, numEntriesU[i], &(U_tmp->x[U_tmp->p[i]]), &(U_tmp->j[U_tmp->p[i]]) );
00269   }
00270 
00271   IFPACK_CHK_ERR(L_->FillComplete());
00272   IFPACK_CHK_ERR(U_->FillComplete());
00273 
00274   int MyNonzeros = L_->NumGlobalNonzeros() + U_->NumGlobalNonzeros();
00275   Comm().SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
00276 
00277   IsComputed_ = true;
00278 
00279   ++NumCompute_;
00280   ComputeTime_ += Time_.ElapsedTime();
00281 
00282   return(0);
00283 
00284 }
00285   
00286 //=============================================================================
00287 int Ifpack_IKLU::ApplyInverse(const Epetra_MultiVector& X, 
00288            Epetra_MultiVector& Y) const
00289 {
00290   if (!IsComputed())
00291     IFPACK_CHK_ERR(-2); // compute preconditioner first
00292 
00293   if (X.NumVectors() != Y.NumVectors()) 
00294     IFPACK_CHK_ERR(-3); // Return error: X and Y not the same size
00295 
00296   Time_.ResetStartTime();
00297 
00298   // NOTE: L_ and U_ are based on SerialMap_, while Xcopy is based
00299   // on A.Map()... which are in general different. However, Solve()
00300   // does not seem to care... which is fine with me.
00301   //
00302   // AztecOO gives X and Y pointing to the same memory location,
00303   // need to create an auxiliary vector, Xcopy and apply permutation.
00304   vector<int> invq( NumMyRows_ );
00305 
00306   for (int i=0; i<NumMyRows_; ++i ) {
00307     csrnN_->perm[ csrnN_->pinv[i] ] = i;
00308     invq[ cssS_->q[i] ] = i;
00309   }
00310 
00311   Teuchos::RefCountPtr<Epetra_MultiVector> Xcopy = Teuchos::rcp( new Epetra_MultiVector(X.Map(),X.NumVectors()), false );
00312   Teuchos::RefCountPtr<Epetra_MultiVector> Ytemp = Teuchos::rcp( new Epetra_MultiVector(Y.Map(),Y.NumVectors()) );
00313 
00314   for (int i=0; i<NumMyRows_; ++i) {
00315     for (int j=0; j<X.NumVectors(); ++j) {
00316       Xcopy->ReplaceMyValue( invq[i], j, (*X(j))[i] );
00317     }
00318   }
00319 
00320   if (!UseTranspose_)
00321   {
00322     // solves LU Y = X 
00323     IFPACK_CHK_ERR(L_->Solve(false,false,false,*Xcopy,*Ytemp));
00324     IFPACK_CHK_ERR(U_->Solve(true,false,false,*Ytemp,*Ytemp));
00325   }
00326   else
00327   {
00328     // solves U(trans) L(trans) Y = X
00329     IFPACK_CHK_ERR(U_->Solve(true,true,false,*Xcopy,*Ytemp));
00330     IFPACK_CHK_ERR(L_->Solve(false,true,false,*Ytemp,*Ytemp));
00331   }
00332 
00333   // Reverse the permutation.
00334   for (int i=0; i<NumMyRows_; ++i) {
00335     for (int j=0; j<Y.NumVectors(); ++j) {
00336       Y.ReplaceMyValue( csrnN_->perm[i], j, (*(*Ytemp)(j))[i] );
00337     }
00338   }
00339 
00340   ++NumApplyInverse_;
00341   ApplyInverseFlops_ += X.NumVectors() * 2 * GlobalNonzeros_;
00342   ApplyInverseTime_ += Time_.ElapsedTime();
00343 
00344   return(0);
00345 
00346 }
00347 //=============================================================================
00348 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00349 int Ifpack_IKLU::Apply(const Epetra_MultiVector& X, 
00350           Epetra_MultiVector& Y) const 
00351 {
00352 
00353   return(-98);
00354 }
00355 
00356 //=============================================================================
00357 double Ifpack_IKLU::Condest(const Ifpack_CondestType CT, 
00358                             const int MaxIters, const double Tol,
00359           Epetra_RowMatrix* Matrix_in)
00360 {
00361   if (!IsComputed()) // cannot compute right now
00362     return(-1.0);
00363 
00364   // NOTE: this is computing the *local* condest
00365   if (Condest_ == -1.0)
00366     Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00367 
00368   return(Condest_);
00369 }
00370 
00371 //=============================================================================
00372 std::ostream&
00373 Ifpack_IKLU::Print(std::ostream& os) const
00374 {
00375   if (!Comm().MyPID()) {
00376     os << endl;
00377     os << "================================================================================" << endl;
00378     os << "Ifpack_IKLU: " << Label() << endl << endl;
00379     os << "Level-of-fill      = " << LevelOfFill() << endl;
00380     os << "Absolute threshold = " << AbsoluteThreshold() << endl;
00381     os << "Relative threshold = " << RelativeThreshold() << endl;
00382     os << "Relax value        = " << RelaxValue() << endl;
00383     os << "Condition number estimate       = " << Condest() << endl;
00384     os << "Global number of rows           = " << A_.NumGlobalRows() << endl;
00385     if (IsComputed_) {
00386       os << "Number of nonzeros in A         = " << A_.NumGlobalNonzeros() << endl;
00387       os << "Number of nonzeros in L + U     = " << NumGlobalNonzeros() 
00388          << " ( = " << 100.0 * NumGlobalNonzeros() / A_.NumGlobalNonzeros() 
00389          << " % of A)" << endl;
00390       os << "nonzeros / rows                 = " 
00391         << 1.0 * NumGlobalNonzeros() / U_->NumGlobalRows() << endl;
00392     }
00393     os << endl;
00394     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00395     os << "-----           -------   --------------       ------------     --------" << endl;
00396     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
00397        << "  " << std::setw(15) << InitializeTime() 
00398        << "               0.0            0.0" << endl;
00399     os << "Compute()       "   << std::setw(5) << NumCompute() 
00400        << "  " << std::setw(15) << ComputeTime()
00401        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops(); 
00402     if (ComputeTime() != 0.0)
00403       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
00404     else
00405       os << "  " << std::setw(15) << 0.0 << endl;
00406     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00407        << "  " << std::setw(15) << ApplyInverseTime()
00408        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
00409     if (ApplyInverseTime() != 0.0) 
00410       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
00411     else
00412       os << "  " << std::setw(15) << 0.0 << endl;
00413     os << "================================================================================" << endl;
00414     os << endl;
00415   }
00416 
00417   return(os);
00418 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines