43 #include "Epetra_ConfigDefs.h" 44 #include "Epetra_Comm.h" 45 #include "Epetra_Map.h" 46 #include "Epetra_CrsGraph.h" 47 #include "Epetra_VbrMatrix.h" 48 #include "Epetra_RowMatrix.h" 49 #include "Epetra_MultiVector.h" 51 #include <Teuchos_ParameterList.hpp> 56 : UserMatrixIsVbr_(
false),
57 UserMatrixIsCrs_(
false),
59 Comm_(Graph_in.Comm()),
63 ValuesInitialized_(
false),
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)
188 Teuchos::RefCountPtr<Epetra_CrsMatrix> OverlapA =
Teuchos::rcp( (Epetra_CrsMatrix *) &
A,
false );
194 EPETRA_CHK_ERR(OverlapA->FillComplete());
198 int MaxNumEntries = OverlapA->MaxNumEntries();
211 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME LONG LONG 228 Teuchos::RefCountPtr<Epetra_VbrMatrix> OverlapA =
Teuchos::rcp( (Epetra_VbrMatrix *) &
A,
false );
234 EPETRA_CHK_ERR(OverlapA->FillComplete());
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());
276 EPETRA_CHK_ERR(
D_->ExtractView(&DV));
283 EPETRA_CHK_ERR(OverlapA.ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]));
291 for (j=0; j< NumIn; j++) {
299 else if (k < 0) {EPETRA_CHK_ERR(-1);}
315 if (DiagFound) NumNonzeroDiags++;
320 EPETRA_CHK_ERR(
L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
323 EPETRA_CHK_ERR(
L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));
329 EPETRA_CHK_ERR(
U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
332 EPETRA_CHK_ERR(
U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));
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;
372 double MinDiagonalValue = Epetra_MinDouble;
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;
408 EPETRA_CHK_ERR(
L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));
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];
431 EPETRA_CHK_ERR(
U_->ExtractMyRowView(j, NumUU, UUV, UUI));
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 456 EPETRA_CHK_ERR(
L_->ReplaceMyValues(i, NumL, LV, LI));
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];
476 EPETRA_CHK_ERR(
U_->ReplaceMyValues(i, NumU, UV, UI));
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;
496 EPETRA_CHK_ERR(
Graph_.
L_Graph().RowMap().Comm().SumAll(¤t_flops, &total_flops, 1));
499 total_flops += (double)
L_->NumGlobalNonzeros();
500 total_flops += (double)
D_->GlobalLength();
501 if (
RelaxValue_!=0.0) total_flops += 2 * (double)
D_->GlobalLength();
503 UpdateFlops(total_flops);
514 Epetra_MultiVector& Y)
const {
520 Teuchos::RefCountPtr<Epetra_MultiVector> X1;
521 Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
522 EPETRA_CHK_ERR(
GenerateXY(Trans, X, Y, &X1, &Y1));
526 bool UnitDiagonal =
true;
528 #ifdef IFPACK_FLOPCOUNTERS 529 Epetra_Flops * counter = this->GetFlopCounter();
531 L_->SetFlopCounter(*counter);
532 Y1->SetFlopCounter(*counter);
533 U_->SetFlopCounter(*counter);
539 EPETRA_CHK_ERR(
L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1));
540 EPETRA_CHK_ERR(Y1->Multiply(1.0, *
D_, *Y1, 0.0));
541 EPETRA_CHK_ERR(
U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1));
545 EPETRA_CHK_ERR(
U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1));
546 EPETRA_CHK_ERR(Y1->Multiply(1.0, *
D_, *Y1, 0.0));
547 EPETRA_CHK_ERR(
L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1));
555 Epetra_MultiVector& Y)
const {
561 Teuchos::RefCountPtr<Epetra_MultiVector> X1;
562 Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
563 EPETRA_CHK_ERR(
GenerateXY(Trans, X, Y, &X1, &Y1));
565 #ifdef IFPACK_FLOPCOUNTERS 566 Epetra_Flops * counter = this->GetFlopCounter();
568 L_->SetFlopCounter(*counter);
569 Y1->SetFlopCounter(*counter);
570 U_->SetFlopCounter(*counter);
575 EPETRA_CHK_ERR(
U_->Multiply(Trans, *X1, *Y1));
576 EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0));
577 EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0));
578 Epetra_MultiVector Y1temp(*Y1);
579 EPETRA_CHK_ERR(
L_->Multiply(Trans, Y1temp, *Y1));
580 EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0));
585 EPETRA_CHK_ERR(
L_->Multiply(Trans, *X1, *Y1));
586 EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0));
587 EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0));
588 Epetra_MultiVector Y1temp(*Y1);
589 EPETRA_CHK_ERR(
U_->Multiply(Trans, Y1temp, *Y1));
590 EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0));
603 Epetra_Vector Ones(
U_->DomainMap());
604 Epetra_Vector OnesResult(
L_->RangeMap());
607 EPETRA_CHK_ERR(
Solve(Trans, Ones, OnesResult));
608 EPETRA_CHK_ERR(OnesResult.Abs(OnesResult));
609 EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate));
616 if (!BG.IndicesAreLocal()) {EPETRA_CHK_ERR(-1);}
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++;
662 int jstart = EPETRA_MAX(0,i-RowDim+1);
664 for (
int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
667 EPETRA_CHK_ERR(PG.InsertMyIndices(i, NumEntries, &tmpIndices[0]));
681 int MaxElementSize = BlockMap.MaxElementSize();
682 int PtNumMyElements = BlockMap.NumMyPoints();
684 std::vector<int> PtMyGlobalElements_int;
685 std::vector<long long> PtMyGlobalElements_LL;
687 int NumMyElements = BlockMap.NumMyElements();
690 if (PtNumMyElements>0) {
691 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 692 if(BlockMap.GlobalIndicesInt()) {
693 PtMyGlobalElements_int.resize(PtNumMyElements);
694 for (
int i=0; i<NumMyElements; i++) {
695 int StartID = BlockMap.GID(i)*MaxElementSize;
696 int ElementSize = BlockMap.ElementSize(i);
697 for (
int j=0; j<ElementSize; j++) PtMyGlobalElements_int[curID++] = StartID+j;
702 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 703 if(BlockMap.GlobalIndicesLongLong()) {
704 PtMyGlobalElements_LL.resize(PtNumMyElements);
705 for (
int i=0; i<NumMyElements; i++) {
706 long long StartID = BlockMap.GID64(i)*MaxElementSize;
707 int ElementSize = BlockMap.ElementSize(i);
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 719 if(BlockMap.GlobalIndicesInt())
720 (*PointMap) =
Teuchos::rcp(
new Epetra_Map(-1, PtNumMyElements, &PtMyGlobalElements_int[0], BlockMap.IndexBase(), BlockMap.Comm()) );
723 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 724 if(BlockMap.GlobalIndicesLongLong())
725 (*PointMap) =
Teuchos::rcp(
new Epetra_Map(-1LL, PtNumMyElements, &PtMyGlobalElements_LL[0], BlockMap.IndexBase64(), BlockMap.Comm()) );
728 throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
730 if (!BlockMap.PointSameAs(*(*PointMap))) {EPETRA_CHK_ERR(-1);}
735 const Epetra_MultiVector& Xin,
const Epetra_MultiVector& Yin,
736 Teuchos::RefCountPtr<Epetra_MultiVector>* Xout,
737 Teuchos::RefCountPtr<Epetra_MultiVector>* Yout)
const {
741 if (Xin.NumVectors()!=Yin.NumVectors()) EPETRA_CHK_ERR(-1);
744 (*Xout) =
Teuchos::rcp( (Epetra_MultiVector *) &Xin,
false );
745 (*Yout) =
Teuchos::rcp( (Epetra_MultiVector *) &Yin,
false );
749 if (
VbrX_!=Teuchos::null) {
750 if (
VbrX_->NumVectors()!=Xin.NumVectors()) {
751 VbrX_ = Teuchos::null;
752 VbrY_ = Teuchos::null;
755 if (
VbrX_==Teuchos::null) {
760 EPETRA_CHK_ERR(
VbrX_->ResetView((*Xout)->Pointers()));
761 EPETRA_CHK_ERR(
VbrY_->ResetView((*Yout)->Pointers()));
770 if (
OverlapX_->NumVectors()!=Xin.NumVectors()) {
780 EPETRA_CHK_ERR(
OverlapX_->Import(*(*Xout),*
U_->Importer(), Insert));
783 EPETRA_CHK_ERR(
OverlapX_->Import(*(*Xout),*
L_->Exporter(), Insert));
802 int LevelFill =
A.Graph().LevelFill();
803 int LevelOverlap =
A.Graph().LevelOverlap();
804 Epetra_CrsMatrix & L = (Epetra_CrsMatrix &)
A.L();
805 Epetra_CrsMatrix & U = (Epetra_CrsMatrix &)
A.U();
806 Epetra_Vector &
D = (Epetra_Vector &)
A.D();
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 = ";
Teuchos::RefCountPtr< const Epetra_Map > L_RangeMap_
int NumMyRows() const
Returns the number of local matrix rows.
int BlockMap2PointMap(const Epetra_BlockMap &BlockMap, Teuchos::RefCountPtr< Epetra_Map > *PointMap)
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
void SetFactored(bool Flag)
double double_params[FIRST_INT_PARAM]
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
Teuchos::RefCountPtr< Epetra_CrsGraph > U_Graph_
Teuchos::RefCountPtr< const Epetra_Map > U_DomainMap_
const Epetra_CrsMatrix & L() const
Returns the address of the L factor associated with this factored matrix.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapY_
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
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 ...
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...
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
int BlockGraph2PointGraph(const Epetra_CrsGraph &BG, Epetra_CrsGraph &PG, bool Upper)
int SetAllocated(bool Flag)
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
void SetValuesInitialized(bool Flag)
std::ostream & operator<<(std::ostream &os, const Ifpack_CrsRiluk &A)
<< operator will work for Ifpack_CrsRiluk.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
Teuchos::RefCountPtr< Epetra_CrsGraph > L_Graph_
Teuchos::RefCountPtr< Epetra_Map > IlukDomainMap_
Epetra_CombineMode overlap_mode
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
int InitAllValues(const Epetra_RowMatrix &A, int MaxNumEntries)
int SetParameters(const Teuchos::ParameterList ¶meterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
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...
const Ifpack_IlukGraph & Graph_
int NumMyCols() const
Returns the number of local matrix columns.
Teuchos::RefCountPtr< Epetra_Map > IlukRowMap_
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapX_
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
Teuchos::RefCountPtr< Epetra_MultiVector > VbrY_
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.
int GenerateXY(bool Trans, const Epetra_MultiVector &Xin, const Epetra_MultiVector &Yin, Teuchos::RefCountPtr< Epetra_MultiVector > *Xout, Teuchos::RefCountPtr< Epetra_MultiVector > *Yout) const
void set_parameters(const Teuchos::ParameterList ¶meterlist, param_struct ¶ms, bool cerr_warning_if_unused)
Teuchos::RefCountPtr< Epetra_Vector > D_
Teuchos::RefCountPtr< Epetra_MultiVector > VbrX_
virtual ~Ifpack_CrsRiluk()
Ifpack_CrsRiluk Destructor.
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
Epetra_CombineMode OverlapMode_
Teuchos::RefCountPtr< Epetra_Map > IlukRangeMap_
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
Ifpack_CrsRiluk(const Ifpack_IlukGraph &Graph_in)
Ifpack_CrsRiluk constuctor with variable number of indices per row.
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.