IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_ILU.cpp
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 
57 using Teuchos::RefCountPtr;
58 using 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 //==============================================================================
92 void Ifpack_ILU::Destroy()
93 {
94  // reset pointers to already allocated stuff
95  U_DomainMap_ = 0;
96  L_RangeMap_ = 0;
97 }
98 
99 //==========================================================================
100 int 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)",
109  LevelOfFill(), RelaxValue(), AbsoluteThreshold(),
110  RelativeThreshold());
111  return(0);
112 }
113 
114 //==========================================================================
115 int Ifpack_ILU::ComputeSetup()
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 
246  Time_.ResetStartTime();
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.
314  Graph_ = rcp(new Ifpack_IlukGraph(*CrsGraph_, LevelOfFill_, 0));
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;
327  NumInitialize_++;
328  InitializeTime_ += Time_.ElapsedTime();
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())
342  IFPACK_CHK_ERR(Initialize());
343 
344  Time_.ResetStartTime();
345  IsComputed_ = false;
346 
347  // convert Matrix() into L and U factors.
348  IFPACK_CHK_ERR(ComputeSetup());
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_++;
492  ComputeTime_ += Time_.ElapsedTime();
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
500 int 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
535 int Ifpack_ILU::Multiply(bool Trans, const Epetra_MultiVector& X,
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 
589  Time_.ResetStartTime();
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 
599  IFPACK_CHK_ERR(Solve(Ifpack_ILU::UseTranspose(), *Xcopy, Y));
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 
607  ++NumApplyInverse_;
608  ApplyInverseTime_ += Time_.ElapsedTime();
609 
610  return(0);
611 
612 }
613 
614 //=============================================================================
615 double Ifpack_ILU::Condest(const Ifpack_CondestType CT,
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 //=============================================================================
633 std::ostream&
634 Ifpack_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 }
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...
Definition: Ifpack_ILU.cpp:575
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack_ILU.h:242
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack_ILU.h:277
double ElapsedTime(void) const
bool UseTranspose() const
Returns the current UseTranspose setting.
Definition: Ifpack_ILU.h:215
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.
Definition: Ifpack_ILU.h:227
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition: Ifpack_ILU.h:224
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack_ILU.h:260
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.
Definition: Ifpack_ILU.h:120
int Compute()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
Definition: Ifpack_ILU.cpp:334
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition: Ifpack_ILU.h:266
int SetParameters(Teuchos::ParameterList &parameterlist)
Set parameters using a Teuchos::ParameterList object.
Definition: Ifpack_ILU.cpp:100
int Initialize()
Initialize the preconditioner, does not touch matrix values.
Definition: Ifpack_ILU.cpp:239
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack_ILU.h:254
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Definition: Ifpack_ILU.h:282
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
Definition: Ifpack_ILU.h:103
double Condest() const
Returns the computed estimated condition number, or -1.0 if not computed.
Definition: Ifpack_ILU.h:181
virtual std::ostream & Print(std::ostream &os) const
Prints on stream basic information about this object.
Definition: Ifpack_ILU.cpp:634
Ifpack_ILU(Epetra_RowMatrix *A)
Constructor.
Definition: Ifpack_ILU.cpp:65
const Epetra_CrsGraph & Graph() const
const char * Label() const
Returns a character string describing the operator.
Definition: Ifpack_ILU.h:199
void ResetStartTime(void)
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition: Ifpack_ILU.h:248
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack_ILU.h:236