43 #include "Ifpack_ConfigDefs.h" 44 #include "Ifpack_SILU.h" 45 #ifdef HAVE_IFPACK_SUPERLU 47 #include "Ifpack_CondestType.h" 48 #include "Epetra_ConfigDefs.h" 49 #include "Epetra_Comm.h" 50 #include "Epetra_Map.h" 51 #include "Epetra_RowMatrix.h" 52 #include "Epetra_Vector.h" 53 #include "Epetra_MultiVector.h" 54 #include "Epetra_CrsGraph.h" 55 #include "Epetra_CrsMatrix.h" 56 #include "Teuchos_ParameterList.hpp" 57 #include "Teuchos_RefCountPtr.hpp" 60 extern "C" {
int dfill_diag(
int n, NCformat *Astore);}
62 using Teuchos::RefCountPtr;
65 #ifdef IFPACK_TEUCHOS_TIME_MONITOR 66 # include "Teuchos_TimeMonitor.hpp" 71 A_(rcp(Matrix_in,false)),
72 Comm_(Matrix_in->Comm()),
80 IsInitialized_(false),
87 ApplyInverseTime_(0.0),
93 Teuchos::ParameterList List;
101 void Ifpack_SILU::Destroy()
107 Destroy_CompCol_Permuted(&SAc_);
108 Destroy_SuperNode_Matrix(&SL_);
109 Destroy_CompCol_Matrix(&SU_);
112 Destroy_SuperMatrix_Store(&SA_);
113 if(SY_.Store) Destroy_SuperMatrix_Store(&SY_);
117 delete [] etree_;etree_=0;
118 delete [] perm_r_;perm_r_=0;
119 delete [] perm_c_;perm_c_=0;
126 DropTol_=List.get(
"fact: drop tolerance",DropTol_);
127 FillTol_=List.get(
"fact: zero pivot threshold",FillTol_);
128 FillFactor_=List.get(
"fact: maximum fill factor",FillFactor_);
129 DropRule_=List.get(
"fact: silu drop rule",DropRule_);
132 sprintf(Label_,
"IFPACK SILU (drop=%d, zpv=%f, ffact=%f, rthr=%f)",
133 DropRule(),FillTol(),FillFactor(),DropTol());
138 template<
typename int_type>
139 int Ifpack_SILU::TInitialize()
142 #ifdef IFPACK_TEUCHOS_TIME_MONITOR 143 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_SILU::Initialize");
146 Time_.ResetStartTime();
151 IsInitialized_ =
false;
153 Epetra_CrsMatrix* CrsMatrix;
154 CrsMatrix =
dynamic_cast<Epetra_CrsMatrix*
>(&*A_);
156 if(CrsMatrix && CrsMatrix->RowMap().SameAs(CrsMatrix->ColMap()) && CrsMatrix->IndicesAreContiguous()){
158 Aover_=rcp(CrsMatrix,
false);
160 else if(CrsMatrix && CrsMatrix->IndicesAreContiguous()){
162 int size = A_->MaxNumEntries();
163 int N=A_->NumMyRows();
164 Aover_ = rcp(
new Epetra_CrsMatrix(Copy,A_->RowMatrixRowMap(), size));
165 std::vector<int_type> Indices(size);
166 std::vector<double> Values(size);
168 int i,j,ct,*rowptr,*colind;
170 IFPACK_CHK_ERR(CrsMatrix->ExtractCrsDataPointers(rowptr,colind,values));
174 for(j=rowptr[i],ct=0;j<rowptr[i+1];j++){
176 Indices[ct]= (int_type) CrsMatrix->GCID64(colind[j]);
177 Values[ct]=values[j];
181 Aover_->InsertGlobalValues((int_type) CrsMatrix->GRID64(i),ct,&Values[0],&Indices[0]);
183 IFPACK_CHK_ERR(Aover_->FillComplete(CrsMatrix->RowMap(),CrsMatrix->RowMap()));
187 int size = A_->MaxNumEntries();
188 Aover_ = rcp(
new Epetra_CrsMatrix(Copy,A_->RowMatrixRowMap(), size));
189 if (Aover_.get() == 0) IFPACK_CHK_ERR(-5);
191 std::vector<int> Indices1(size);
192 std::vector<int_type> Indices2(size);
193 std::vector<double> Values1(size),Values2(size);
197 int N=A_->NumMyRows();
198 for (
int i = 0 ; i < N ; ++i) {
200 int_type GlobalRow = (int_type) A_->RowMatrixRowMap().GID64(i);
201 IFPACK_CHK_ERR(A_->ExtractMyRowCopy(i, size, NumEntries,
202 &Values1[0], &Indices1[0]));
206 for (
int j=0; j < NumEntries ; ++j) {
208 Indices2[ct] = (int_type) A_->RowMatrixColMap().GID64(Indices1[j]);
209 Values2[ct]=Values1[j];
213 IFPACK_CHK_ERR(Aover_->InsertGlobalValues(GlobalRow,ct,
214 &Values2[0],&Indices2[0]));
216 IFPACK_CHK_ERR(Aover_->FillComplete(A_->RowMatrixRowMap(),A_->RowMatrixRowMap()));
220 Aover_->OptimizeStorage();
221 Graph_=rcp(const_cast<Epetra_CrsGraph*>(&Aover_->Graph()),
false);
223 IsInitialized_ =
true;
225 InitializeTime_ += Time_.ElapsedTime();
231 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 232 if(A_->RowMatrixRowMap().GlobalIndicesInt()) {
233 return TInitialize<int>();
237 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 238 if(A_->RowMatrixRowMap().GlobalIndicesLongLong()) {
239 return TInitialize<long long>();
243 throw "Ifpack_SILU::Initialize: GlobalIndices type unknown for A_";
251 #ifdef IFPACK_TEUCHOS_TIME_MONITOR 252 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_SILU::Compute");
258 Time_.ResetStartTime();
263 ilu_set_default_options(&options_);
264 options_.ILU_DropTol=DropTol_;
265 options_.ILU_FillTol=FillTol_;
266 options_.ILU_DropRule=DropRule_;
267 options_.ILU_FillFactor=FillFactor_;
272 IFPACK_CHK_ERR(Aover_->ExtractCrsDataPointers(rowptr,colind,values));
273 int N=Aover_->NumMyRows();
276 dCreate_CompCol_Matrix(&SA_,N,N,Aover_->NumMyNonzeros(),
277 values,colind,rowptr,SLU_NC,SLU_D,SLU_GE);
281 dfill_diag(N, (NCformat*)SA_.Store);
289 int permc_spec=options_.ColPerm;
290 if ( permc_spec != MY_PERMC && options_.Fact == DOFACT )
291 get_perm_c(permc_spec,&SA_,perm_c_);
294 sp_preorder(&options_, &SA_, perm_c_, etree_, &SAc_);
297 int panel_size = sp_ienv(1);
298 int relax = sp_ienv(2);
300 dgsitrf(&options_,&SAc_,relax,panel_size,etree_,NULL,0,perm_c_,perm_r_,&SL_,&SU_,
301 #ifdef HAVE_IFPACK_SUPERLU5_API
305 if(info<0) IFPACK_CHK_ERR(info);
309 ComputeTime_ += Time_.ElapsedTime();
315 int Ifpack_SILU::Solve(
bool Trans,
const Epetra_MultiVector& X,
316 Epetra_MultiVector& Y)
const 319 #ifdef IFPACK_TEUCHOS_TIME_MONITOR 320 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_SILU::ApplyInverse - Solve");
325 int nrhs=X.NumVectors();
326 int N=A_->NumMyRows();
333 if(SY_.Store) Destroy_SuperMatrix_Store(&SY_);
335 dCreate_Dense_Matrix(&SY_,N,nrhs,Y[0],N,SLU_DN,SLU_D,SLU_GE);
338 dgstrs(TRANS,&SL_,&SU_,perm_c_,perm_r_,&SY_,&stat_,&info);
339 if(!info) IFPACK_CHK_ERR(info);
346 int Ifpack_SILU::Multiply(
bool Trans,
const Epetra_MultiVector& X,
347 Epetra_MultiVector& Y)
const 359 Epetra_MultiVector& Y)
const 362 #ifdef IFPACK_TEUCHOS_TIME_MONITOR 363 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_SILU::ApplyInverse");
369 if (X.NumVectors() != Y.NumVectors())
372 Time_.ResetStartTime();
376 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
377 if (X.Pointers()[0] == Y.Pointers()[0])
378 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
380 Xcopy = Teuchos::rcp( &X,
false );
385 ApplyInverseTime_ += Time_.ElapsedTime();
393 const int MaxIters,
const double Tol,
394 Epetra_RowMatrix* Matrix_in)
397 #ifdef IFPACK_TEUCHOS_TIME_MONITOR 398 TEUCHOS_FUNC_TIME_MONITOR(
"Ifpack_SILU::Condest");
404 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
415 if (!
Comm().MyPID()) {
417 os <<
"================================================================================" << endl;
418 os <<
"Ifpack_SILU: " <<
Label() << endl << endl;
419 os <<
"Dropping rule = "<< DropRule() << endl;
420 os <<
"Zero pivot thresh = "<< FillTol() << endl;
421 os <<
"Max fill factor = "<< FillFactor() << endl;
422 os <<
"Drop tolerance = "<< DropTol() << endl;
423 os <<
"Condition number estimate = " <<
Condest() << endl;
424 os <<
"Global number of rows = " << A_->NumGlobalRows64() << endl;
428 if(SL_.Store) fnnz+=((SCformat*)SL_.Store)->nnz;
429 if(SU_.Store) fnnz+=((NCformat*)SU_.Store)->nnz;
430 int lufill=fnnz - A_->NumMyRows();
431 os <<
"No. of nonzeros in L+U = "<< lufill<<endl;
432 os <<
"Fill ratio: nnz(F)/nnz(A) = "<<(double)lufill / (
double)A_->NumMyNonzeros()<<endl;
435 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
436 os <<
"----- ------- -------------- ------------ --------" << endl;
439 <<
" 0.0 0.0" << endl;
440 os <<
"Compute() " << std::setw(5) <<
NumCompute()
446 os <<
" " << std::setw(15) << 0.0 << endl;
453 os <<
" " << std::setw(15) << 0.0 << endl;
454 os <<
"================================================================================" << endl;
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Ifpack_SILU(Epetra_RowMatrix *A)
Constructor.
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
virtual double ComputeTime() const
Returns the time spent in Compute().
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y...
double Condest() const
Returns the computed estimated condition number, or -1.0 if not computed.
const char * Label() const
Returns a character string describing the operator.
virtual int NumCompute() const
Returns the number of calls to Compute().
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
int Compute()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
int Initialize()
Initialize the preconditioner, does not touch matrix values.
int SetParameters(Teuchos::ParameterList ¶meterlist)
Set parameters using a Teuchos::ParameterList object.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
bool UseTranspose() const
Returns the current UseTranspose setting.
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual std::ostream & Print(std::ostream &os) const
Prints on stream basic information about this object.