IFPACK Development
Ifpack_SPARSKIT.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 #ifdef HAVE_IFPACK_SPARSKIT
00032 #include "Ifpack_Preconditioner.h"
00033 #include "Ifpack_SPARSKIT.h"
00034 #include "Ifpack_Condest.h"
00035 #include "Ifpack_Utils.h"
00036 #include "Epetra_Comm.h"
00037 #include "Epetra_Map.h"
00038 #include "Epetra_RowMatrix.h"
00039 #include "Epetra_Vector.h"
00040 #include "Epetra_MultiVector.h"
00041 #include "Epetra_Util.h"
00042 #include "Teuchos_ParameterList.hpp"
00043 
00044 #define F77_ILUT  F77_FUNC(ilut, ILUT)
00045 #define F77_ILUD  F77_FUNC(ilud, ILUD)
00046 #define F77_ILUTP F77_FUNC(ilutp, ILUTP)
00047 #define F77_ILUDP F77_FUNC(iludp, ILUDP)
00048 #define F77_ILUK  F77_FUNC(iluk, ILUK)
00049 #define F77_ILU0  F77_FUNC(ilu0, ILU0)
00050 #define F77_LUSOL F77_FUNC(lusol, LUSOL)
00051 
00052 extern "C" {
00053   void F77_ILUT(int *, double*, int*, int*, int*, double*,
00054                 double*, int*, int*, int*, double*, int*, int*);
00055   void F77_ILUD(int *, double*, int*, int*, double*, double*,
00056                 double*, int*, int*, int*, double*, int*, int*);
00057   void F77_ILUTP(int *, double*, int*, int*, int*, double*, double*, int*,
00058                  double*, int*, int*, int*, double*, int*, int*, int*);
00059   void F77_ILUDP(int *, double*, int*, int*, double*, double*, double*, int*,
00060                  double*, int*, int*, int*, double*, int*, int*, int*);
00061   void F77_ILUK(int *, double*, int*, int*, int*, 
00062                 double*, int*, int*, int*, int*, double*, int*, int*);
00063   void F77_ILU0(int*, double*, int*, int*, double*, int*, int*, int*, int*);
00064   void F77_LUSOL(int *, double*, double*, double*, int*, int*);
00065 }
00066 
00067 
00068 //==============================================================================
00069 Ifpack_SPARSKIT::Ifpack_SPARSKIT(Epetra_RowMatrix* A) :
00070   A_(*A),
00071   Comm_(A->Comm()),
00072   UseTranspose_(false),
00073   lfil_(0),
00074   droptol_(0.0),
00075   tol_(0.0),
00076   permtol_(0.1),
00077   alph_(0.0),
00078   mbloc_(-1),
00079   Type_("ILUT"),
00080   Condest_(-1.0),
00081   IsInitialized_(false),
00082   IsComputed_(false),
00083   NumInitialize_(0),
00084   NumCompute_(0),
00085   NumApplyInverse_(0),
00086   InitializeTime_(0.0),
00087   ComputeTime_(0),
00088   ApplyInverseTime_(0),
00089   ComputeFlops_(0.0),
00090   ApplyInverseFlops_(0.0)
00091 {
00092   Teuchos::ParameterList List;
00093   SetParameters(List);
00094 }
00095 
00096 //==============================================================================
00097 Ifpack_SPARSKIT::~Ifpack_SPARSKIT()
00098 {
00099 }
00100 
00101 //==========================================================================
00102 int Ifpack_SPARSKIT::SetParameters(Teuchos::ParameterList& List)
00103 {
00104   lfil_    = List.get("fact: sparskit: lfil", lfil_);
00105   tol_     = List.get("fact: sparskit: tol", tol_);
00106   droptol_ = List.get("fact: sparskit: droptol", droptol_);
00107   permtol_ = List.get("fact: sparskit: permtol", permtol_);
00108   alph_    = List.get("fact: sparskit: alph", alph_);
00109   mbloc_   = List.get("fact: sparskit: mbloc", mbloc_);
00110   Type_    = List.get("fact: sparskit: type", Type_);
00111 
00112   // set label
00113   Label_ = "IFPACK SPARSKIT (Type=" + Type_ + ", fill=" +
00114     Ifpack_toString(lfil_) + ")";
00115 
00116   return(0);
00117 }
00118 
00119 //==========================================================================
00120 int Ifpack_SPARSKIT::Initialize()
00121 {
00122   IsInitialized_ = true;
00123   IsComputed_    = false;
00124   return(0);
00125 }
00126 
00127 //==========================================================================
00128 int Ifpack_SPARSKIT::Compute() 
00129 {
00130   if (!IsInitialized()) 
00131     IFPACK_CHK_ERR(Initialize());
00132 
00133   IsComputed_ = false;
00134 
00135   // convert the matrix into SPARSKIT format. The matrix is then
00136   // free'd after method Compute() returns.
00137 
00138   // convert the matrix into CSR format. Note that nnz is an over-estimate,
00139   // since it contains the nonzeros corresponding to external nodes as well.
00140   int n   = Matrix().NumMyRows();
00141   int nnz = Matrix().NumMyNonzeros();
00142 
00143   vector<double> a(nnz);
00144   vector<int>    ja(nnz);
00145   vector<int>    ia(n + 1);
00146 
00147   const int MaxNumEntries = Matrix().MaxNumEntries();
00148 
00149   vector<double> Values(MaxNumEntries);
00150   vector<int>    Indices(MaxNumEntries);
00151 
00152   int count = 0;
00153 
00154   ia[0] = 1;
00155   for (int i = 0 ; i < n ; ++i)
00156   {
00157     int NumEntries;
00158     int NumMyEntries = 0;
00159     Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0], 
00160                               &Indices[0]);
00161 
00162     // NOTE: There might be some issues here with the ILU(0) if the column indices aren't sorted.
00163     // The other factorizations are probably OK.
00164 
00165     for (int j = 0 ; j < NumEntries ; ++j)
00166     {
00167       if (Indices[j] < n) // skip non-local columns
00168       {
00169         a[count]  = Values[j];
00170         ja[count] = Indices[j] + 1; // SPARSKIT is FORTRAN
00171         ++count;
00172         ++NumMyEntries;
00173       }
00174     }
00175     ia[i + 1] = ia[i] + NumMyEntries;
00176   }
00177 
00178   if (mbloc_ == -1) mbloc_ = n;
00179 
00180   int iwk;
00181 
00182   if (Type_ == "ILUT" || Type_ == "ILUTP" || Type_ == "ILUD" ||
00183       Type_ == "ILUDP")
00184     iwk = nnz + 2 * lfil_ * n;
00185   else if (Type_ == "ILUK")
00186     iwk = (2 * lfil_ + 1) * nnz + 1;
00187   else if (Type_ == "ILU0")
00188     iwk = nnz+2;
00189 
00190   int ierr = 0;
00191 
00192   alu_.resize(iwk);
00193   jlu_.resize(iwk);
00194   ju_.resize(n + 1);
00195 
00196   vector<int>    jw(n + 1);
00197   vector<double> w(n + 1);
00198 
00199   if (Type_ == "ILUT")
00200   {
00201     jw.resize(2 * n);
00202     F77_ILUT(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_,
00203              &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
00204   }
00205   else if (Type_ == "ILUD")
00206   {
00207     jw.resize(2 * n);
00208     F77_ILUD(&n, &a[0], &ja[0], &ia[0], &alph_, &tol_,
00209              &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
00210   }
00211   else if (Type_ == "ILUTP")
00212   {
00213     jw.resize(2 * n);
00214     iperm_.resize(2 * n);
00215     F77_ILUTP(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_, &permtol_, 
00216               &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0],
00217               &iperm_[0], &ierr);
00218     for (int i = 0 ; i < n ; ++i)
00219       iperm_[i]--;
00220   }
00221   else if (Type_ == "ILUDP")
00222   {
00223     jw.resize(2 * n);
00224     iperm_.resize(2 * n);
00225     F77_ILUDP(&n, &a[0], &ja[0], &ia[0], &alph_, &droptol_, &permtol_, 
00226               &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &n, &w[0], &jw[0],
00227               &iperm_[0], &ierr);
00228     for (int i = 0 ; i < n ; ++i)
00229       iperm_[i]--;
00230   }
00231   else if (Type_ == "ILUK")
00232   {
00233     vector<int> levs(iwk);
00234     jw.resize(3 * n);
00235     F77_ILUK(&n, &a[0], &ja[0], &ia[0], &lfil_, 
00236              &alu_[0], &jlu_[0], &ju_[0], &levs[0], &iwk, &w[0], &jw[0], &ierr);
00237   }
00238   else if (Type_ == "ILU0")
00239   {
00240     // here w is only of size n
00241     jw.resize(2 * n);
00242     F77_ILU0(&n, &a[0], &ja[0], &ia[0], 
00243              &alu_[0], &jlu_[0], &ju_[0], &jw[0], &ierr);
00244   }
00245   IFPACK_CHK_ERR(ierr);
00246 
00247   IsComputed_ = true;
00248   return(0);
00249 }
00250 
00251 //=============================================================================
00252 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00253 int Ifpack_SPARSKIT::ApplyInverse(const Epetra_MultiVector& X, 
00254                                   Epetra_MultiVector& Y) const
00255 {
00256   if (!IsComputed())
00257     IFPACK_CHK_ERR(-3); // compute preconditioner first
00258 
00259   if (X.NumVectors() != Y.NumVectors()) 
00260     IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
00261 
00262   int n = Matrix().NumMyRows();
00263 
00264   for (int i = 0 ; i < X.NumVectors() ; ++i)
00265     F77_LUSOL(&n, (double*)X(i)->Values(), Y(i)->Values(), (double*)&alu_[0], 
00266                 (int*)&jlu_[0], (int*)&ju_[0]);
00267 
00268   // still need to fix support for permutation
00269   if (Type_ == "ILUTP" || Type_ == "ILUDP")
00270   {
00271     vector<double> tmp(n);
00272     for (int j = 0 ; j < n ; ++j)
00273       tmp[iperm_[j]] = Y[0][j];
00274     for (int j = 0 ; j < n ; ++j)
00275       Y[0][j] = tmp[j];
00276   }
00277 
00278   ++NumApplyInverse_;
00279   return(0);
00280 
00281 }
00282 
00283 //=============================================================================
00284 double Ifpack_SPARSKIT::Condest(const Ifpack_CondestType CT, 
00285                                 const int MaxIters, const double Tol,
00286                                 Epetra_RowMatrix* Matrix)
00287 {
00288   if (!IsComputed()) // cannot compute right now
00289     return(-1.0);
00290 
00291   if (Condest_ == -1.0)
00292     Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);
00293 
00294   return(Condest_);
00295 }
00296 
00297 //=============================================================================
00298 std::ostream&
00299 Ifpack_SPARSKIT::Print(std::ostream& os) const
00300 {
00301   if (!Comm().MyPID()) {
00302     os << endl;
00303     os << "================================================================================" << endl;
00304     os << "Ifpack_SPARSKIT: " << Label() << endl << endl;
00305     os << "Condition number estimate = " << Condest() << endl;
00306     os << "Global number of rows            = " << A_.NumGlobalRows() << endl;
00307     os << endl;
00308     os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
00309     os << "-----           -------   --------------       ------------     --------" << endl;
00310     os << "Initialize()    "   << std::setw(5) << NumInitialize() 
00311       << "  " << std::setw(15) << InitializeTime() 
00312       << "               0.0            0.0" << endl;
00313     os << "Compute()       "   << std::setw(5) << NumCompute() 
00314       << "  " << std::setw(15) << ComputeTime()
00315       << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
00316     if (ComputeTime() != 0.0)
00317       os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
00318     else
00319       os << "  " << std::setw(15) << 0.0 << endl;
00320     os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
00321       << "  " << std::setw(15) << ApplyInverseTime()
00322       << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
00323     if (ApplyInverseTime() != 0.0) 
00324       os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
00325     else
00326       os << "  " << std::setw(15) << 0.0 << endl;
00327     os << "================================================================================" << endl;
00328     os << endl;
00329   }
00330 
00331 
00332   return(os);
00333 } 
00334 
00335 #endif // HAVE_IFPACK_SPARSKIT
 All Classes Files Functions Variables Enumerations Friends