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     for (int j = 0 ; j < NumEntries ; ++j)
00163     {
00164       if (Indices[j] < n) // skip non-local columns
00165       {
00166         a[count]  = Values[j];
00167         ja[count] = Indices[j] + 1; // SPARSKIT is FORTRAN
00168         ++count;
00169         ++NumMyEntries;
00170       }
00171     }
00172     ia[i + 1] = ia[i] + NumMyEntries;
00173   }
00174 
00175   if (mbloc_ == -1) mbloc_ = n;
00176 
00177   int iwk;
00178 
00179   if (Type_ == "ILUT" || Type_ == "ILUTP" || Type_ == "ILUD" ||
00180       Type_ == "ILUDP")
00181     iwk = nnz + 2 * lfil_ * n;
00182   else if (Type_ == "ILUK")
00183     iwk = (2 * lfil_ + 1) * nnz + 1;
00184   else if (Type_ == "ILU0")
00185     iwk = nnz;
00186 
00187   int ierr = 0;
00188 
00189   alu_.resize(iwk);
00190   jlu_.resize(iwk);
00191   ju_.resize(n + 1);
00192 
00193   vector<int>    jw(n + 1);
00194   vector<double> w(n + 1);
00195 
00196   if (Type_ == "ILUT")
00197   {
00198     jw.resize(2 * n);
00199     F77_ILUT(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_,
00200              &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
00201   }
00202   else if (Type_ == "ILUD")
00203   {
00204     jw.resize(2 * n);
00205     F77_ILUD(&n, &a[0], &ja[0], &ia[0], &alph_, &tol_,
00206              &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
00207   }
00208   else if (Type_ == "ILUTP")
00209   {
00210     jw.resize(2 * n);
00211     iperm_.resize(2 * n);
00212     F77_ILUTP(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_, &permtol_, 
00213               &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0],
00214               &iperm_[0], &ierr);
00215     for (int i = 0 ; i < n ; ++i)
00216       iperm_[i]--;
00217   }
00218   else if (Type_ == "ILUDP")
00219   {
00220     jw.resize(2 * n);
00221     iperm_.resize(2 * n);
00222     F77_ILUDP(&n, &a[0], &ja[0], &ia[0], &alph_, &droptol_, &permtol_, 
00223               &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &n, &w[0], &jw[0],
00224               &iperm_[0], &ierr);
00225     for (int i = 0 ; i < n ; ++i)
00226       iperm_[i]--;
00227   }
00228   else if (Type_ == "ILUK")
00229   {
00230     vector<int> levs(iwk);
00231     jw.resize(3 * n);
00232     F77_ILUK(&n, &a[0], &ja[0], &ia[0], &lfil_, 
00233              &alu_[0], &jlu_[0], &ju_[0], &levs[0], &iwk, &w[0], &jw[0], &ierr);
00234   }
00235   else if (Type_ == "ILU0")
00236   {
00237     // here w is only of size n
00238     jw.resize(2 * n);
00239     F77_ILU0(&n, &a[0], &ja[0], &ia[0], 
00240              &alu_[0], &jlu_[0], &ju_[0], &jw[0], &ierr);
00241   }
00242 
00243   IFPACK_CHK_ERR(ierr);
00244 
00245   IsComputed_ = true;
00246   return(0);
00247 }
00248 
00249 //=============================================================================
00250 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
00251 int Ifpack_SPARSKIT::ApplyInverse(const Epetra_MultiVector& X, 
00252                                   Epetra_MultiVector& Y) const
00253 {
00254   if (!IsComputed())
00255     IFPACK_CHK_ERR(-3); // compute preconditioner first
00256 
00257   if (X.NumVectors() != Y.NumVectors()) 
00258     IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
00259 
00260   int n = Matrix().NumMyRows();
00261 
00262   for (int i = 0 ; i < X.NumVectors() ; ++i)
00263       F77_LUSOL(&n, (double*)Y(i)->Values(), X(i)->Values(), (double*)&alu_[0], 
00264                 (int*)&jlu_[0], (int*)&ju_[0]);
00265 
00266 #if 0 
00267   // still need to fix support for permutation
00268   if (Type_ == "ILUTP" || Type_ == "ILUDP")
00269   {
00270     vector<double> tmp(n);
00271     for (int j = 0 ; j < n ; ++j)
00272       tmp[j] = Y[0][iperm_[j]];
00273     for (int j = 0 ; j < n ; ++j)
00274       Y[0][j] = tmp[j];
00275   }
00276 #endif
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

Generated on Tue Oct 20 12:48:54 2009 for IFPACK by doxygen 1.4.7