Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_ILU.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
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"
56
57using Teuchos::RefCountPtr;
58using Teuchos::rcp;
59
60#ifdef IFPACK_TEUCHOS_TIME_MONITOR
61# include "Teuchos_TimeMonitor.hpp"
62#endif
63
64//==============================================================================
66 A_(rcp(Matrix_in,false)),
67 Comm_(Matrix_in->Comm()),
68 UseTranspose_(false),
69 NumMyDiagonals_(0),
70 RelaxValue_(0.0),
71 Athresh_(0.0),
72 Rthresh_(1.0),
73 Condest_(-1.0),
74 LevelOfFill_(0),
75 IsInitialized_(false),
76 IsComputed_(false),
77 NumInitialize_(0),
78 NumCompute_(0),
79 NumApplyInverse_(0),
80 InitializeTime_(0.0),
81 ComputeTime_(0.0),
82 ApplyInverseTime_(0.0),
83 ComputeFlops_(0.0),
84 ApplyInverseFlops_(0.0),
85 Time_(Comm())
86{
87 Teuchos::ParameterList List;
88 SetParameters(List);
89}
90
91//==============================================================================
93{
94 // reset pointers to already allocated stuff
95 U_DomainMap_ = 0;
96 L_RangeMap_ = 0;
97}
98
99//==========================================================================
100int Ifpack_ILU::SetParameters(Teuchos::ParameterList& List)
101{
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_);
106
107 // set label
108 sprintf(Label_, "IFPACK ILU (fill=%d, relax=%f, athr=%f, rthr=%f)",
111 return(0);
112}
113
114//==========================================================================
116{
117
118#ifdef IFPACK_TEUCHOS_TIME_MONITOR
119 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_ILU::ComputeSetup");
120#endif
121
122 L_ = rcp(new Epetra_CrsMatrix(Copy, Graph().L_Graph()));
123 U_ = rcp(new Epetra_CrsMatrix(Copy, Graph().U_Graph()));
124 D_ = rcp(new Epetra_Vector(Graph().L_Graph().RowMap()));
125 if ((L_.get() == 0) || (U_.get() == 0) || (D_.get() == 0))
126 IFPACK_CHK_ERR(-5);
127
128 // Get Maximun Row length
129 int MaxNumEntries = Matrix().MaxNumEntries();
130
131 // Set L range map and U domain map
132 U_DomainMap_ = &(Matrix().OperatorDomainMap());
133 L_RangeMap_ = &(Matrix().OperatorRangeMap());
134
135 // this is the old InitAllValues()
136 int ierr = 0;
137 int i, j;
138 int NumIn, NumL, NumU;
139 bool DiagFound;
140 int NumNonzeroDiags = 0;
141
142 std::vector<int> InI(MaxNumEntries); // Allocate temp space
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);
148
149 bool ReplaceValues = (L_->StaticGraph() || L_->IndicesAreLocal()); // Check if values should be inserted or replaced
150
151 if (ReplaceValues) {
152 L_->PutScalar(0.0); // Zero out L and U matrices
153 U_->PutScalar(0.0);
154 }
155
156 D_->PutScalar(0.0); // Set diagonal values to zero
157 double *DV;
158 IFPACK_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
159
160 // First we copy the user's matrix into L and U, regardless of fill level
161
162 for (i=0; i< NumMyRows(); i++) {
163
164 IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0])); // Get Values and Indices
165
166 // Split into L and U (we don't assume that indices are ordered).
167
168 NumL = 0;
169 NumU = 0;
170 DiagFound = false;
171
172 for (j=0; j< NumIn; j++) {
173 int k = InI[j];
174
175 if (k==i) {
176 DiagFound = true;
177 DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
178 }
179
180 else if (k < 0) {IFPACK_CHK_ERR(-4);} // Out of range
181
182 else if (k < i) {
183 LI[NumL] = k;
184 LV[NumL] = InV[j];
185 NumL++;
186 }
187 else if (k<NumMyRows()) {
188 UI[NumU] = k;
189 UV[NumU] = InV[j];
190 NumU++;
191 }
192 }
193
194 // Check in things for this row of L and U
195
196 if (DiagFound) NumNonzeroDiags++;
197 else DV[i] = Athresh_;
198
199 if (NumL) {
200 if (ReplaceValues) {
201 (L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
202 }
203 else {
204 IFPACK_CHK_ERR(L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));
205 }
206 }
207
208 if (NumU) {
209 if (ReplaceValues) {
210 (U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
211 }
212 else {
213 IFPACK_CHK_ERR(U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));
214 }
215 }
216
217 }
218
219 if (!ReplaceValues) {
220 // The domain of L and the range of U are exactly their own row maps (there is no communication).
221 // The domain of U and the range of L must be the same as those of the original matrix,
222 // However if the original matrix is a VbrMatrix, these two latter maps are translation from
223 // a block map to a point map.
224 IFPACK_CHK_ERR(L_->FillComplete((L_->RowMatrixColMap()), *L_RangeMap_));
225 IFPACK_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));
226 }
227
228 // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
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; // Diagonals are not right, warn user
233
234 IFPACK_CHK_ERR(ierr);
235 return(ierr);
236}
237
238//==========================================================================
240{
241
242#ifdef IFPACK_TEUCHOS_TIME_MONITOR
243 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_ILU::Initialize");
244#endif
245
247 IsInitialized_ = false;
248
249 // reset this object
250 Destroy();
251
252 Epetra_CrsMatrix* CrsMatrix;
253 CrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(&*A_);
254 if (CrsMatrix == 0) {
255 // this means that we have to create
256 // the graph from a given Epetra_RowMatrix. Note
257 // that at this point we are ignoring any possible
258 // graph coming from VBR matrices.
259 int size = A_->MaxNumEntries();
260 CrsGraph_ = rcp(new Epetra_CrsGraph(Copy,A_->RowMatrixRowMap(), size));
261 if (CrsGraph_.get() == 0)
262 IFPACK_CHK_ERR(-5); // memory allocation error
263
264 std::vector<double> Values(size);
265
266#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
267 if(A_->RowMatrixRowMap().GlobalIndicesInt()) {
268 std::vector<int> Indices(size);
269 // extract each row at-a-time, and insert it into
270 // the graph, ignore all off-process entries
271 for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
272 int NumEntries;
273 int GlobalRow = A_->RowMatrixRowMap().GID(i);
274 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries,
275 &Values[0], &Indices[0]));
276 // convert to global indices
277 for (int j = 0 ; j < NumEntries ; ++j) {
278 Indices[j] = A_->RowMatrixColMap().GID(Indices[j]);
279 }
280 IFPACK_CHK_ERR(CrsGraph_->InsertGlobalIndices(GlobalRow,NumEntries,
281 &Indices[0]));
282 }
283 }
284 else
285#endif
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);
290 // extract each row at-a-time, and insert it into
291 // the graph, ignore all off-process entries
292 for (int i = 0 ; i < A_->NumMyRows() ; ++i) {
293 int NumEntries;
294 long long GlobalRow = A_->RowMatrixRowMap().GID64(i);
295 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries,
296 &Values[0], &Indices_local[0]));
297 // convert to global indices
298 for (int j = 0 ; j < NumEntries ; ++j) {
299 Indices[j] = A_->RowMatrixColMap().GID64(Indices_local[j]);
300 }
301 IFPACK_CHK_ERR(CrsGraph_->InsertGlobalIndices(GlobalRow,NumEntries,
302 &Indices[0]));
303 }
304 }
305 else
306#endif
307 throw "Ifpack_ILU::Initialize: GlobalIndices type unknown";
308
309 IFPACK_CHK_ERR(CrsGraph_->FillComplete(A_->RowMatrixRowMap(),
310 A_->RowMatrixRowMap()));
311
312 // always overlap zero, wider overlap will be handled
313 // by the AdditiveSchwarz preconditioner.
315
316 }
317 else {
318 // see comment above for the overlap.
319 Graph_ = rcp(new Ifpack_IlukGraph(CrsMatrix->Graph(), LevelOfFill_, 0));
320 }
321
322 if (Graph_.get() == 0)
323 IFPACK_CHK_ERR(-5); // memory allocation error
324 IFPACK_CHK_ERR(Graph_->ConstructFilledGraph());
325
326 IsInitialized_ = true;
329
330 return(0);
331}
332
333//==========================================================================
335{
336
337#ifdef IFPACK_TEUCHOS_TIME_MONITOR
338 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_ILU::Compute");
339#endif
340
341 if (!IsInitialized())
343
345 IsComputed_ = false;
346
347 // convert Matrix() into L and U factors.
349
350 // MinMachNum should be officially defined, for now pick something a little
351 // bigger than IEEE underflow value
352
353 double MinDiagonalValue = Epetra_MinDouble;
354 double MaxDiagonalValue = 1.0/MinDiagonalValue;
355
356 int ierr = 0;
357 int i, j, k;
358 int *LI, *UI;
359 double *LV, *UV;
360 int NumIn, NumL, NumU;
361
362 // Get Maximun Row length
363 int MaxNumEntries = L_->MaxNumEntries() + U_->MaxNumEntries() + 1;
364
365 std::vector<int> InI(MaxNumEntries+1); // Allocate temp space, pad by one to
366 std::vector<double> InV(MaxNumEntries+1); // to avoid debugger complaints for pathological cases
367 std::vector<int> colflag(NumMyCols());
368
369 double *DV;
370 ierr = D_->ExtractView(&DV); // Get view of diagonal
371
372#ifdef IFPACK_FLOPCOUNTERS
373 int current_madds = 0; // We will count multiply-add as they happen
374#endif
375
376 // =========================== //
377 // Now start the factorization //
378 // =========================== //
379
380 // Need some integer workspace and pointers
381 int NumUU;
382 int * UUI;
383 double * UUV;
384 for (j = 0; j < NumMyCols(); ++j) colflag[j] = - 1;
385
386 for (i = 0; i < NumMyRows(); ++i) {
387
388 // Fill InV, InI with current row of L, D and U combined
389
390 NumIn = MaxNumEntries;
391 IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));
392 LV = &InV[0];
393 LI = &InI[0];
394
395 InV[NumL] = DV[i]; // Put in diagonal
396 InI[NumL] = i;
397
398 IFPACK_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
399 NumIn = NumL+NumU+1;
400 UV = &InV[NumL+1];
401 UI = &InI[NumL+1];
402
403 // Set column flags
404 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
405
406 double diagmod = 0.0; // Off-diagonal accumulator
407
408 for (int jj=0; jj<NumL; jj++) {
409 j = InI[jj];
410 double multiplier = InV[jj]; // current_mults++;
411
412 InV[jj] *= DV[j];
413
414 IFPACK_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
415
416 if (RelaxValue_==0.0) {
417 for (k=0; k<NumUU; k++) {
418 int kk = colflag[UUI[k]];
419 if (kk>-1) {
420 InV[kk] -= multiplier*UUV[k];
421#ifdef IFPACK_FLOPCOUNTERS
422 current_madds++;
423#endif
424 }
425 }
426 }
427 else {
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
433 current_madds++;
434#endif
435 }
436 }
437 }
438 if (NumL) {
439 IFPACK_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI)); // Replace current row of L
440 }
441
442 DV[i] = InV[NumL]; // Extract Diagonal value
443
444 if (RelaxValue_!=0.0) {
445 DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
446 // current_madds++;
447 }
448
449 if (fabs(DV[i]) > MaxDiagonalValue) {
450 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
451 else DV[i] = MinDiagonalValue;
452 }
453 else
454 DV[i] = 1.0/DV[i]; // Invert diagonal value
455
456 for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
457
458 if (NumU) {
459 IFPACK_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI)); // Replace current row of L and U
460 }
461
462 // Reset column flags
463 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
464 }
465
466 // Validate that the L and U factors are actually lower and upper triangular
467
468 if (!L_->LowerTriangular())
469 IFPACK_CHK_ERR(-4);
470 if (!U_->UpperTriangular())
471 IFPACK_CHK_ERR(-4);
472
473#ifdef IFPACK_FLOPCOUNTERS
474 // Add up flops
475
476 double current_flops = 2 * current_madds;
477 double total_flops = 0;
478
479 IFPACK_CHK_ERR(Graph().L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1)); // Get total madds across all PEs
480
481 // Now count the rest
482 total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
483 total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
484 if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
485
486 // add to this object's counter
487 ComputeFlops_ += total_flops;
488#endif
489
490 IsComputed_ = true;
491 NumCompute_++;
493
494 return(ierr);
495
496}
497
498//=============================================================================
499// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
500int Ifpack_ILU::Solve(bool Trans, const Epetra_MultiVector& X,
501 Epetra_MultiVector& Y) const
502{
503
504#ifdef IFPACK_TEUCHOS_TIME_MONITOR
505 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_ILU::ApplyInverse - Solve");
506#endif
507
508 // in this function the overlap is always zero
509 bool Upper = true;
510 bool Lower = false;
511 bool UnitDiagonal = true;
512
513 if (!Trans) {
514
515 IFPACK_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, X, Y));
516 // y = D*y (D_ has inverse of diagonal)
517 IFPACK_CHK_ERR(Y.Multiply(1.0, *D_, Y, 0.0));
518 // Solve Uy = y
519 IFPACK_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, Y, Y));
520 }
521 else {
522 // Solve Uy = y
523 IFPACK_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, X, Y));
524 // y = D*y (D_ has inverse of diagonal)
525 IFPACK_CHK_ERR(Y.Multiply(1.0, *D_, Y, 0.0));
526 IFPACK_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, Y, Y));
527 }
528
529
530 return(0);
531}
532
533//=============================================================================
534// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
536 Epetra_MultiVector& Y) const
537{
538
539#ifdef IFPACK_TEUCHOS_TIME_MONITOR
540 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_ILU::Multiply");
541#endif
542
543 if (!IsComputed())
544 IFPACK_CHK_ERR(-3);
545
546 if (!Trans) {
547 IFPACK_CHK_ERR(U_->Multiply(Trans, X, Y));
548 // Y1 = Y1 + X1 (account for implicit unit diagonal)
549 IFPACK_CHK_ERR(Y.Update(1.0, X, 1.0));
550 // y = D*y (D_ has inverse of diagonal)
551 IFPACK_CHK_ERR(Y.ReciprocalMultiply(1.0, *D_, Y, 0.0));
552 Epetra_MultiVector Y1temp(Y); // Need a temp copy of Y1
553 IFPACK_CHK_ERR(L_->Multiply(Trans, Y1temp, Y));
554 // (account for implicit unit diagonal)
555 IFPACK_CHK_ERR(Y.Update(1.0, Y1temp, 1.0));
556 }
557 else {
558
559 IFPACK_CHK_ERR(L_->Multiply(Trans, X, Y));
560 // Y1 = Y1 + X1 (account for implicit unit diagonal)
561 IFPACK_CHK_ERR(Y.Update(1.0, X, 1.0));
562 // y = D*y (D_ has inverse of diagonal)
563 IFPACK_CHK_ERR(Y.ReciprocalMultiply(1.0, *D_, Y, 0.0));
564 Epetra_MultiVector Y1temp(Y); // Need a temp copy of Y1
565 IFPACK_CHK_ERR(U_->Multiply(Trans, Y1temp, Y));
566 // (account for implicit unit diagonal)
567 IFPACK_CHK_ERR(Y.Update(1.0, Y1temp, 1.0));
568 }
569
570 return(0);
571}
572
573//=============================================================================
574// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
576 Epetra_MultiVector& Y) const
577{
578
579#ifdef IFPACK_TEUCHOS_TIME_MONITOR
580 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_ILU::ApplyInverse");
581#endif
582
583 if (!IsComputed())
584 IFPACK_CHK_ERR(-3);
585
586 if (X.NumVectors() != Y.NumVectors())
587 IFPACK_CHK_ERR(-2);
588
590
591 // AztecOO gives X and Y pointing to the same memory location,
592 // need to create an auxiliary vector, Xcopy
593 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
594 if (X.Pointers()[0] == Y.Pointers()[0])
595 Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
596 else
597 Xcopy = Teuchos::rcp( &X, false );
598
600
601 // approx is the number of nonzeros in L and U
602#ifdef IFPACK_FLOPCOUNTERS
603 ApplyInverseFlops_ += X.NumVectors() * 4 *
604 (L_->NumGlobalNonzeros() + U_->NumGlobalNonzeros());
605#endif
606
609
610 return(0);
611
612}
613
614//=============================================================================
616 const int MaxIters, const double Tol,
617 Epetra_RowMatrix* Matrix_in)
618{
619
620#ifdef IFPACK_TEUCHOS_TIME_MONITOR
621 TEUCHOS_FUNC_TIME_MONITOR("Ifpack_ILU::Condest");
622#endif
623
624 if (!IsComputed()) // cannot compute right now
625 return(-1.0);
626
627 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
628
629 return(Condest_);
630}
631
632//=============================================================================
633std::ostream&
634Ifpack_ILU::Print(std::ostream& os) const
635{
636 using std::endl;
637
638 if (!Comm().MyPID()) {
639 os << endl;
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;
648 if (IsComputed_) {
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;
653 }
654 os << endl;
655 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
656 os << "----- ------- -------------- ------------ --------" << endl;
657 os << "Initialize() " << std::setw(5) << NumInitialize()
658 << " " << std::setw(15) << InitializeTime()
659 << " 0.0 0.0" << endl;
660 os << "Compute() " << std::setw(5) << NumCompute()
661 << " " << std::setw(15) << ComputeTime()
662 << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
663 if (ComputeTime() != 0.0)
664 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
665 else
666 os << " " << std::setw(15) << 0.0 << endl;
667 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
668 << " " << std::setw(15) << ApplyInverseTime()
669 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
670 if (ApplyInverseTime() != 0.0)
671 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
672 else
673 os << " " << std::setw(15) << 0.0 << endl;
674 os << "================================================================================" << endl;
675 os << endl;
676 }
677
678 return(os);
679}
#define EPETRA_SGN(x)
const double Epetra_MinDouble
Copy
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
int MaxNumEntries() const
int NumVectors() const
double ** Pointers() const
void ResetStartTime(void)
double ElapsedTime(void) const
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition Ifpack_ILU.h:254
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
Definition Ifpack_ILU.h:103
double ComputeFlops_
Contains the number of flops for Compute().
Definition Ifpack_ILU.h:443
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Definition Ifpack_ILU.h:282
const Epetra_Map * L_RangeMap_
Definition Ifpack_ILU.h:398
double Athresh_
absolute threshold
Definition Ifpack_ILU.h:417
double ComputeTime_
Contains the time for all successful calls to Compute().
Definition Ifpack_ILU.h:439
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
Contains the U factors.
Definition Ifpack_ILU.h:403
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Definition Ifpack_ILU.h:277
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition Ifpack_ILU.h:260
double RelaxValue_
Relaxation value.
Definition Ifpack_ILU.h:415
Teuchos::RefCountPtr< Ifpack_IlukGraph > Graph_
Definition Ifpack_ILU.h:392
int SetParameters(Teuchos::ParameterList &parameterlist)
Set parameters using a Teuchos::ParameterList object.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition Ifpack_ILU.h:248
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
Definition Ifpack_ILU.h:227
const Epetra_Map * U_DomainMap_
Definition Ifpack_ILU.h:397
int Compute()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
int Solve(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_ILU forward/back solve on a Epetra_MultiVector X in Y.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition Ifpack_ILU.h:224
virtual std::ostream & Print(std::ostream &os) const
Prints on stream basic information about this object.
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the Epetra_RowMatrix to factorize.
Definition Ifpack_ILU.h:391
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
Definition Ifpack_ILU.h:435
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
Contains the L factors.
Definition Ifpack_ILU.h:401
int ComputeSetup()
double AbsoluteThreshold() const
Get absolute threshold value.
Definition Ifpack_ILU.h:330
bool IsComputed_
If true, the preconditioner has been successfully computed.
Definition Ifpack_ILU.h:427
double RelativeThreshold() const
Get relative threshold value.
Definition Ifpack_ILU.h:333
int NumMyRows() const
Returns the number of local matrix rows.
Definition Ifpack_ILU.h:358
int NumMyCols() const
Returns the number of local matrix columns.
Definition Ifpack_ILU.h:361
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition Ifpack_ILU.h:266
bool UseTranspose() const
Returns the current UseTranspose setting.
Definition Ifpack_ILU.h:215
int NumCompute_
Contains the number of successful call to Compute().
Definition Ifpack_ILU.h:433
long long NumGlobalNonzeros64() const
Definition Ifpack_ILU.h:353
double Condest() const
Returns the computed estimated condition number, or -1.0 if not computed.
Definition Ifpack_ILU.h:181
Teuchos::RefCountPtr< Epetra_CrsGraph > CrsGraph_
Definition Ifpack_ILU.h:393
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition Ifpack_ILU.h:236
int LevelOfFill_
Level of fill.
Definition Ifpack_ILU.h:423
double Rthresh_
relative threshold
Definition Ifpack_ILU.h:419
Teuchos::RefCountPtr< Epetra_Vector > D_
Diagonal of factors.
Definition Ifpack_ILU.h:407
void Destroy()
Destroys all internal data.
const Ifpack_IlukGraph & Graph() const
Returns the address of the Ifpack_IlukGraph associated with this factored matrix.
Definition Ifpack_ILU.h:379
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.
int Initialize()
Initialize the preconditioner, does not touch matrix values.
double Condest_
condition number estimate
Definition Ifpack_ILU.h:421
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition Ifpack_ILU.h:242
double RelaxValue() const
Get ILU(k) relaxation parameter.
Definition Ifpack_ILU.h:327
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
char Label_[160]
Label of this object.
Definition Ifpack_ILU.h:429
Epetra_Time Time_
Used for timing issues.
Definition Ifpack_ILU.h:447
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
Definition Ifpack_ILU.h:120
double InitializeTime_
Contains the time for all successful calls to Initialize().
Definition Ifpack_ILU.h:437
int LevelOfFill() const
Returns the level of fill.
Definition Ifpack_ILU.h:324
Ifpack_ILU(Epetra_RowMatrix *A)
Constructor.
const char * Label() const
Returns a character string describing the operator.
Definition Ifpack_ILU.h:199
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Definition Ifpack_ILU.h:441
int NumInitialize_
Contains the number of successful calls to Initialize().
Definition Ifpack_ILU.h:431
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
Definition Ifpack_ILU.h:425
int NumMyDiagonals_
Definition Ifpack_ILU.h:410
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
Definition Ifpack_ILU.h:445
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
#define false