45 #if defined(HAVE_IFPACK_HIPS) && defined(HAVE_MPI) 53 #include "Epetra_MpiComm.h" 54 #include "Epetra_IntVector.h" 55 #include "Epetra_Import.h" 56 #include "Teuchos_ParameterList.hpp" 57 #include "Teuchos_RefCountPtr.hpp" 58 #ifdef IFPACK_NODE_AWARE_CODE 59 extern int ML_NODE_ID;
62 using Teuchos::RefCountPtr;
67 Ifpack_HIPS::Ifpack_HIPS(Epetra_RowMatrix*
A):
71 IsInitialized_(
false),
79 void Ifpack_HIPS::Destroy(){
87 int Ifpack_HIPS::Initialize(){
88 if(Comm().NumProc() != 1) IsParallel_ =
true;
89 else IsParallel_ =
false;
98 HIPS_id=List_.
get(
"hips: id",-1);
102 HIPS_SetDefaultOptions(HIPS_id,List_.get(
"hips: strategy",HIPS_ITERATIVE));
105 const Epetra_MpiComm* MpiComm=
dynamic_cast<const Epetra_MpiComm*
>(&A_->Comm());
107 HIPS_SetCommunicator(HIPS_id,MpiComm->GetMpiComm());
110 HIPS_SetOptionINT(HIPS_id,HIPS_SYMMETRIC,List_.get(
"hips: symmetric",0));
111 HIPS_SetOptionINT(HIPS_id,HIPS_VERBOSE,List_.get(
"hips: setup output",1));
112 HIPS_SetOptionINT(HIPS_id,HIPS_LOCALLY,List_.get(
"hips: fill",0));
116 HIPS_SetOptionINT(HIPS_id,HIPS_REORDER,List_.get(
"hips: reorder",1));
117 HIPS_SetOptionINT(HIPS_id,HIPS_GRAPH_SYM,List_.get(
"hips: graph symmetric",0));
118 HIPS_SetOptionINT(HIPS_id,HIPS_FORTRAN_NUMBERING,List_.get(
"hips: fortran numbering",0));
119 HIPS_SetOptionINT(HIPS_id,HIPS_DOF,List_.get(
"hips: dof per node",1));
123 HIPS_SetOptionINT(HIPS_id,HIPS_ITMAX,-1);
125 HIPS_SetOptionINT(HIPS_id,HIPS_KRYLOV_RESTART,-1);
128 HIPS_SetOptionREAL(HIPS_id,HIPS_PREC,1e-100);
131 if(List_.get(
"hips: strategy",HIPS_ITERATIVE)==HIPS_ITERATIVE){
132 HIPS_SetOptionREAL(HIPS_id, HIPS_DROPTOL0, List_.get(
"hips: drop tolerance",1e-2));
133 HIPS_SetOptionREAL(HIPS_id, HIPS_DROPTOL1, List_.get(
"hips: drop tolerance",1e-2));
134 HIPS_SetOptionREAL(HIPS_id, HIPS_DROPTOLE, List_.get(
"hips: drop tolerance",1e-2));
141 int Ifpack_HIPS::Compute(){
143 int N=A_->NumMyRows(),
nnz=A_->NumMyNonzeros();
144 const Epetra_Comm &Comm=A_->Comm();
145 int mypid=Comm.MyPID();
148 int *rowptr,*colind,ierr,maxnr,Nr;
150 Epetra_CrsMatrix *Acrs=
dynamic_cast<Epetra_CrsMatrix*
>(&*A_);
151 const Epetra_Map &RowMap=A_->RowMatrixRowMap();
152 const Epetra_Map &ColMap=A_->RowMatrixColMap();
153 if(Acrs) Acrs->ExtractCrsDataPointers(rowptr,colind,values);
155 maxnr=A_->MaxNumEntries();
156 colind=
new int[maxnr];
157 values=
new double[maxnr];
161 RowMap0_=
rcp(
new Epetra_Map(-1,
N,0,Comm));
162 Epetra_IntVector RowGIDs(View,RowMap,RowMap0_->MyGlobalElements());
163 Epetra_IntVector ColGIDs(ColMap);
164 Epetra_Import RowToCol(ColMap,RowMap);
165 ColGIDs.Import(RowGIDs,RowToCol,Insert);
166 ColMap0_=
rcp(
new Epetra_Map(-1,ColMap.NumMyElements(),ColGIDs.Values(),0,Comm));
171 gcolind=
new int[
nnz];
172 for(
int j=0;j<
nnz;j++) gcolind[j]=RowMap0_->GID(colind[j]);
173 ierr = HIPS_GraphDistrCSR(HIPS_id,A_->NumGlobalRows(),A_->NumMyRows(),RowMap0_->MyGlobalElements(),
178 ierr=HIPS_GraphBegin(HIPS_id,A_->NumGlobalRows(),
nnz);
182 for(
int i=0;i<
N;i++){
183 A_->ExtractMyRowCopy(i,maxnr,Nr,values,colind);
184 for(
int j=0;j<Nr;j++){
185 ierr=HIPS_GraphEdge(HIPS_id,RowMap0_->GID(i),ColMap0_->GID(colind[j]));
189 ierr=HIPS_GraphEnd(HIPS_id);
197 Epetra_Map OnePerProcMap(-1,1,0,Comm);
198 Epetra_IntVector RowsPerProc(OnePerProcMap);
199 Epetra_IntVector RowGID(View,*RowMap0_,RowMap0_->MyGlobalElements());
202 Comm.ScanSum(&
N,&(RowsPerProc[0]),1);
205 int OPP_els=0,RPP_els=0;
206 if(!mypid){OPP_els=Comm.NumProc(); RPP_els=A_->NumGlobalRows();}
207 Epetra_Map OPPMap_0(-1,OPP_els,0,Comm);
208 Epetra_Map RPPMap_0(-1,RPP_els,0,Comm);
209 Epetra_Import OPP_importer(OPPMap_0,OnePerProcMap);
210 Epetra_Import RPP_importer(RPPMap_0,*RowMap0_);
213 Epetra_IntVector OPP_0(OPPMap_0);
214 Epetra_IntVector RPP_0(RPPMap_0);
215 OPP_0.Import(RowsPerProc,OPP_importer,Add);
216 RPP_0.Import(RowGID,RPP_importer,Add);
221 mapptr=
new int[Comm.NumProc()+1];
223 for(
int i=0;i<Comm.NumProc();i++)
224 mapptr[i+1]=OPP_0[i];
227 ierr=HIPS_SetPartition(HIPS_id,A_->Comm().NumProc(),mapptr,RPP_0.Values());
228 HIPS_ExitOnError(ierr);
233 ierr = HIPS_MatrixDistrCSR(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),rowptr,gcolind,values,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_FOOL,0);
236 ierr = HIPS_AssemblyBegin(HIPS_id, nnz, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_FOOL,0);
240 for(
int i=0;i<
N;i++){
241 A_->ExtractMyRowCopy(i,maxnr,Nr,values,colind);
242 for(
int j=0;j<Nr;j++){
243 ierr = HIPS_AssemblySetValue(HIPS_id,RowMap0_->GID(i),ColMap0_->GID(colind[j]), values[j]);
247 ierr = HIPS_AssemblyEnd(HIPS_id);
253 double *X=
new double[A_->NumMyRows()];
254 double *Y=
new double[A_->NumMyRows()];
255 for(
int i=0;i<A_->NumMyRows();i++) X[i]=1.0;
258 ierr=HIPS_SetRHS(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),X,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_FOOL);
261 ierr=HIPS_GetSolution(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),Y,HIPS_ASSEMBLY_FOOL);
262 if(ierr!=HIPS_SUCCESS) {
263 HIPS_PrintError(ierr);
268 HIPS_SetOptionINT(HIPS_id,HIPS_VERBOSE,List_.get(
"hips: iteration output",0));
272 HIPS_GetInfoINT(HIPS_id,HIPS_INFO_NNZ,&nnzP);
273 if(nnzP>0) sprintf(Label_,
"Ifpack_HIPS [dt=%4.1e fill=%4.2f]",List_.get(
"hips: drop tolerance",1e-2),(double)nnzP/(
double)A_->NumGlobalNonzeros());
274 else sprintf(Label_,
"Ifpack_HIPS [dt=%4.1e]",List_.get(
"hips: drop tolerance",1e-2));
295 int Ifpack_HIPS::ApplyInverse(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
300 if (X.NumVectors() != Y.NumVectors())
307 Teuchos::RefCountPtr< Epetra_MultiVector > X2;
308 if (X.Pointers()[0] == Y.Pointers()[0])
313 Time_.ResetStartTime();
316 rv=HIPS_SetRHS(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),(*X2)[0],HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_FOOL);
319 rv=HIPS_GetSolution(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),Y[0],HIPS_ASSEMBLY_FOOL);
320 if(rv!=HIPS_SUCCESS) {
326 ApplyInverseTime_ += Time_.ElapsedTime();
332 std::ostream& Ifpack_HIPS::Print(std::ostream& os)
const{
333 os<<
"Need to add meaningful output"<<endl;
341 Epetra_RowMatrix* Matrix_in){
UF_long CHOLMOD() nnz(cholmod_sparse *A, cholmod_common *Common)
T & get(ParameterList &l, const std::string &name)
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define IFPACK_CHK_ERR(ifpack_err)