51#include <Teuchos_ParameterList.hpp>
56 : UserMatrixIsVbr_(
false),
57 UserMatrixIsCrs_(
false),
59 Comm_(Graph_in.Comm()),
63 ValuesInitialized_(
false),
72 IsOverlapped_ = (Graph_in.LevelOverlap()>0 && Graph_in.DomainMap().DistributedGlobal());
77 : UserMatrixIsVbr_(FactoredMatrix.UserMatrixIsVbr_),
78 UserMatrixIsCrs_(FactoredMatrix.UserMatrixIsCrs_),
79 IsOverlapped_(FactoredMatrix.IsOverlapped_),
80 Graph_(FactoredMatrix.Graph_),
81 IlukRowMap_(FactoredMatrix.IlukRowMap_),
82 IlukDomainMap_(FactoredMatrix.IlukDomainMap_),
83 IlukRangeMap_(FactoredMatrix.IlukRangeMap_),
84 Comm_(FactoredMatrix.Comm_),
85 UseTranspose_(FactoredMatrix.UseTranspose_),
86 NumMyDiagonals_(FactoredMatrix.NumMyDiagonals_),
87 Allocated_(FactoredMatrix.Allocated_),
88 ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
89 Factored_(FactoredMatrix.Factored_),
90 RelaxValue_(FactoredMatrix.RelaxValue_),
91 Athresh_(FactoredMatrix.Athresh_),
92 Rthresh_(FactoredMatrix.Rthresh_),
93 Condest_(FactoredMatrix.Condest_),
94 OverlapMode_(FactoredMatrix.OverlapMode_)
136 if (
Graph().LevelFill()) {
163 bool cerr_warning_if_unused)
167 params.double_params[Ifpack::absolute_threshold] =
Athresh_;
168 params.double_params[Ifpack::relative_threshold] =
Rthresh_;
173 RelaxValue_ = params.double_params[Ifpack::relax_value];
174 Athresh_ = params.double_params[Ifpack::absolute_threshold];
175 Rthresh_ = params.double_params[Ifpack::relative_threshold];
188 Teuchos::RefCountPtr<Epetra_CrsMatrix> OverlapA = Teuchos::rcp( (
Epetra_CrsMatrix *) &A,
false );
198 int MaxNumEntries = OverlapA->MaxNumEntries();
211#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
228 Teuchos::RefCountPtr<Epetra_VbrMatrix> OverlapA = Teuchos::rcp( (
Epetra_VbrMatrix *) &A,
false );
240 int MaxNumEntries = OverlapA->MaxNumNonzeros();
255 int NumIn, NumL, NumU;
257 int NumNonzeroDiags = 0;
260 std::vector<int> InI(MaxNumEntries);
261 std::vector<int> LI(MaxNumEntries);
262 std::vector<int> UI(MaxNumEntries);
263 std::vector<double> InV(MaxNumEntries);
264 std::vector<double> LV(MaxNumEntries);
265 std::vector<double> UV(MaxNumEntries);
267 bool ReplaceValues = (
L_->StaticGraph() ||
L_->IndicesAreLocal());
283 EPETRA_CHK_ERR(OverlapA.ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]));
291 for (j=0; j< NumIn; j++) {
315 if (DiagFound) NumNonzeroDiags++;
338 if (!ReplaceValues) {
352 int TotalNonzeroDiags = 0;
353 EPETRA_CHK_ERR(
Graph_.L_Graph().RowMap().Comm().SumAll(&NumNonzeroDiags, &TotalNonzeroDiags, 1));
355 if (NumNonzeroDiags !=
NumMyRows()) ierr = 1;
373 double MaxDiagonalValue = 1.0/MinDiagonalValue;
377 int * LI=0, * UI = 0;
378 double * LV=0, * UV = 0;
379 int NumIn, NumL, NumU;
382 int MaxNumEntries =
L_->MaxNumEntries() +
U_->MaxNumEntries() + 1;
384 std::vector<int> InI(MaxNumEntries);
385 std::vector<double> InV(MaxNumEntries);
389 ierr =
D_->ExtractView(&DV);
391#ifdef IFPACK_FLOPCOUNTERS
392 int current_madds = 0;
401 for (j=0; j<
NumMyCols(); j++) colflag[j] = - 1;
407 NumIn = MaxNumEntries;
415 EPETRA_CHK_ERR(
U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
421 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
423 double diagmod = 0.0;
425 for (
int jj=0; jj<NumL; jj++) {
427 double multiplier = InV[jj];
434 for (k=0; k<NumUU; k++) {
435 int kk = colflag[UUI[k]];
437 InV[kk] -= multiplier*UUV[k];
438#ifdef IFPACK_FLOPCOUNTERS
445 for (k=0; k<NumUU; k++) {
446 int kk = colflag[UUI[k]];
447 if (kk>-1) InV[kk] -= multiplier*UUV[k];
448 else diagmod -= multiplier*UUV[k];
449#ifdef IFPACK_FLOPCOUNTERS
466 if (fabs(DV[i]) > MaxDiagonalValue) {
467 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
468 else DV[i] = MinDiagonalValue;
473 for (j=0; j<NumU; j++) UV[j] *= DV[i];
480 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
485 if( !
L_->LowerTriangular() )
487 if( !
U_->UpperTriangular() )
490#ifdef IFPACK_FLOPCOUNTERS
493 double current_flops = 2 * current_madds;
494 double total_flops = 0;
499 total_flops += (double)
L_->NumGlobalNonzeros();
500 total_flops += (double)
D_->GlobalLength();
501 if (
RelaxValue_!=0.0) total_flops += 2 * (double)
D_->GlobalLength();
520 Teuchos::RefCountPtr<Epetra_MultiVector> X1;
521 Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
526 bool UnitDiagonal =
true;
528#ifdef IFPACK_FLOPCOUNTERS
531 L_->SetFlopCounter(*counter);
532 Y1->SetFlopCounter(*counter);
533 U_->SetFlopCounter(*counter);
561 Teuchos::RefCountPtr<Epetra_MultiVector> X1;
562 Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
565#ifdef IFPACK_FLOPCOUNTERS
568 L_->SetFlopCounter(*counter);
569 Y1->SetFlopCounter(*counter);
570 U_->SetFlopCounter(*counter);
618 int * ColFirstPointInElementList = BG.RowMap().FirstPointInElementList();
619 int * ColElementSizeList = BG.RowMap().ElementSizeList();
620 if (BG.Importer()!=0) {
621 ColFirstPointInElementList = BG.ImportMap().FirstPointInElementList();
622 ColElementSizeList = BG.ImportMap().ElementSizeList();
625 int Length = (BG.MaxNumIndices()+1) * BG.ImportMap().MaxMyElementSize();
626 std::vector<int> tmpIndices(Length);
628 int BlockRow, BlockOffset, NumEntries;
632 int NumMyRows_tmp = PG.NumMyRows();
634 for (
int i=0; i<NumMyRows_tmp; i++) {
635 EPETRA_CHK_ERR(BG.RowMap().FindLocalElementID(i, BlockRow, BlockOffset));
636 EPETRA_CHK_ERR(BG.ExtractMyRowView(BlockRow, NumBlockEntries, BlockIndices));
638 int * ptr = &tmpIndices[0];
640 int RowDim = BG.RowMap().ElementSize(BlockRow);
647 int jstop =
EPETRA_MIN(NumMyRows_tmp,i+RowDim-BlockOffset);
648 for (
int j= jstart; j< jstop; j++) {*ptr++ = j; NumEntries++;}
651 for (
int j=0; j<NumBlockEntries; j++) {
652 int ColDim = ColElementSizeList[BlockIndices[j]];
653 NumEntries += ColDim;
654 assert(NumEntries<=Length);
655 int Index = ColFirstPointInElementList[BlockIndices[j]];
656 for (
int k=0; k < ColDim; k++) *ptr++ = Index++;
664 for (
int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
667 EPETRA_CHK_ERR(PG.InsertMyIndices(i, NumEntries, &tmpIndices[0]));
684 std::vector<int> PtMyGlobalElements_int;
685 std::vector<long long> PtMyGlobalElements_LL;
690 if (PtNumMyElements>0) {
691#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
693 PtMyGlobalElements_int.resize(PtNumMyElements);
694 for (
int i=0; i<NumMyElements; i++) {
695 int StartID = BlockMap.
GID(i)*MaxElementSize;
697 for (
int j=0; j<ElementSize; j++) PtMyGlobalElements_int[curID++] = StartID+j;
702#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
704 PtMyGlobalElements_LL.resize(PtNumMyElements);
705 for (
int i=0; i<NumMyElements; i++) {
706 long long StartID = BlockMap.
GID64(i)*MaxElementSize;
708 for (
int j=0; j<ElementSize; j++) PtMyGlobalElements_LL[curID++] = StartID+j;
713 throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
716 assert(curID==PtNumMyElements);
718#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
720 (*PointMap) = Teuchos::rcp(
new Epetra_Map(-1, PtNumMyElements, &PtMyGlobalElements_int[0], BlockMap.
IndexBase(), BlockMap.
Comm()) );
723#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
725 (*PointMap) = Teuchos::rcp(
new Epetra_Map(-1LL, PtNumMyElements, &PtMyGlobalElements_LL[0], BlockMap.
IndexBase64(), BlockMap.
Comm()) );
728 throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
736 Teuchos::RefCountPtr<Epetra_MultiVector>* Xout,
737 Teuchos::RefCountPtr<Epetra_MultiVector>* Yout)
const {
749 if (
VbrX_!=Teuchos::null) {
750 if (
VbrX_->NumVectors()!=Xin.NumVectors()) {
751 VbrX_ = Teuchos::null;
752 VbrY_ = Teuchos::null;
755 if (
VbrX_==Teuchos::null) {
770 if (
OverlapX_->NumVectors()!=Xin.NumVectors()) {
802 int LevelFill = A.
Graph().LevelFill();
803 int LevelOverlap = A.
Graph().LevelOverlap();
810 os <<
" Level of Fill = "; os << LevelFill;
813 os <<
" Level of Overlap = "; os << LevelOverlap;
817 os <<
" Lower Triangle = ";
823 os <<
" Inverse of Diagonal = ";
829 os <<
" Upper Triangle = ";
#define EPETRA_CHK_ERR(a)
const double Epetra_MinDouble
std::ostream & operator<<(std::ostream &os, const Ifpack_CrsRiluk &A)
<< operator will work for Ifpack_CrsRiluk.
long long IndexBase64() const
int MaxElementSize() const
bool PointSameAs(const Epetra_BlockMap &Map) const
long long GID64(int LID) const
bool GlobalIndicesInt() const
const Epetra_Comm & Comm() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
Epetra_Flops * GetFlopCounter() const
void UpdateFlops(int Flops_in) const
const Epetra_Map & RangeMap() const
const Epetra_Map & DomainMap() const
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
int BlockGraph2PointGraph(const Epetra_CrsGraph &BG, Epetra_CrsGraph &PG, bool Upper)
int GenerateXY(bool Trans, const Epetra_MultiVector &Xin, const Epetra_MultiVector &Yin, Teuchos::RefCountPtr< Epetra_MultiVector > *Xout, Teuchos::RefCountPtr< Epetra_MultiVector > *Yout) const
const Ifpack_IlukGraph & Graph_
virtual ~Ifpack_CrsRiluk()
Ifpack_CrsRiluk Destructor.
Teuchos::RefCountPtr< Epetra_Map > IlukRowMap_
Teuchos::RefCountPtr< Epetra_MultiVector > VbrX_
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of multiplying U, D and L in that order on an Epetra_MultiVector X in Y.
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
int SetParameters(const Teuchos::ParameterList ¶meterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
int NumMyRows() const
Returns the number of local matrix rows.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapY_
const Epetra_CrsMatrix & U() const
Returns the address of the L factor associated with this factored matrix.
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
Teuchos::RefCountPtr< const Epetra_Map > L_RangeMap_
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapX_
int InitAllValues(const Epetra_RowMatrix &A, int MaxNumEntries)
Teuchos::RefCountPtr< Epetra_Map > IlukDomainMap_
Teuchos::RefCountPtr< const Epetra_Map > U_DomainMap_
Ifpack_CrsRiluk(const Ifpack_IlukGraph &Graph_in)
Ifpack_CrsRiluk constuctor with variable number of indices per row.
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
void SetFactored(bool Flag)
Teuchos::RefCountPtr< Epetra_Vector > D_
void SetValuesInitialized(bool Flag)
Teuchos::RefCountPtr< Epetra_CrsGraph > L_Graph_
Teuchos::RefCountPtr< Epetra_CrsGraph > U_Graph_
Epetra_CombineMode OverlapMode_
Teuchos::RefCountPtr< Epetra_MultiVector > VbrY_
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
int BlockMap2PointMap(const Epetra_BlockMap &BlockMap, Teuchos::RefCountPtr< Epetra_Map > *PointMap)
Teuchos::RefCountPtr< Epetra_Map > IlukRangeMap_
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
const Epetra_CrsMatrix & L() const
Returns the address of the L factor associated with this factored matrix.
int NumMyCols() const
Returns the number of local matrix columns.
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors.
int Solve(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_CrsRiluk forward/back solve on a Epetra_MultiVector X in Y (works for ...
int SetAllocated(bool Flag)
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix.
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]