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" 51 #include <Teuchos_ParameterList.hpp> 52 #include <ifp_parameters.h> 60 ValuesInitialized_(false),
70 int ierr = Allocate();
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());
90 D_ =
new Epetra_Vector(Graph_.
L_Graph().RowMap());
95 int Ifpack_CrsRick::Allocate() {
98 U_ =
new Epetra_CrsMatrix(Copy, Graph_.
U_Graph());
99 D_ =
new Epetra_Vector(Graph_.
L_Graph().RowMap());
112 if (OverlapX_!=0)
delete OverlapX_;
113 if (OverlapY_!=0)
delete OverlapY_;
115 ValuesInitialized_ =
false;
122 bool cerr_warning_if_unused)
125 params.double_params[Ifpack::relax_value] = RelaxValue_;
126 params.double_params[Ifpack::absolute_threshold] = Athresh_;
127 params.double_params[Ifpack::relative_threshold] = Rthresh_;
128 params.overlap_mode = OverlapMode_;
130 Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
132 RelaxValue_ = params.double_params[Ifpack::relax_value];
133 Athresh_ = params.double_params[Ifpack::absolute_threshold];
134 Rthresh_ = params.double_params[Ifpack::relative_threshold];
135 OverlapMode_ = params.overlap_mode;
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_;
157 OverlapA =
new Epetra_CrsMatrix(Copy, *Graph_.
OverlapGraph());
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++) {
190 DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_;
193 else if (k < 0)
return(-1);
203 if (DiagFound) NumNonzeroDiags++;
204 if (NumU) U_->ReplaceMyValues(i, NumU, UV, UI);
213 if (Graph_.
LevelOverlap()>0 && Graph_.
L_Graph().DomainMap().DistributedGlobal())
delete OverlapA;
220 SetValuesInitialized(
true);
223 int TotalNonzeroDiags = 0;
224 Graph_.
U_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1);
240 SetValuesInitialized(
false);
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];
304 IFPACK_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI));
306 if (RelaxValue_==0.0) {
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 330 IFPACK_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));
334 if (RelaxValue_!=0.0) {
335 DV[i] += RelaxValue_*diagmod;
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];
349 IFPACK_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI));
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;
399 OverlapX_ =
new Epetra_Vector(Graph_.
OverlapGraph()->RowMap());
400 OverlapY_ =
new Epetra_Vector(Graph_.
OverlapGraph()->RowMap());
403 X1 = (Epetra_Vector *) OverlapX_;
404 Y1 = (Epetra_Vector *) OverlapY_;
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()) {
457 delete OverlapX_; OverlapX_ = 0;
458 delete OverlapY_; OverlapY_ = 0;
462 OverlapX_ =
new Epetra_MultiVector(Graph_.
OverlapGraph()->RowMap(), X.NumVectors());
463 OverlapY_ =
new Epetra_MultiVector(Graph_.
OverlapGraph()->RowMap(), Y.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()) {
518 delete OverlapX_; OverlapX_ = 0;
519 delete OverlapY_; OverlapY_ = 0;
523 OverlapX_ =
new Epetra_MultiVector(Graph_.
OverlapGraph()->RowMap(), X.NumVectors());
524 OverlapY_ =
new Epetra_MultiVector(Graph_.
OverlapGraph()->RowMap(), Y.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);
567 ConditionNumberEstimate = Condest_;
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));
578 Condest_ = ConditionNumberEstimate;
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 = ";
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
int NumMyRows() const
Returns the number of local matrix rows.
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...
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.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
const Epetra_CrsMatrix & U() const
Returns the address of the U factor associated with this factored matrix.
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.
Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a giv...
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.
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...
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.