43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_CondestType.h"
45 #include "Ifpack_ILU.h"
46 #include "Epetra_ConfigDefs.h"
47 #include "Epetra_Comm.h"
48 #include "Epetra_Map.h"
49 #include "Epetra_RowMatrix.h"
50 #include "Epetra_Vector.h"
51 #include "Epetra_MultiVector.h"
52 #include "Epetra_CrsGraph.h"
53 #include "Epetra_CrsMatrix.h"
54 #include "Teuchos_ParameterList.hpp"
55 #include "Teuchos_RefCountPtr.hpp"
57 using Teuchos::RefCountPtr;
60 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
61 # include "Teuchos_TimeMonitor.hpp"
66 A_(rcp(Matrix_in,false)),
67 Comm_(Matrix_in->Comm()),
75 IsInitialized_(false),
82 ApplyInverseTime_(0.0),
84 ApplyInverseFlops_(0.0),
87 Teuchos::ParameterList List;
92 void Ifpack_ILU::Destroy()
102 RelaxValue_ = List.get(
"fact: relax value", RelaxValue_);
103 Athresh_ = List.get(
"fact: absolute threshold", Athresh_);
104 Rthresh_ = List.get(
"fact: relative threshold", Rthresh_);
105 LevelOfFill_ = List.get(
"fact: level-of-fill", LevelOfFill_);
108 sprintf(Label_,
"IFPACK ILU (fill=%d, relax=%f, athr=%f, rthr=%f)",
109 LevelOfFill(), RelaxValue(), AbsoluteThreshold(),
110 RelativeThreshold());
115 int Ifpack_ILU::ComputeSetup()
118 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
119 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_ILU::ComputeSetup");
125 if ((L_.get() == 0) || (U_.get() == 0) || (D_.get() == 0))
138 int NumIn, NumL, NumU;
140 int NumNonzeroDiags = 0;
142 std::vector<int> InI(MaxNumEntries);
143 std::vector<int> LI(MaxNumEntries);
144 std::vector<int> UI(MaxNumEntries);
145 std::vector<double> InV(MaxNumEntries);
146 std::vector<double> LV(MaxNumEntries);
147 std::vector<double> UV(MaxNumEntries);
149 bool ReplaceValues = (L_->StaticGraph() || L_->IndicesAreLocal());
158 IFPACK_CHK_ERR(D_->ExtractView(&DV));
162 for (i=0; i< NumMyRows(); i++) {
164 IFPACK_CHK_ERR(
Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]));
172 for (j=0; j< NumIn; j++) {
177 DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_;
180 else if (k < 0) {IFPACK_CHK_ERR(-4);}
187 else if (k<NumMyRows()) {
196 if (DiagFound) NumNonzeroDiags++;
197 else DV[i] = Athresh_;
201 (L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
204 IFPACK_CHK_ERR(L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));
210 (U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
213 IFPACK_CHK_ERR(U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));
219 if (!ReplaceValues) {
224 IFPACK_CHK_ERR(L_->FillComplete((L_->RowMatrixColMap()), *L_RangeMap_));
225 IFPACK_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));
229 int TotalNonzeroDiags = 0;
230 IFPACK_CHK_ERR(Graph().L_Graph().RowMap().
Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1));
231 NumMyDiagonals_ = NumNonzeroDiags;
232 if (NumNonzeroDiags != NumMyRows()) ierr = 1;
234 IFPACK_CHK_ERR(ierr);
242 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
243 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_ILU::Initialize");
247 IsInitialized_ =
false;
254 if (CrsMatrix == 0) {
259 int size = A_->MaxNumEntries();
261 if (CrsGraph_.get() == 0)
264 std::vector<double> Values(size);
266 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
267 if(A_->RowMatrixRowMap().GlobalIndicesInt()) {
268 std::vector<int> Indices(size);
271 for (
int i = 0 ; i < A_->NumMyRows() ; ++i) {
273 int GlobalRow = A_->RowMatrixRowMap().GID(i);
274 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries,
275 &Values[0], &Indices[0]));
277 for (
int j = 0 ; j < NumEntries ; ++j) {
278 Indices[j] = A_->RowMatrixColMap().GID(Indices[j]);
280 IFPACK_CHK_ERR(CrsGraph_->InsertGlobalIndices(GlobalRow,NumEntries,
286 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
287 if(A_->RowMatrixRowMap().GlobalIndicesLongLong()) {
288 std::vector<int> Indices_local(size);
289 std::vector<long long> Indices(size);
292 for (
int i = 0 ; i < A_->NumMyRows() ; ++i) {
294 long long GlobalRow = A_->RowMatrixRowMap().GID64(i);
295 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries,
296 &Values[0], &Indices_local[0]));
298 for (
int j = 0 ; j < NumEntries ; ++j) {
299 Indices[j] = A_->RowMatrixColMap().GID64(Indices_local[j]);
301 IFPACK_CHK_ERR(CrsGraph_->InsertGlobalIndices(GlobalRow,NumEntries,
307 throw "Ifpack_ILU::Initialize: GlobalIndices type unknown";
309 IFPACK_CHK_ERR(CrsGraph_->FillComplete(A_->RowMatrixRowMap(),
310 A_->RowMatrixRowMap()));
322 if (Graph_.get() == 0)
324 IFPACK_CHK_ERR(Graph_->ConstructFilledGraph());
326 IsInitialized_ =
true;
337 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
338 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_ILU::Compute");
348 IFPACK_CHK_ERR(ComputeSetup());
353 double MinDiagonalValue = Epetra_MinDouble;
354 double MaxDiagonalValue = 1.0/MinDiagonalValue;
360 int NumIn, NumL, NumU;
363 int MaxNumEntries = L_->MaxNumEntries() + U_->MaxNumEntries() + 1;
365 std::vector<int> InI(MaxNumEntries+1);
366 std::vector<double> InV(MaxNumEntries+1);
367 std::vector<int> colflag(NumMyCols());
370 ierr = D_->ExtractView(&DV);
372 #ifdef IFPACK_FLOPCOUNTERS
373 int current_madds = 0;
384 for (j = 0; j < NumMyCols(); ++j) colflag[j] = - 1;
386 for (i = 0; i < NumMyRows(); ++i) {
390 NumIn = MaxNumEntries;
391 IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));
398 IFPACK_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
404 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
406 double diagmod = 0.0;
408 for (
int jj=0; jj<NumL; jj++) {
410 double multiplier = InV[jj];
414 IFPACK_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI));
416 if (RelaxValue_==0.0) {
417 for (k=0; k<NumUU; k++) {
418 int kk = colflag[UUI[k]];
420 InV[kk] -= multiplier*UUV[k];
421 #ifdef IFPACK_FLOPCOUNTERS
428 for (k=0; k<NumUU; k++) {
429 int kk = colflag[UUI[k]];
430 if (kk>-1) InV[kk] -= multiplier*UUV[k];
431 else diagmod -= multiplier*UUV[k];
432 #ifdef IFPACK_FLOPCOUNTERS
439 IFPACK_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));
444 if (RelaxValue_!=0.0) {
445 DV[i] += RelaxValue_*diagmod;
449 if (fabs(DV[i]) > MaxDiagonalValue) {
450 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
451 else DV[i] = MinDiagonalValue;
456 for (j=0; j<NumU; j++) UV[j] *= DV[i];
459 IFPACK_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI));
463 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
468 if (!L_->LowerTriangular())
470 if (!U_->UpperTriangular())
473 #ifdef IFPACK_FLOPCOUNTERS
476 double current_flops = 2 * current_madds;
477 double total_flops = 0;
479 IFPACK_CHK_ERR(Graph().L_Graph().RowMap().
Comm().SumAll(¤t_flops, &total_flops, 1));
482 total_flops += (double) L_->NumGlobalNonzeros();
483 total_flops += (double) D_->GlobalLength();
484 if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength();
487 ComputeFlops_ += total_flops;
504 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
505 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_ILU::ApplyInverse - Solve");
511 bool UnitDiagonal =
true;
515 IFPACK_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, X, Y));
517 IFPACK_CHK_ERR(Y.Multiply(1.0, *D_, Y, 0.0));
519 IFPACK_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, Y, Y));
523 IFPACK_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, X, Y));
525 IFPACK_CHK_ERR(Y.Multiply(1.0, *D_, Y, 0.0));
526 IFPACK_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, Y, Y));
539 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
540 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_ILU::Multiply");
547 IFPACK_CHK_ERR(U_->Multiply(Trans, X, Y));
549 IFPACK_CHK_ERR(Y.Update(1.0, X, 1.0));
551 IFPACK_CHK_ERR(Y.ReciprocalMultiply(1.0, *D_, Y, 0.0));
553 IFPACK_CHK_ERR(L_->Multiply(Trans, Y1temp, Y));
555 IFPACK_CHK_ERR(Y.Update(1.0, Y1temp, 1.0));
559 IFPACK_CHK_ERR(L_->Multiply(Trans, X, Y));
561 IFPACK_CHK_ERR(Y.Update(1.0, X, 1.0));
563 IFPACK_CHK_ERR(Y.ReciprocalMultiply(1.0, *D_, Y, 0.0));
565 IFPACK_CHK_ERR(U_->Multiply(Trans, Y1temp, Y));
567 IFPACK_CHK_ERR(Y.Update(1.0, Y1temp, 1.0));
579 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
580 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_ILU::ApplyInverse");
586 if (X.NumVectors() != Y.NumVectors())
593 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
594 if (X.Pointers()[0] == Y.Pointers()[0])
597 Xcopy = Teuchos::rcp( &X,
false );
602 #ifdef IFPACK_FLOPCOUNTERS
603 ApplyInverseFlops_ += X.NumVectors() * 4 *
604 (L_->NumGlobalNonzeros() + U_->NumGlobalNonzeros());
616 const int MaxIters,
const double Tol,
620 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
621 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_ILU::Condest");
627 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
638 if (!
Comm().MyPID()) {
640 os <<
"================================================================================" << endl;
641 os <<
"Ifpack_ILU: " <<
Label() << endl << endl;
642 os <<
"Level-of-fill = " << LevelOfFill() << endl;
643 os <<
"Absolute threshold = " << AbsoluteThreshold() << endl;
644 os <<
"Relative threshold = " << RelativeThreshold() << endl;
645 os <<
"Relax value = " << RelaxValue() << endl;
646 os <<
"Condition number estimate = " <<
Condest() << endl;
647 os <<
"Global number of rows = " << A_->NumGlobalRows64() << endl;
649 os <<
"Number of rows of L, D, U = " << L_->NumGlobalRows64() << endl;
650 os <<
"Number of nonzeros of L + U = " << NumGlobalNonzeros64() << endl;
651 os <<
"nonzeros / rows = "
652 << 1.0 * NumGlobalNonzeros64() / U_->NumGlobalRows64() << endl;
655 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
656 os <<
"----- ------- -------------- ------------ --------" << endl;
659 <<
" 0.0 0.0" << endl;
660 os <<
"Compute() " << std::setw(5) <<
NumCompute()
666 os <<
" " << std::setw(15) << 0.0 << endl;
673 os <<
" " << std::setw(15) << 0.0 << endl;
674 os <<
"================================================================================" << endl;
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y...
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
double ElapsedTime(void) const
bool UseTranspose() const
Returns the current UseTranspose setting.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
virtual const Epetra_Map & OperatorDomainMap() const =0
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual int MaxNumEntries() const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
int Compute()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
int SetParameters(Teuchos::ParameterList ¶meterlist)
Set parameters using a Teuchos::ParameterList object.
int Initialize()
Initialize the preconditioner, does not touch matrix values.
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
double Condest() const
Returns the computed estimated condition number, or -1.0 if not computed.
virtual std::ostream & Print(std::ostream &os) const
Prints on stream basic information about this object.
Ifpack_ILU(Epetra_RowMatrix *A)
Constructor.
const Epetra_CrsGraph & Graph() const
const char * Label() const
Returns a character string describing the operator.
void ResetStartTime(void)
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual int NumInitialize() const
Returns the number of calls to Initialize().