43 #include "Ifpack_ConfigDefs.h"
44 #ifdef HAVE_IFPACK_SPARSKIT
45 #include "Ifpack_Preconditioner.h"
46 #include "Ifpack_SPARSKIT.h"
47 #include "Ifpack_Condest.h"
49 #include "Epetra_Comm.h"
50 #include "Epetra_Map.h"
51 #include "Epetra_RowMatrix.h"
52 #include "Epetra_Vector.h"
53 #include "Epetra_MultiVector.h"
54 #include "Epetra_Util.h"
55 #include "Teuchos_ParameterList.hpp"
57 #define F77_ILUT F77_FUNC(ilut, ILUT)
58 #define F77_ILUD F77_FUNC(ilud, ILUD)
59 #define F77_ILUTP F77_FUNC(ilutp, ILUTP)
60 #define F77_ILUDP F77_FUNC(iludp, ILUDP)
61 #define F77_ILUK F77_FUNC(iluk, ILUK)
62 #define F77_ILU0 F77_FUNC(ilu0, ILU0)
63 #define F77_LUSOL F77_FUNC(lusol, LUSOL)
66 void F77_ILUT(
int *,
double*,
int*,
int*,
int*,
double*,
67 double*,
int*,
int*,
int*,
double*,
int*,
int*);
68 void F77_ILUD(
int *,
double*,
int*,
int*,
double*,
double*,
69 double*,
int*,
int*,
int*,
double*,
int*,
int*);
70 void F77_ILUTP(
int *,
double*,
int*,
int*,
int*,
double*,
double*,
int*,
71 double*,
int*,
int*,
int*,
double*,
int*,
int*,
int*);
72 void F77_ILUDP(
int *,
double*,
int*,
int*,
double*,
double*,
double*,
int*,
73 double*,
int*,
int*,
int*,
double*,
int*,
int*,
int*);
74 void F77_ILUK(
int *,
double*,
int*,
int*,
int*,
75 double*,
int*,
int*,
int*,
int*,
double*,
int*,
int*);
76 void F77_ILU0(
int*,
double*,
int*,
int*,
double*,
int*,
int*,
int*,
int*);
77 void F77_LUSOL(
int *,
double*,
double*,
double*,
int*,
int*);
94 IsInitialized_(false),
101 ApplyInverseTime_(0),
103 ApplyInverseFlops_(0.0)
105 Teuchos::ParameterList List;
110 Ifpack_SPARSKIT::~Ifpack_SPARSKIT()
115 int Ifpack_SPARSKIT::SetParameters(Teuchos::ParameterList& List)
117 lfil_ = List.get(
"fact: sparskit: lfil", lfil_);
118 tol_ = List.get(
"fact: sparskit: tol", tol_);
119 droptol_ = List.get(
"fact: sparskit: droptol", droptol_);
120 permtol_ = List.get(
"fact: sparskit: permtol", permtol_);
121 alph_ = List.get(
"fact: sparskit: alph", alph_);
122 mbloc_ = List.get(
"fact: sparskit: mbloc", mbloc_);
123 Type_ = List.get(
"fact: sparskit: type", Type_);
126 Label_ =
"IFPACK SPARSKIT (Type=" + Type_ +
", fill=" +
133 int Ifpack_SPARSKIT::Initialize()
135 IsInitialized_ =
true;
141 int Ifpack_SPARSKIT::Compute()
143 if (!IsInitialized())
144 IFPACK_CHK_ERR(Initialize());
153 int n = Matrix().NumMyRows();
154 int nnz = Matrix().NumMyNonzeros();
156 std::vector<double> a(nnz);
157 std::vector<int> ja(nnz);
158 std::vector<int> ia(n + 1);
160 const int MaxNumEntries = Matrix().MaxNumEntries();
162 std::vector<double> Values(MaxNumEntries);
163 std::vector<int> Indices(MaxNumEntries);
168 for (
int i = 0 ; i < n ; ++i)
171 int NumMyEntries = 0;
172 Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumEntries, &Values[0],
178 for (
int j = 0 ; j < NumEntries ; ++j)
182 a[count] = Values[j];
183 ja[count] = Indices[j] + 1;
188 ia[i + 1] = ia[i] + NumMyEntries;
191 if (mbloc_ == -1) mbloc_ = n;
195 if (Type_ ==
"ILUT" || Type_ ==
"ILUTP" || Type_ ==
"ILUD" ||
197 iwk = nnz + 2 * lfil_ * n;
198 else if (Type_ ==
"ILUK")
199 iwk = (2 * lfil_ + 1) * nnz + 1;
200 else if (Type_ ==
"ILU0")
209 std::vector<int> jw(n + 1);
210 std::vector<double> w(n + 1);
215 F77_ILUT(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_,
216 &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
218 else if (Type_ ==
"ILUD")
221 F77_ILUD(&n, &a[0], &ja[0], &ia[0], &alph_, &tol_,
222 &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0], &ierr);
224 else if (Type_ ==
"ILUTP")
227 iperm_.resize(2 * n);
228 F77_ILUTP(&n, &a[0], &ja[0], &ia[0], &lfil_, &droptol_, &permtol_,
229 &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &iwk, &w[0], &jw[0],
231 for (
int i = 0 ; i < n ; ++i)
234 else if (Type_ ==
"ILUDP")
237 iperm_.resize(2 * n);
238 F77_ILUDP(&n, &a[0], &ja[0], &ia[0], &alph_, &droptol_, &permtol_,
239 &mbloc_, &alu_[0], &jlu_[0], &ju_[0], &n, &w[0], &jw[0],
241 for (
int i = 0 ; i < n ; ++i)
244 else if (Type_ ==
"ILUK")
246 std::vector<int> levs(iwk);
248 F77_ILUK(&n, &a[0], &ja[0], &ia[0], &lfil_,
249 &alu_[0], &jlu_[0], &ju_[0], &levs[0], &iwk, &w[0], &jw[0], &ierr);
251 else if (Type_ ==
"ILU0")
255 F77_ILU0(&n, &a[0], &ja[0], &ia[0],
256 &alu_[0], &jlu_[0], &ju_[0], &jw[0], &ierr);
258 IFPACK_CHK_ERR(ierr);
272 if (X.NumVectors() != Y.NumVectors())
275 int n = Matrix().NumMyRows();
277 for (
int i = 0 ; i < X.NumVectors() ; ++i)
278 F77_LUSOL(&n, (
double*)X(i)->Values(), Y(i)->Values(), (
double*)&alu_[0],
279 (
int*)&jlu_[0], (
int*)&ju_[0]);
282 if (Type_ ==
"ILUTP" || Type_ ==
"ILUDP")
284 std::vector<double> tmp(n);
285 for (
int j = 0 ; j < n ; ++j)
286 tmp[iperm_[j]] = Y[0][j];
287 for (
int j = 0 ; j < n ; ++j)
297 double Ifpack_SPARSKIT::Condest(
const Ifpack_CondestType CT,
298 const int MaxIters,
const double Tol,
304 if (Condest_ == -1.0)
305 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix);
312 Ifpack_SPARSKIT::Print(std::ostream& os)
const
315 if (!Comm().MyPID()) {
317 os <<
"================================================================================" << endl;
318 os <<
"Ifpack_SPARSKIT: " << Label() << endl << endl;
319 os <<
"Condition number estimate = " << Condest() << endl;
320 os <<
"Global number of rows = " << A_.NumGlobalRows() << endl;
322 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
323 os <<
"----- ------- -------------- ------------ --------" << endl;
324 os <<
"Initialize() " << std::setw(5) << NumInitialize()
325 <<
" " << std::setw(15) << InitializeTime()
326 <<
" 0.0 0.0" << endl;
327 os <<
"Compute() " << std::setw(5) << NumCompute()
328 <<
" " << std::setw(15) << ComputeTime()
329 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops();
330 if (ComputeTime() != 0.0)
331 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
333 os <<
" " << std::setw(15) << 0.0 << endl;
334 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse()
335 <<
" " << std::setw(15) << ApplyInverseTime()
336 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
337 if (ApplyInverseTime() != 0.0)
338 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
340 os <<
" " << std::setw(15) << 0.0 << endl;
341 os <<
"================================================================================" << endl;
349 #endif // HAVE_IFPACK_SPARSKIT
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.