45 #ifdef HAVE_IFPACK_SUPERLU
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"
72 Comm_(Matrix_in->Comm()),
80 IsInitialized_(
false),
87 ApplyInverseTime_(0.0),
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
146 Time_.ResetStartTime();
151 IsInitialized_ =
false;
158 Aover_=
rcp(CrsMatrix,
false);
162 int size = A_->MaxNumEntries();
163 int N=A_->NumMyRows();
165 std::vector<int_type> Indices(size);
166 std::vector<double> Values(size);
168 int i,j,ct,*rowptr,*colind;
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]);
187 int size = A_->MaxNumEntries();
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);
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];
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();
230 int Ifpack_SILU::Initialize() {
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_";
248 int Ifpack_SILU::Compute()
251 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
255 if (!IsInitialized())
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
309 ComputeTime_ += Time_.ElapsedTime();
319 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
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);
362 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
369 if (X.NumVectors() != Y.NumVectors())
372 Time_.ResetStartTime();
376 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
377 if (X.Pointers()[0] == Y.Pointers()[0])
385 ApplyInverseTime_ += Time_.ElapsedTime();
393 const int MaxIters,
const double Tol,
397 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
411 Ifpack_SILU::Print(std::ostream& os)
const
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;
437 os <<
"Initialize() " << std::setw(5) << NumInitialize()
438 <<
" " << std::setw(15) << InitializeTime()
439 <<
" 0.0 0.0" << endl;
440 os <<
"Compute() " << std::setw(5) << NumCompute()
441 <<
" " << std::setw(15) << ComputeTime()
442 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops();
443 if (ComputeTime() != 0.0)
444 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
446 os <<
" " << std::setw(15) << 0.0 << endl;
447 os <<
"ApplyInverse() " << std::setw(5) << NumApplyInverse()
448 <<
" " << std::setw(15) << ApplyInverseTime()
449 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
450 if (ApplyInverseTime() != 0.0)
451 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
453 os <<
" " << std::setw(15) << 0.0 << endl;
454 os <<
"================================================================================" << endl;
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
bool SameAs(const Epetra_BlockMap &Map) const
T & get(ParameterList &l, const std::string &name)
long long GRID64(int LRID_in) const
const Epetra_Map & ColMap() const
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
const Epetra_Map & RowMap() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool IndicesAreContiguous() const
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
int Solve(int, TYPE *, TYPE *, TYPE *)
#define IFPACK_CHK_ERR(ifpack_err)
long long GCID64(int LCID_in) const
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const