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" 51 #include <Teuchos_ParameterList.hpp> 60 ValuesInitialized_(
false),
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_),
87 OverlapMode_(FactoredMatrix.OverlapMode_)
89 U_ =
new Epetra_CrsMatrix(FactoredMatrix.
U());
122 bool cerr_warning_if_unused)
147 int * InI=0, * LI=0, * UI = 0;
148 double * InV=0, * LV=0, * UV = 0;
149 int NumIn, NumL, NumU;
151 int NumNonzeroDiags = 0;
153 Epetra_CrsMatrix * OverlapA = (Epetra_CrsMatrix *) &
A_;
162 int MaxNumEntries = OverlapA->MaxNumEntries();
164 InI =
new int[MaxNumEntries];
165 UI =
new int[MaxNumEntries];
166 InV =
new double[MaxNumEntries];
167 UV =
new double[MaxNumEntries];
170 ierr =
D_->ExtractView(&DV);
177 OverlapA->ExtractMyRowCopy(i, MaxNumEntries, NumIn, InV, InI);
185 for (j=0; j< NumIn; j++) {
193 else if (k < 0)
return(-1);
203 if (DiagFound) NumNonzeroDiags++;
204 if (NumU)
U_->ReplaceMyValues(i, NumU, UV, UI);
223 int TotalNonzeroDiags = 0;
224 Graph_.
U_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1);
245 double MinDiagonalValue = Epetra_MinDouble;
246 double MaxDiagonalValue = 1.0/MinDiagonalValue;
255 int MaxNumEntries =
U_->MaxNumEntries() + 1;
257 int * InI =
new int[MaxNumEntries];
258 double * InV =
new double[MaxNumEntries];
262 ierr =
D_->ExtractView(&DV);
264 #ifdef IFPACK_FLOPCOUNTERS 265 int current_madds = 0;
274 for (j=0; j<
NumMyCols(); j++) colflag[j] = - 1;
280 NumIn = MaxNumEntries;
281 IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, InV, InI)==0);
288 IFPACK_CHK_ERR(
U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, InV+NumL+1, InI+NumL+1));
294 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
296 double diagmod = 0.0;
298 for (
int jj=0; jj<NumL; jj++) {
300 double multiplier = InV[jj];
307 for (k=0; k<NumUU; k++) {
308 int kk = colflag[UUI[k]];
310 InV[kk] -= multiplier*UUV[k];
311 #ifdef IFPACK_FLOPCOUNTERS 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 339 if (fabs(DV[i]) > MaxDiagonalValue) {
340 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
341 else DV[i] = MinDiagonalValue;
346 for (j=0; j<NumU; j++) UV[j] *= DV[i];
353 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
357 #ifdef IFPACK_FLOPCOUNTERS 360 double current_flops = 2 * current_madds;
361 double total_flops = 0;
363 Graph_.
L_Graph().RowMap().Comm().SumAll(¤t_flops, &total_flops, 1);
366 total_flops += (double) L_->NumGlobalNonzeros();
367 total_flops += (double)
D_->GlobalLength();
368 if (
RelaxValue_!=0.0) total_flops += 2 * (double)
D_->GlobalLength();
370 UpdateFlops(total_flops);
385 Epetra_Vector& Y)
const {
392 bool UnitDiagonal =
true;
394 Epetra_Vector * X1 = (Epetra_Vector *) &X;
395 Epetra_Vector * Y1 = (Epetra_Vector *) &Y;
407 #ifdef IFPACK_FLOPCOUNTERS 408 Epetra_Flops * counter = this->GetFlopCounter();
410 L_->SetFlopCounter(*counter);
411 Y1->SetFlopCounter(*counter);
412 U_->SetFlopCounter(*counter);
418 L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
419 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
420 U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1);
424 U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1);
425 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
426 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
439 Epetra_MultiVector& Y)
const {
444 if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1);
448 bool UnitDiagonal =
true;
450 Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
451 Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
456 if (
OverlapX_->NumVectors()!=X.NumVectors()) {
470 #ifdef IFPACK_FLOPCOUNTERS 471 Epetra_Flops * counter = this->GetFlopCounter();
473 L_->SetFlopCounter(*counter);
474 Y1->SetFlopCounter(*counter);
475 U_->SetFlopCounter(*counter);
481 L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
482 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
483 U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1);
487 U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1);
488 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
489 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
500 Epetra_MultiVector& Y)
const {
505 if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1);
509 bool UnitDiagonal =
true;
511 Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
512 Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
517 if (
OverlapX_->NumVectors()!=X.NumVectors()) {
531 #ifdef IFPACK_FLOPCOUNTERS 532 Epetra_Flops * counter = this->GetFlopCounter();
534 L_->SetFlopCounter(*counter);
535 Y1->SetFlopCounter(*counter);
536 U_->SetFlopCounter(*counter);
542 L_->Multiply(Trans, *X1, *Y1);
543 Y1->Update(1.0, *X1, 1.0);
544 Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0);
545 Epetra_MultiVector Y1temp(*Y1);
546 U_->Multiply(Trans, Y1temp, *Y1);
547 Y1->Update(1.0, Y1temp, 1.0);
550 U_->Multiply(Trans, *X1, *Y1);
551 Y1->Update(1.0, *X1, 1.0);
552 Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0);
553 Epetra_MultiVector Y1temp(*Y1);
554 L_->Multiply(Trans, Y1temp, *Y1);
555 Y1->Update(1.0, Y1temp, 1.0);
571 Epetra_Vector Ones(
A_.RowMap());
572 Epetra_Vector OnesResult(Ones);
575 EPETRA_CHK_ERR(
Solve(Trans, Ones, OnesResult));
576 EPETRA_CHK_ERR(OnesResult.Abs(OnesResult));
577 EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate));
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();
597 os <<
" Level of Fill = "; os << LevelFill;
600 os <<
" Level of Overlap = "; os << LevelOverlap;
604 os <<
" Lower Triangle = ";
610 os <<
" Inverse of Diagonal = ";
616 os <<
" Upper Triangle = ";
double double_params[FIRST_INT_PARAM]
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
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.
Epetra_MultiVector * OverlapX_
int NumMyRows() const
Returns the number of local matrix rows.
std::ostream & operator<<(std::ostream &os, const Ifpack_CrsRick &A)
<< operator will work for Ifpack_CrsRick.
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int InitValues()
Initialize L and U with values from user matrix A.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
const Epetra_CrsMatrix & A_
const Epetra_CrsMatrix & U() const
Returns the address of the U factor associated with this factored matrix.
const Ifpack_IlukGraph & Graph_
int SetParameters(const Teuchos::ParameterList ¶meterlist, 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.
Epetra_CombineMode overlap_mode
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_
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
virtual int NumGlobalDiagonals() const
Returns the number of diagonal entries found in the global input graph.
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...
Ifpack_CrsRick(const Epetra_CrsMatrix &A, const Ifpack_IlukGraph &Graph)
Ifpack_CrsRick constuctor with variable number of indices per row.
virtual int NumMyDiagonals() const
Returns the number of diagonal entries found in the local input graph.
int NumMyCols() const
Returns the number of local matrix columns.
virtual ~Ifpack_CrsRick()
Ifpack_CrsRick Destructor.
void SetFactored(bool Flag)
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 ¶meterlist, param_struct ¶ms, bool cerr_warning_if_unused)
#define IFPACK_CHK_ERR(ifpack_err)
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local IC set of factors...
int NumGlobalRows() const
Returns the number of global matrix rows.
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.