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_)
122 bool cerr_warning_if_unused)
126 params.double_params[Ifpack::absolute_threshold] =
Athresh_;
127 params.double_params[Ifpack::relative_threshold] =
Rthresh_;
132 RelaxValue_ = params.double_params[Ifpack::relax_value];
133 Athresh_ = params.double_params[Ifpack::absolute_threshold];
134 Rthresh_ = params.double_params[Ifpack::relative_threshold];
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;
155 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal()) {
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];
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++;
213 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal())
delete OverlapA;
223 int TotalNonzeroDiags = 0;
224 Graph_.U_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1);
225 if (
Graph_.LevelOverlap()==0 &&
246 double MaxDiagonalValue = 1.0/MinDiagonalValue;
257 int * InI =
new int[MaxNumEntries];
258 double * InV =
new double[MaxNumEntries];
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);
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();
392 bool UnitDiagonal =
true;
397 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal()) {
407#ifdef IFPACK_FLOPCOUNTERS
410 L_->SetFlopCounter(*counter);
411 Y1->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);
431 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal())
448 bool UnitDiagonal =
true;
453 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal()) {
470#ifdef IFPACK_FLOPCOUNTERS
473 L_->SetFlopCounter(*counter);
474 Y1->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);
494 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal())
509 bool UnitDiagonal =
true;
514 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal()) {
531#ifdef IFPACK_FLOPCOUNTERS
534 L_->SetFlopCounter(*counter);
535 Y1->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);
547 Y1->Update(1.0, Y1temp, 1.0);
551 Y1->Update(1.0, *X1, 1.0);
552 Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0);
554 L_->Multiply(Trans, Y1temp, *Y1);
555 Y1->Update(1.0, Y1temp, 1.0);
559 if (
Graph_.LevelOverlap()>0 &&
Graph_.L_Graph().DomainMap().DistributedGlobal())
589 int LevelFill = A.
Graph().LevelFill();
590 int LevelOverlap = A.
Graph().LevelOverlap();
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 = ";
#define EPETRA_CHK_ERR(a)
const double Epetra_MinDouble
#define IFPACK_CHK_ERR(ifpack_err)
std::ostream & operator<<(std::ostream &os, const Ifpack_CrsRick &A)
<< operator will work for Ifpack_CrsRick.
adjacency_list< vecS, vecS, undirectedS, no_property, property< edge_weight_t, double > > Graph
Epetra_Flops * GetFlopCounter() const
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
void UpdateFlops(int Flops_in) const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int FillComplete(bool OptimizeDataStorage=true)
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int MaxNumEntries() const
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int ExtractView(double **V) const
Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a giv...
int NumMyCols() const
Returns the number of local matrix columns.
const Epetra_CrsMatrix & A_
const Ifpack_IlukGraph & Graph_
Ifpack_CrsRick(const Epetra_CrsMatrix &A, const Ifpack_IlukGraph &Graph)
Ifpack_CrsRick constuctor with variable number of indices per row.
virtual int NumGlobalDiagonals() const
Returns the number of diagonal entries found in the global input graph.
Epetra_MultiVector * OverlapX_
Epetra_CombineMode OverlapMode_
virtual int NumMyDiagonals() const
Returns the number of diagonal entries found in the local input graph.
int SetParameters(const Teuchos::ParameterList ¶meterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
int Factor()
Compute IC factor L using the specified graph, diagonal perturbation thresholds and relaxation parame...
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
Epetra_MultiVector * OverlapY_
int NumMyRows() const
Returns the number of local matrix rows.
int InitValues()
Initialize L and U with values from user matrix A.
int NumGlobalRows() const
Returns the number of global matrix rows.
const Ifpack_IlukGraph & Graph() const
Returns the address of the Ifpack_IlukGraph associated with this factored matrix.
virtual ~Ifpack_CrsRick()
Ifpack_CrsRick Destructor.
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
void SetFactored(bool Flag)
const Epetra_CrsMatrix & U() const
Returns the address of the U factor associated with this factored matrix.
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.
int SetAllocated(bool Flag)
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local IC set of factors.
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 SetValuesInitialized(bool Flag)
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
void set_parameters(const Teuchos::ParameterList ¶meterlist, param_struct ¶ms, bool cerr_warning_if_unused)
double double_params[FIRST_INT_PARAM]