44 #include "Ifpack_HIPS.h"
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;
71 IsInitialized_(false),
79 void Ifpack_HIPS::Destroy(){
87 int Ifpack_HIPS::Initialize(){
88 if(Comm().NumProc() != 1) IsParallel_ =
true;
89 else IsParallel_ =
false;
94 int Ifpack_HIPS::SetParameters(Teuchos::ParameterList& parameterlist){
98 HIPS_id=List_.get(
"hips: id",-1);
99 if(HIPS_id==-1) IFPACK_CHK_ERR(-1);
102 HIPS_SetDefaultOptions(HIPS_id,List_.get(
"hips: strategy",HIPS_ITERATIVE));
106 if(!MpiComm) IFPACK_CHK_ERR(-2);
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(){
142 if(HIPS_id==-1) IFPACK_CHK_ERR(-1);
143 int N=A_->NumMyRows(), nnz=A_->NumMyNonzeros();
145 int mypid=Comm.
MyPID();
148 int *rowptr,*colind,ierr,maxnr,Nr;
151 const Epetra_Map &RowMap=A_->RowMatrixRowMap();
152 const Epetra_Map &ColMap=A_->RowMatrixColMap();
155 maxnr=A_->MaxNumEntries();
156 colind=
new int[maxnr];
157 values=
new double[maxnr];
165 ColGIDs.Import(RowGIDs,RowToCol,Insert);
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);
179 if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-2);
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]));
186 if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-3);
189 ierr=HIPS_GraphEnd(HIPS_id);
191 if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-4);
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();}
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);
237 if(ierr!=HIPS_SUCCESS){HIPS_PrintError(ierr);IFPACK_CHK_ERR(-5);}
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]);
244 if(ierr!=HIPS_SUCCESS){HIPS_PrintError(ierr);IFPACK_CHK_ERR(-6);}
247 ierr = HIPS_AssemblyEnd(HIPS_id);
249 if(ierr!=HIPS_SUCCESS){HIPS_PrintError(ierr);IFPACK_CHK_ERR(-7);}
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);
259 if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-11);
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));
300 if (X.NumVectors() != Y.NumVectors())
304 if(X.NumVectors()!=1) IFPACK_CHK_ERR(-42);
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);
317 if(rv!=HIPS_SUCCESS) IFPACK_CHK_ERR(-11);
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;
338 double Ifpack_HIPS::Condest(
const Ifpack_CondestType CT,
int MyGlobalElements(int *MyGlobalElementList) const
virtual int MyPID() const =0
int NumMyElements() const
virtual int NumProc() const =0
virtual int ScanSum(double *MyVals, double *ScanSums, int Count) const =0
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const