Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_CrsRick.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_CrsRick.h"
44 #include "Epetra_Comm.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_CrsGraph.h"
47 #include "Epetra_CrsMatrix.h"
48 #include "Epetra_Vector.h"
49 #include "Epetra_MultiVector.h"
50 
51 #include <Teuchos_ParameterList.hpp>
52 #include <ifp_parameters.h>
53 
54 //==============================================================================
56  : A_(A),
57  Graph_(Graph),
58  UseTranspose_(false),
59  Allocated_(false),
60  ValuesInitialized_(false),
61  Factored_(false),
62  RelaxValue_(0.0),
63  Condest_(-1.0),
64  Athresh_(0.0),
65  Rthresh_(1.0),
66  OverlapX_(0),
67  OverlapY_(0),
68  OverlapMode_(Zero)
69 {
70  int ierr = Allocate();
71 }
72 
73 //==============================================================================
75  : A_(FactoredMatrix.A_),
76  Graph_(FactoredMatrix.Graph_),
77  UseTranspose_(FactoredMatrix.UseTranspose_),
78  Allocated_(FactoredMatrix.Allocated_),
79  ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
80  Factored_(FactoredMatrix.Factored_),
81  RelaxValue_(FactoredMatrix.RelaxValue_),
82  Condest_(FactoredMatrix.Condest_),
83  Athresh_(FactoredMatrix.Athresh_),
84  Rthresh_(FactoredMatrix.Rthresh_),
85  OverlapX_(0),
86  OverlapY_(0),
87  OverlapMode_(FactoredMatrix.OverlapMode_)
88 {
89  U_ = new Epetra_CrsMatrix(FactoredMatrix.U());
91 
92 }
93 
94 //==============================================================================
96 
97  // Allocate Epetra_CrsMatrix using ILUK graphs
100 
101 
102  SetAllocated(true);
103  return(0);
104 }
105 //==============================================================================
107 
108 
109  delete U_;
110  delete D_; // Diagonal is stored separately. We store the inverse.
111 
112  if (OverlapX_!=0) delete OverlapX_;
113  if (OverlapY_!=0) delete OverlapY_;
114 
115  ValuesInitialized_ = false;
116  Factored_ = false;
117  Allocated_ = false;
118 }
119 
120 //==========================================================================
122  bool cerr_warning_if_unused)
123 {
124  Ifpack::param_struct params;
128  params.overlap_mode = OverlapMode_;
129 
130  Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
131 
135  OverlapMode_ = params.overlap_mode;
136 
137  return(0);
138 }
139 
140 //==========================================================================
142 
143  // if (!Allocated()) return(-1); // This test is not needed at this time. All constructors allocate.
144 
145  int ierr = 0;
146  int i, j;
147  int * InI=0, * LI=0, * UI = 0;
148  double * InV=0, * LV=0, * UV = 0;
149  int NumIn, NumL, NumU;
150  bool DiagFound;
151  int NumNonzeroDiags = 0;
152 
153  Epetra_CrsMatrix * OverlapA = (Epetra_CrsMatrix *) &A_;
154 
156 
157  OverlapA = new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph());
158  OverlapA->Import(A_, *Graph_.OverlapImporter(), Insert);
159  }
160 
161  // Get Maximun Row length
162  int MaxNumEntries = OverlapA->MaxNumEntries();
163 
164  InI = new int[MaxNumEntries]; // Allocate temp space
165  UI = new int[MaxNumEntries];
166  InV = new double[MaxNumEntries];
167  UV = new double[MaxNumEntries];
168 
169  double *DV;
170  ierr = D_->ExtractView(&DV); // Get view of diagonal
171 
172 
173  // First we copy the user's matrix into diagonal vector and U, regardless of fill level
174 
175  for (i=0; i< NumMyRows(); i++) {
176 
177  OverlapA->ExtractMyRowCopy(i, MaxNumEntries, NumIn, InV, InI); // Get Values and Indices
178 
179  // Split into L and U (we don't assume that indices are ordered).
180 
181  NumL = 0;
182  NumU = 0;
183  DiagFound = false;
184 
185  for (j=0; j< NumIn; j++) {
186  int k = InI[j];
187 
188  if (k==i) {
189  DiagFound = true;
190  DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
191  }
192 
193  else if (k < 0) return(-1); // Out of range
194  else if (k<NumMyRows()) {
195  UI[NumU] = k;
196  UV[NumU] = InV[j];
197  NumU++;
198  }
199  }
200 
201  // Check in things for this row of L and U
202 
203  if (DiagFound) NumNonzeroDiags++;
204  if (NumU) U_->ReplaceMyValues(i, NumU, UV, UI);
205 
206  }
207 
208  delete [] UI;
209  delete [] UV;
210  delete [] InI;
211  delete [] InV;
212 
213  if (Graph_.LevelOverlap()>0 && Graph_.L_Graph().DomainMap().DistributedGlobal()) delete OverlapA;
214 
215 
216  U_->FillComplete();
217 
218  // At this point L and U have the values of A in the structure of L and U, and diagonal vector D
219 
220  SetValuesInitialized(true);
221  SetFactored(false);
222 
223  int TotalNonzeroDiags = 0;
224  Graph_.U_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1);
225  if (Graph_.LevelOverlap()==0 &&
226  ((TotalNonzeroDiags!=NumGlobalRows()) ||
227  (TotalNonzeroDiags!=NumGlobalDiagonals()))) ierr = 1;
228  if (NumNonzeroDiags != NumMyDiagonals()) ierr = 1; // Diagonals are not right, warn user
229 
230  return(ierr);
231 }
232 
233 //==========================================================================
235 
236  // if (!Allocated()) return(-1); // This test is not needed at this time. All constructors allocate.
237  if (!ValuesInitialized()) return(-2); // Must have values initialized.
238  if (Factored()) return(-3); // Can't have already computed factors.
239 
240  SetValuesInitialized(false);
241 
242  // MinMachNum should be officially defined, for now pick something a little
243  // bigger than IEEE underflow value
244 
245  double MinDiagonalValue = Epetra_MinDouble;
246  double MaxDiagonalValue = 1.0/MinDiagonalValue;
247 
248  int ierr = 0;
249  int i, j, k;
250  int * UI = 0;
251  double * UV = 0;
252  int NumIn, NumU;
253 
254  // Get Maximun Row length
255  int MaxNumEntries = U_->MaxNumEntries() + 1;
256 
257  int * InI = new int[MaxNumEntries]; // Allocate temp space
258  double * InV = new double[MaxNumEntries];
259  int * colflag = new int[NumMyCols()];
260 
261  double *DV;
262  ierr = D_->ExtractView(&DV); // Get view of diagonal
263 
264 #ifdef IFPACK_FLOPCOUNTERS
265  int current_madds = 0; // We will count multiply-add as they happen
266 #endif
267 
268  // Now start the factorization.
269 
270  // Need some integer workspace and pointers
271  int NumUU;
272  int * UUI;
273  double * UUV;
274  for (j=0; j<NumMyCols(); j++) colflag[j] = - 1;
275 
276  for(i=0; i<NumMyRows(); i++) {
277 
278  // Fill InV, InI with current row of L, D and U combined
279 
280  NumIn = MaxNumEntries;
281  IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, InV, InI)==0);
282  LV = InV;
283  LI = InI;
284 
285  InV[NumL] = DV[i]; // Put in diagonal
286  InI[NumL] = i;
287 
288  IFPACK_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, InV+NumL+1, InI+NumL+1));
289  NumIn = NumL+NumU+1;
290  UV = InV+NumL+1;
291  UI = InI+NumL+1;
292 
293  // Set column flags
294  for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
295 
296  double diagmod = 0.0; // Off-diagonal accumulator
297 
298  for (int jj=0; jj<NumL; jj++) {
299  j = InI[jj];
300  double multiplier = InV[jj]; // current_mults++;
301 
302  InV[jj] *= DV[j];
303 
304  IFPACK_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above
305 
306  if (RelaxValue_==0.0) {
307  for (k=0; k<NumUU; k++) {
308  int kk = colflag[UUI[k]];
309  if (kk>-1) {
310  InV[kk] -= multiplier*UUV[k];
311 #ifdef IFPACK_FLOPCOUNTERS
312  current_madds++;
313 #endif
314 #endif
315  }
316  }
317  }
318  else {
319  for (k=0; k<NumUU; k++) {
320  int kk = colflag[UUI[k]];
321  if (kk>-1) InV[kk] -= multiplier*UUV[k];
322  else diagmod -= multiplier*UUV[k];
323 #ifdef IFPACK_FLOPCOUNTERS
324  current_madds++;
325 #endif
326  }
327  }
328  }
329  if (NumL)
330  IFPACK_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI)); // Replace current row of L
331 
332  DV[i] = InV[NumL]; // Extract Diagonal value
333 
334  if (RelaxValue_!=0.0) {
335  DV[i] += RelaxValue_*diagmod; // Add off diagonal modifications
336  // current_madds++;
337  }
338 
339  if (fabs(DV[i]) > MaxDiagonalValue) {
340  if (DV[i] < 0) DV[i] = - MinDiagonalValue;
341  else DV[i] = MinDiagonalValue;
342  }
343  else
344  DV[i] = 1.0/DV[i]; // Invert diagonal value
345 
346  for (j=0; j<NumU; j++) UV[j] *= DV[i]; // Scale U by inverse of diagonal
347 
348  if (NumU)
349  IFPACK_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI)); // Replace current row of L and U
350 
351 
352  // Reset column flags
353  for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
354  }
355 
356 
357 #ifdef IFPACK_FLOPCOUNTERS
358  // Add up flops
359 
360  double current_flops = 2 * current_madds;
361  double total_flops = 0;
362 
363  Graph_.L_Graph().RowMap().Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs
364 
365  // Now count the rest
366  total_flops += (double) L_->NumGlobalNonzeros(); // Accounts for multiplier above
367  total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
368  if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->GlobalLength(); // Accounts for relax update of diag
369 
370  UpdateFlops(total_flops); // Update flop count
371 #endif
372 
373  delete [] InI;
374  delete [] InV;
375  delete [] colflag;
376 
377  SetFactored(true);
378 
379  return(ierr);
380 
381 }
382 
383 //=============================================================================
384 int Ifpack_CrsRick::Solve(bool Trans, const Epetra_Vector& X,
385  Epetra_Vector& Y) const {
386 //
387 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for a single RHS
388 //
389 
390  bool Upper = true;
391  bool Lower = false;
392  bool UnitDiagonal = true;
393 
394  Epetra_Vector * X1 = (Epetra_Vector *) &X;
395  Epetra_Vector * Y1 = (Epetra_Vector *) &Y;
396 
398  if (OverlapX_==0) { // Need to allocate space for overlap X and Y
401  }
402  OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
403  X1 = (Epetra_Vector *) OverlapX_;
404  Y1 = (Epetra_Vector *) OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
405  }
406 
407 #ifdef IFPACK_FLOPCOUNTERS
408  Epetra_Flops * counter = this->GetFlopCounter();
409  if (counter!=0) {
410  L_->SetFlopCounter(*counter);
411  Y1->SetFlopCounter(*counter);
412  U_->SetFlopCounter(*counter);
413  }
414 #endif
415 
416  if (!Trans) {
417 
418  L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
419  Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
420  U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1); // Solve Uy = y
421  }
422  else
423  {
424  U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1); // Solve Uy = y
425  Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
426  L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
427 
428  }
429 
430  // Export computed Y values as directed
433  return(0);
434 }
435 
436 
437 //=============================================================================
438 int Ifpack_CrsRick::Solve(bool Trans, const Epetra_MultiVector& X,
439  Epetra_MultiVector& Y) const {
440 //
441 // This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
442 //
443 
444  if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
445 
446  bool Upper = true;
447  bool Lower = false;
448  bool UnitDiagonal = true;
449 
452 
454  // Make sure the number of vectors in the multivector is the same as before.
455  if (OverlapX_!=0) {
456  if (OverlapX_->NumVectors()!=X.NumVectors()) {
457  delete OverlapX_; OverlapX_ = 0;
458  delete OverlapY_; OverlapY_ = 0;
459  }
460  }
461  if (OverlapX_==0) { // Need to allocate space for overlap X and Y
462  OverlapX_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), X.NumVectors());
463  OverlapY_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), Y.NumVectors());
464  }
465  OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
466  X1 = OverlapX_;
467  Y1 = OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
468  }
469 
470 #ifdef IFPACK_FLOPCOUNTERS
471  Epetra_Flops * counter = this->GetFlopCounter();
472  if (counter!=0) {
473  L_->SetFlopCounter(*counter);
474  Y1->SetFlopCounter(*counter);
475  U_->SetFlopCounter(*counter);
476  }
477 #endif
478 
479  if (!Trans) {
480 
481  L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
482  Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
483  U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1); // Solve Uy = y
484  }
485  else
486  {
487  U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1); // Solve Uy = y
488  Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
489  L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
490 
491  }
492 
493  // Export computed Y values as directed
496  return(0);
497 }
498 //=============================================================================
500  Epetra_MultiVector& Y) const {
501 //
502 // This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
503 //
504 
505  if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size
506 
507  bool Upper = true;
508  bool Lower = false;
509  bool UnitDiagonal = true;
510 
513 
515  // Make sure the number of vectors in the multivector is the same as before.
516  if (OverlapX_!=0) {
517  if (OverlapX_->NumVectors()!=X.NumVectors()) {
518  delete OverlapX_; OverlapX_ = 0;
519  delete OverlapY_; OverlapY_ = 0;
520  }
521  }
522  if (OverlapX_==0) { // Need to allocate space for overlap X and Y
523  OverlapX_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), X.NumVectors());
524  OverlapY_ = new Epetra_MultiVector(Graph_.OverlapGraph()->RowMap(), Y.NumVectors());
525  }
526  OverlapX_->Import(X,*Graph_.OverlapImporter(), Insert); // Import X values for solve
527  X1 = OverlapX_;
528  Y1 = OverlapY_; // Set pointers for X1 and Y1 to point to overlap space
529  }
530 
531 #ifdef IFPACK_FLOPCOUNTERS
532  Epetra_Flops * counter = this->GetFlopCounter();
533  if (counter!=0) {
534  L_->SetFlopCounter(*counter);
535  Y1->SetFlopCounter(*counter);
536  U_->SetFlopCounter(*counter);
537  }
538 #endif
539 
540  if (Trans) {
541 
542  L_->Multiply(Trans, *X1, *Y1);
543  Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
544  Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
545  Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
546  U_->Multiply(Trans, Y1temp, *Y1);
547  Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
548  }
549  else {
550  U_->Multiply(Trans, *X1, *Y1); //
551  Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
552  Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
553  Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
554  L_->Multiply(Trans, Y1temp, *Y1);
555  Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
556  }
557 
558  // Export computed Y values as directed
561  return(0);
562 }
563 //=============================================================================
564 int Ifpack_CrsRick::Condest(bool Trans, double & ConditionNumberEstimate) const {
565 
566  if (Condest_>=0.0) {
567  ConditionNumberEstimate = Condest_;
568  return(0);
569  }
570  // Create a vector with all values equal to one
571  Epetra_Vector Ones(A_.RowMap());
572  Epetra_Vector OnesResult(Ones);
573  Ones.PutScalar(1.0);
574 
575  EPETRA_CHK_ERR(Solve(Trans, Ones, OnesResult)); // Compute the effect of the solve on the vector of ones
576  EPETRA_CHK_ERR(OnesResult.Abs(OnesResult)); // Make all values non-negative
577  EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate)); // Get the maximum value across all processors
578  Condest_ = ConditionNumberEstimate; // Save value for possible later calls
579  return(0);
580 }
581 //=============================================================================
582 // Non-member functions
583 
584 std::ostream& operator << (std::ostream& os, const Ifpack_CrsRick& A)
585 {
586 /* Epetra_fmtflags olda = os.setf(ios::right,ios::adjustfield);
587  Epetra_fmtflags oldf = os.setf(ios::scientific,ios::floatfield);
588  int oldp = os.precision(12); */
589  int LevelFill = A.Graph().LevelFill();
590  int LevelOverlap = A.Graph().LevelOverlap();
591  Epetra_CrsMatrix & L = (Epetra_CrsMatrix &) A.L();
592  Epetra_CrsMatrix & U = (Epetra_CrsMatrix &) A.U();
593  Epetra_Vector & D = (Epetra_Vector &) A.D();
594 
595  os.width(14);
596  os << endl;
597  os << " Level of Fill = "; os << LevelFill;
598  os << endl;
599  os.width(14);
600  os << " Level of Overlap = "; os << LevelOverlap;
601  os << endl;
602 
603  os.width(14);
604  os << " Lower Triangle = ";
605  os << endl;
606  os << L; // Let Epetra_CrsMatrix handle the rest.
607  os << endl;
608 
609  os.width(14);
610  os << " Inverse of Diagonal = ";
611  os << endl;
612  os << D; // Let Epetra_Vector handle the rest.
613  os << endl;
614 
615  os.width(14);
616  os << " Upper Triangle = ";
617  os << endl;
618  os << U; // Let Epetra_CrsMatrix handle the rest.
619  os << endl;
620 
621  // Reset os flags
622 
623 /* os.setf(olda,ios::adjustfield);
624  os.setf(oldf,ios::floatfield);
625  os.precision(oldp); */
626 
627  return os;
628 }
bool DistributedGlobal() const
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
Epetra_CrsMatrix * U_
const Epetra_CrsMatrix & U() const
Returns the address of the U factor associated with this factored matrix.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
double double_params[FIRST_INT_PARAM]
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Epetra_MultiVector * OverlapY_
int SetAllocated(bool Flag)
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
Epetra_MultiVector * OverlapX_
#define EPETRA_CHK_ERR(a)
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int NumMyRows() const
Returns the number of local matrix rows.
const Epetra_BlockMap & DomainMap() const
const Ifpack_IlukGraph & Graph() const
Returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
int InitValues()
Initialize L and U with values from user matrix A.
int FillComplete(bool OptimizeDataStorage=true)
const Epetra_CrsMatrix & A_
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
const Epetra_Map & RowMap() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
const Ifpack_IlukGraph & Graph_
int SetParameters(const Teuchos::ParameterList &parameterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
int NumGlobalRows() const
Returns the number of global matrix rows.
Epetra_CombineMode overlap_mode
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
#define EPETRA_SGN(x)
int MaxNumEntries() const
virtual int NumGlobalDiagonals() const
Returns the number of diagonal entries found in the global input graph.
const Epetra_BlockMap & RowMap() const
Epetra_Vector * D_
adjacency_list< vecS, vecS, undirectedS, no_property, property< edge_weight_t, double > > Graph
Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a giv...
Epetra_CombineMode OverlapMode_
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
#define false
std::ostream & operator<<(std::ostream &os, const Ifpack_Container &obj)
const Epetra_Comm & Comm() const
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.
Ifpack_CrsRick(const Epetra_CrsMatrix &A, const Ifpack_IlukGraph &Graph)
Ifpack_CrsRick constuctor with variable number of indices per row.
const double Epetra_MinDouble
virtual int NumMyDiagonals() const
Returns the number of diagonal entries found in the local input graph.
virtual ~Ifpack_CrsRick()
Ifpack_CrsRick Destructor.
void SetFactored(bool Flag)
int NumMyCols() const
Returns the number of local matrix columns.
int Solve(bool Trans, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Ifpack_CrsRick forward/back solve on a Epetra_Vector x in y...
void SetValuesInitialized(bool Flag)
int Factor()
Compute IC factor L using the specified graph, diagonal perturbation thresholds and relaxation parame...
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of multiplying U, D and U^T in that order on an Epetra_MultiVector X in Y...
void set_parameters(const Teuchos::ParameterList &parameterlist, param_struct &params, bool cerr_warning_if_unused)
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
#define IFPACK_CHK_ERR(ifpack_err)
int ExtractView(double **V) const
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local IC set of factors...
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.