00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "Ifpack_ConfigDefs.h"
00031 #if defined(HAVE_IFPACK_TEUCHOS) && defined(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
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
00136
00137
00138
00139
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)
00165 {
00166 a[count] = Values[j];
00167 ja[count] = Indices[j] + 1;
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
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
00251 int Ifpack_SPARSKIT::ApplyInverse(const Epetra_MultiVector& X,
00252 Epetra_MultiVector& Y) const
00253 {
00254 if (!IsComputed())
00255 IFPACK_CHK_ERR(-3);
00256
00257 if (X.NumVectors() != Y.NumVectors())
00258 IFPACK_CHK_ERR(-2);
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
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())
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 #endif // HAVE_IFPACK_TEUCHOS