Ifpack_IKLU.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 
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 
00184   // nothing else to do here
00185   IsInitialized_ = true;
00186   ++NumInitialize_;
00187   InitializeTime_ += Time_.ElapsedTime();
00188 
00189   return(0);
00190 }
00191 
00192 //==========================================================================
00193 class Ifpack_AbsComp 
00194 {
00195  public:
00196   inline bool operator()(const double& x, const double& y) 
00197   {
00198     return(IFPACK_ABS(x) > IFPACK_ABS(y));
00199   }
00200 };
00201 
00202 //==========================================================================
00203 
00204 int Ifpack_IKLU::Compute() 
00205 {
00206   if (!IsInitialized()) 
00207     IFPACK_CHK_ERR(Initialize());
00208 
00209   Time_.ResetStartTime();
00210   IsComputed_ = false;
00211 
00212   NumMyRows_ = A_.NumMyRows();
00213   int Length = A_.MaxNumEntries();
00214 
00215   bool distributed = (Comm().NumProc() > 1)?true:false;
00216   if (distributed)
00217   {
00218     SerialComm_ = rcp(new Epetra_SerialComm);
00219     SerialMap_ = rcp(new Epetra_Map(NumMyRows_, 0, *SerialComm_));
00220     assert (SerialComm_.get() != 0);
00221     assert (SerialMap_.get() != 0);
00222   }
00223   else
00224     SerialMap_ = rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()), false);
00225 
00226   int RowNnz;
00227   vector<int>    RowIndices(Length);
00228   vector<double> RowValues(Length);
00229 
00230   // copy the values from A_ into csrA_
00231   int count = 0;
00232   for (int i = 0; i < NumMyRows_; ++i ) {
00233 
00234     IFPACK_CHK_ERR(A_.ExtractMyRowCopy(i,Length,RowNnz,
00235                        &RowValues[0],&RowIndices[0]));
00236     // make sure each row has the same number of nonzeros
00237     if (RowNnz != (csrA_->p[i+1]-csrA_->p[i])) {
00238       cout << "The number of nonzeros for this row does not math the expected number of nonzeros!!!" << endl;
00239     }
00240     for (int j = 0 ; j < RowNnz ; ++j) {
00241       
00242       csrA_->x[count++] = RowValues[j];
00243       //cout << "Row = " << i << ", Column = " << RowIndices[j] << ", Value = " << RowValues[j] << endl;
00244     }
00245   }
00246   
00247   // compute the lu factors
00248   double tol = 0.1;
00249   csrnN_ = csr_lu( &*csrA_, &*cssS_, tol );
00250 
00251   // Create L and U as a view of the information stored in csrnN_->L and csrnN_->U
00252   csr* L_tmp = csrnN_->L;
00253   csr* U_tmp = csrnN_->U;
00254   vector<int> numEntriesL( NumMyRows_ ), numEntriesU( NumMyRows_ );
00255   for (int i=0; i < NumMyRows_; ++i) {
00256     numEntriesL[i] = ( L_tmp->p[i+1] - L_tmp->p[i] );
00257     numEntriesU[i] = ( U_tmp->p[i+1] - U_tmp->p[i] );
00258   }
00259   L_ = rcp(new Epetra_CrsMatrix(View, *SerialMap_, &numEntriesL[0]));
00260   U_ = rcp(new Epetra_CrsMatrix(View, *SerialMap_, &numEntriesU[0]));
00261 
00262   // Insert the values into L and U
00263   for (int i=0; i < NumMyRows_; ++i) {
00264     L_->InsertGlobalValues( i, numEntriesL[i], &(L_tmp->x[L_tmp->p[i]]), &(L_tmp->j[L_tmp->p[i]]) );
00265     U_->InsertGlobalValues( i, numEntriesU[i], &(U_tmp->x[U_tmp->p[i]]), &(U_tmp->j[U_tmp->p[i]]) );
00266   }
00267 
00268   IFPACK_CHK_ERR(L_->FillComplete());
00269   IFPACK_CHK_ERR(U_->FillComplete());
00270 
00271   int MyNonzeros = L_->NumGlobalNonzeros() + U_->NumGlobalNonzeros();
00272   Comm().SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
00273 
00274   IsComputed_ = true;
00275 
00276   ++NumCompute_;
00277   ComputeTime_ += Time_.ElapsedTime();
00278 
00279   return(0);
00280 
00281 }
00282   
00283 //=============================================================================
00284 int Ifpack_IKLU::ApplyInverse(const Epetra_MultiVector& X, 
00285                  Epetra_MultiVector& Y) const
00286 {
00287   if (!IsComputed())
00288     IFPACK_CHK_ERR(-2); // compute preconditioner first
00289 
00290   if (X.NumVectors() != Y.NumVectors()) 
00291     IFPACK_CHK_ERR(-3); // Return error: X and Y not the same size
00292 
00293   Time_.ResetStartTime();
00294 
00295   // NOTE: L_ and U_ are based on SerialMap_, while Xcopy is based
00296   // on A.Map()... which are in general different. However, Solve()
00297   // does not seem to care... which is fine with me.
00298   //
00299   // AztecOO gives X and Y pointing to the same memory location,
00300   // need to create an auxiliary vector, Xcopy and apply permutation.
00301   vector<int> invq( NumMyRows_ );
00302 
00303   for (int i=0; i<NumMyRows_; ++i ) {
00304     csrnN_->perm[ csrnN_->pinv[i] ] = i;
00305     invq[ cssS_->q[i] ] = i;
00306   }
00307 
00308   Teuchos::RefCountPtr<Epetra_MultiVector> Xcopy = Teuchos::rcp( new Epetra_MultiVector(X.Map(),X.NumVectors()), false );
00309   Teuchos::RefCountPtr<Epetra_MultiVector> Ytemp = Teuchos::rcp( new Epetra_MultiVector(Y.Map(),Y.NumVectors()) );
00310 
00311   for (int i=0; i<NumMyRows_; ++i) {
00312     for (int j=0; j<X.NumVectors(); ++j) {
00313       Xcopy->ReplaceMyValue( invq[i], j, (*X(j))[i] );
00314     }
00315   }
00316 
00317   if (!UseTranspose_)
00318   {
00319     // solves LU Y = X 
00320     IFPACK_CHK_ERR(L_->Solve(false,false,false,*Xcopy,*Ytemp));
00321     IFPACK_CHK_ERR(U_->Solve(true,false,false,*Ytemp,*Ytemp));
00322   }
00323   else
00324   {
00325     // solves U(trans) L(trans) Y = X
00326     IFPACK_CHK_ERR(U_->Solve(true,true,false,*Xcopy,*Ytemp));
00327     IFPACK_CHK_ERR(L_->Solve(false,true,false,*Ytemp,*Ytemp));
00328   }
00329 
00330   // Reverse the permutation.
00331   for (int i=0; i<NumMyRows_; ++i) {
00332     for (int j=0; j<Y.NumVectors(); ++j) {
00333       Y.ReplaceMyValue( csrnN_->perm[i], j, (*(*Ytemp)(j))[i] );
00334     }
00335   }
00336 
00337   ++NumApplyInverse_;
00338   ApplyInverseFlops_ += X.NumVectors() * 2 * GlobalNonzeros_;
00339   ApplyInverseTime_ += Time_.ElapsedTime();
00340 
00341   return(0);
00342 
00343 }
00344 //=============================================================================
00345 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00346 int Ifpack_IKLU::Apply(const Epetra_MultiVector& X, 
00347               Epetra_MultiVector& Y) const 
00348 {
00349 
00350   return(-98);
00351 }
00352 
00353 //=============================================================================
00354 double Ifpack_IKLU::Condest(const Ifpack_CondestType CT, 
00355                             const int MaxIters, const double Tol,
00356                 Epetra_RowMatrix* Matrix_in)
00357 {
00358   if (!IsComputed()) // cannot compute right now
00359     return(-1.0);
00360 
00361   // NOTE: this is computing the *local* condest
00362   if (Condest_ == -1.0)
00363     Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
00364 
00365   return(Condest_);
00366 }
00367 
00368 //=============================================================================
00369 std::ostream&
00370 Ifpack_IKLU::Print(std::ostream& os) const
00371 {
00372   if (!Comm().MyPID()) {
00373     os << endl;
00374     os << "================================================================================" << endl;
00375     os << "Ifpack_IKLU: " << Label() << endl << endl;
00376     os << "Level-of-fill      = " << LevelOfFill() << endl;
00377     os << "Absolute threshold = " << AbsoluteThreshold() << endl;
00378     os << "Relative threshold = " << RelativeThreshold() << endl;
00379     os << "Relax value        = " << RelaxValue() << endl;
00380     os << "Condition number estimate       = " << Condest() << endl;
00381     os << "Global number of rows           = " << A_.NumGlobalRows() << endl;
00382     if (IsComputed_) {
00383       os << "Number of nonzeros in A         = " << A_.NumGlobalNonzeros() << endl;
00384       os << "Number of nonzeros in L + U     = " << NumGlobalNonzeros() 
00385          << " ( = " << 100.0 * NumGlobalNonzeros() / A_.NumGlobalNonzeros() 
00386          << " % of A)" << endl;
00387       os << "nonzeros / rows                 = " 
00388         << 1.0 * NumGlobalNonzeros() / U_->NumGlobalRows() << endl;
00389     }
00390     os << endl;
00391     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00392     os << "-----           -------   --------------       ------------     --------" << endl;
00393     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
00394        << "  " << std::setw(15) << InitializeTime() 
00395        << "               0.0            0.0" << endl;
00396     os << "Compute()       "   << std::setw(5) << NumCompute() 
00397        << "  " << std::setw(15) << ComputeTime()
00398        << "  " << std::setw(15) << 1.0e-6 * ComputeFlops(); 
00399     if (ComputeTime() != 0.0)
00400       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
00401     else
00402       os << "  " << std::setw(15) << 0.0 << endl;
00403     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00404        << "  " << std::setw(15) << ApplyInverseTime()
00405        << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
00406     if (ApplyInverseTime() != 0.0) 
00407       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
00408     else
00409       os << "  " << std::setw(15) << 0.0 << endl;
00410     os << "================================================================================" << endl;
00411     os << endl;
00412   }
00413 
00414   return(os);
00415 }

Generated on Wed May 12 21:46:03 2010 for IFPACK by  doxygen 1.4.7