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_)
95int Ifpack_CrsRick::Allocate() {
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;
164 InI =
new int[MaxNumEntries];
165 UI =
new int[MaxNumEntries];
166 InV =
new double[MaxNumEntries];
167 UV =
new double[MaxNumEntries];
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++;
220 SetValuesInitialized(
true);
223 int TotalNonzeroDiags = 0;
240 SetValuesInitialized(
false);
245 double MinDiagonalValue = Epetra_MinDouble;
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);
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];
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];
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;
366 total_flops += (double) L_->NumGlobalNonzeros();
368 if (RelaxValue_!=0.0) total_flops += 2 * (double)D_->
GlobalLength();
392 bool UnitDiagonal =
true;
407#ifdef IFPACK_FLOPCOUNTERS
410 L_->SetFlopCounter(*counter);
418 L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
420 U_->
Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1);
424 U_->
Solve(Upper, Trans, UnitDiagonal, *X1, *Y1);
426 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
448 bool UnitDiagonal =
true;
457 delete OverlapX_; OverlapX_ = 0;
458 delete OverlapY_; OverlapY_ = 0;
470#ifdef IFPACK_FLOPCOUNTERS
473 L_->SetFlopCounter(*counter);
481 L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
483 U_->
Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1);
487 U_->
Solve(Upper, Trans, UnitDiagonal, *X1, *Y1);
489 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
509 bool UnitDiagonal =
true;
518 delete OverlapX_; OverlapX_ = 0;
519 delete OverlapY_; OverlapY_ = 0;
531#ifdef IFPACK_FLOPCOUNTERS
534 L_->SetFlopCounter(*counter);
542 L_->Multiply(Trans, *X1, *Y1);
543 Y1->
Update(1.0, *X1, 1.0);
547 Y1->
Update(1.0, Y1temp, 1.0);
551 Y1->
Update(1.0, *X1, 1.0);
554 L_->Multiply(Trans, Y1temp, *Y1);
555 Y1->
Update(1.0, Y1temp, 1.0);
567 ConditionNumberEstimate = Condest_;
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;
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 = ";
bool DistributedGlobal() const
const Epetra_Comm & Comm() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
Epetra_Flops * GetFlopCounter() const
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
void UpdateFlops(int Flops_in) const
const Epetra_BlockMap & DomainMap() const
const Epetra_BlockMap & RowMap() const
const Epetra_Map & RowMap() 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 Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Abs(const Epetra_MultiVector &A)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int MaxValue(double *Result) const
int PutScalar(double ScalarConstant)
int ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
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.
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.
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.
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.
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 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.
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...
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.