IFPACK Development
Ifpack_HIPS.cpp
00001 #include "Ifpack_HIPS.h"
00002 #if defined(HAVE_IFPACK_HIPS) && defined(HAVE_MPI)
00003 
00004 
00005 #include "Ifpack_Utils.h"
00006 extern "C" {
00007 #include "hips.h"
00008 }
00009 //#include "io.h"
00010 #include "Epetra_MpiComm.h"
00011 #include "Epetra_IntVector.h"
00012 #include "Epetra_Import.h"
00013 #include "Teuchos_ParameterList.hpp"
00014 #include "Teuchos_RefCountPtr.hpp"
00015 #ifdef IFPACK_NODE_AWARE_CODE
00016 extern int ML_NODE_ID;
00017 #endif
00018 
00019 using Teuchos::RefCountPtr;
00020 using Teuchos::rcp;
00021 
00022 
00023 
00024 Ifpack_HIPS::Ifpack_HIPS(Epetra_RowMatrix* A):
00025   A_(rcp(A,false)),
00026   HIPS_id(-1), // Assumes user will initialze HIPS outside 
00027   IsParallel_(false),
00028   IsInitialized_(false),
00029   IsComputed_(false),
00030   Label_(),
00031   Time_(A_->Comm())
00032 {
00033   
00034 }
00035 
00036 void Ifpack_HIPS::Destroy(){
00037   //NTS: Assume user will call HIPS_Finalize elsewhere - HIPS_Clean never needs
00038   //to be called if HIPS_Finalize is called at the end, unless you want to reuse
00039   //a slot.
00040 }
00041 
00042 
00043 
00044 int Ifpack_HIPS::Initialize(){
00045   if(Comm().NumProc() != 1) IsParallel_ = true;
00046   else IsParallel_ = false;
00047   IsInitialized_=true;
00048   return 0;
00049 }
00050 
00051 int Ifpack_HIPS::SetParameters(Teuchos::ParameterList& parameterlist){
00052   List_=parameterlist;
00053 
00054   // Grab the hips ID
00055   HIPS_id=List_.get("hips: id",-1);
00056   if(HIPS_id==-1) IFPACK_CHK_ERR(-1);
00057 
00058   // Set Defaults
00059   HIPS_SetDefaultOptions(HIPS_id,List_.get("hips: strategy",HIPS_ITERATIVE));  
00060   
00061   // Set the communicator
00062   const Epetra_MpiComm* MpiComm=dynamic_cast<const Epetra_MpiComm*>(&A_->Comm());
00063   if(!MpiComm) IFPACK_CHK_ERR(-2);  
00064   HIPS_SetCommunicator(HIPS_id,MpiComm->GetMpiComm());
00065   
00066   // Options
00067   HIPS_SetOptionINT(HIPS_id,HIPS_SYMMETRIC,List_.get("hips: symmetric",0));
00068   HIPS_SetOptionINT(HIPS_id,HIPS_VERBOSE,List_.get("hips: setup output",1));
00069   HIPS_SetOptionINT(HIPS_id,HIPS_LOCALLY,List_.get("hips: fill",0));
00070   // Sadly, this fill option doesn't work for HIPS_ITERATIVE mode, meaning the
00071   // only way to control fill-in is via the drop tolerance. 
00072 
00073   HIPS_SetOptionINT(HIPS_id,HIPS_REORDER,List_.get("hips: reorder",1));
00074   HIPS_SetOptionINT(HIPS_id,HIPS_GRAPH_SYM,List_.get("hips: graph symmetric",0));
00075   HIPS_SetOptionINT(HIPS_id,HIPS_FORTRAN_NUMBERING,List_.get("hips: fortran numbering",0));
00076   HIPS_SetOptionINT(HIPS_id,HIPS_DOF,List_.get("hips: dof per node",1));
00077 
00078   // This will disable the GMRES wrap around HIPS.
00079   //  HIPS_SetOptionINT(HIPS_id,HIPS_ITMAX,1);
00080   HIPS_SetOptionINT(HIPS_id,HIPS_ITMAX,-1);  
00081   //  HIPS_SetOptionINT(HIPS_id,HIPS_KRYLOV_METHOD,List_.get("hips: krylov",0));
00082   HIPS_SetOptionINT(HIPS_id,HIPS_KRYLOV_RESTART,-1);
00083   
00084   // Make sure the ILU always runs, by setting the internal tolerance really, really low.
00085   HIPS_SetOptionREAL(HIPS_id,HIPS_PREC,1e-100);
00086 
00087   // Options for Iterative only
00088   if(List_.get("hips: strategy",HIPS_ITERATIVE)==HIPS_ITERATIVE){  
00089     HIPS_SetOptionREAL(HIPS_id, HIPS_DROPTOL0, List_.get("hips: drop tolerance",1e-2));
00090     HIPS_SetOptionREAL(HIPS_id, HIPS_DROPTOL1, List_.get("hips: drop tolerance",1e-2));
00091     HIPS_SetOptionREAL(HIPS_id, HIPS_DROPTOLE, List_.get("hips: drop tolerance",1e-2));    
00092   }
00093   // NTS: This is only a subset of the actual HIPS options. 
00094   return 0;
00095 }
00096 
00097 
00098 int Ifpack_HIPS::Compute(){
00099   if(HIPS_id==-1) IFPACK_CHK_ERR(-1);
00100   int N=A_->NumMyRows(), nnz=A_->NumMyNonzeros();
00101   const Epetra_Comm &Comm=A_->Comm();
00102   int mypid=Comm.MyPID();
00103   
00104   // Pull the column indices, if possible
00105   int *rowptr,*colind,ierr,maxnr,Nr;
00106   double *values;  
00107   Epetra_CrsMatrix *Acrs=dynamic_cast<Epetra_CrsMatrix*>(&*A_);
00108   const Epetra_Map &RowMap=A_->RowMatrixRowMap();
00109   const Epetra_Map &ColMap=A_->RowMatrixColMap();
00110   if(Acrs) Acrs->ExtractCrsDataPointers(rowptr,colind,values);
00111   else{
00112     maxnr=A_->MaxNumEntries();
00113     colind=new int[maxnr];
00114     values=new double[maxnr];
00115   }
00116 
00117   // Create 0-to-N-1 consistent maps for rows and columns
00118   RowMap0_=rcp(new Epetra_Map(-1,N,0,Comm));
00119   Epetra_IntVector RowGIDs(View,RowMap,RowMap0_->MyGlobalElements());
00120   Epetra_IntVector ColGIDs(ColMap);
00121   Epetra_Import RowToCol(ColMap,RowMap);
00122   ColGIDs.Import(RowGIDs,RowToCol,Insert);
00123   ColMap0_=rcp(new Epetra_Map(-1,ColMap.NumMyElements(),ColGIDs.Values(),0,Comm));
00124 
00125   int *gcolind=0;
00126   if(Acrs){
00127     //  Global CIDs
00128     gcolind=new int[nnz];
00129     for(int j=0;j<nnz;j++) gcolind[j]=RowMap0_->GID(colind[j]);        
00130     ierr =  HIPS_GraphDistrCSR(HIPS_id,A_->NumGlobalRows(),A_->NumMyRows(),RowMap0_->MyGlobalElements(),
00131                                rowptr,gcolind);  
00132   }
00133   else{
00134     // Do things the hard way
00135     ierr=HIPS_GraphBegin(HIPS_id,A_->NumGlobalRows(),nnz);
00136     if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-2);
00137 
00138     // Graph insert - RM mode
00139     for(int i=0;i<N;i++){
00140       A_->ExtractMyRowCopy(i,maxnr,Nr,values,colind);
00141       for(int j=0;j<Nr;j++){
00142         ierr=HIPS_GraphEdge(HIPS_id,RowMap0_->GID(i),ColMap0_->GID(colind[j]));
00143         if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-3);        
00144       }
00145     }
00146     ierr=HIPS_GraphEnd(HIPS_id);
00147   }      
00148   if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-4);  
00149   
00150   /*Have processor 0 send in the partition*/
00151   // NTS: This is really, really annoying.  Look at all this import/export
00152   // stuff.  This is mind-numbingly unnecessary.
00153   
00154   Epetra_Map OnePerProcMap(-1,1,0,Comm);
00155   Epetra_IntVector RowsPerProc(OnePerProcMap);
00156   Epetra_IntVector RowGID(View,*RowMap0_,RowMap0_->MyGlobalElements());
00157   
00158   // Get the RPP partial sums
00159   Comm.ScanSum(&N,&(RowsPerProc[0]),1); 
00160 
00161   // Build the maps for xfer to proc 0
00162   int OPP_els=0,RPP_els=0;
00163   if(!mypid){OPP_els=Comm.NumProc(); RPP_els=A_->NumGlobalRows();}
00164   Epetra_Map OPPMap_0(-1,OPP_els,0,Comm);
00165   Epetra_Map RPPMap_0(-1,RPP_els,0,Comm);
00166   Epetra_Import OPP_importer(OPPMap_0,OnePerProcMap);
00167   Epetra_Import RPP_importer(RPPMap_0,*RowMap0_);
00168   
00169   // Pull the vectors to proc 0
00170   Epetra_IntVector OPP_0(OPPMap_0);
00171   Epetra_IntVector RPP_0(RPPMap_0);
00172   OPP_0.Import(RowsPerProc,OPP_importer,Add);
00173   RPP_0.Import(RowGID,RPP_importer,Add);
00174  
00175   // Setup the partition
00176   if(!mypid){    
00177     int *mapptr=0;
00178     mapptr=new int[Comm.NumProc()+1];
00179     mapptr[0]=0;
00180     for(int i=0;i<Comm.NumProc();i++)
00181       mapptr[i+1]=OPP_0[i];
00182 
00183     // Call is only necessary on proc 0
00184     ierr=HIPS_SetPartition(HIPS_id,A_->Comm().NumProc(),mapptr,RPP_0.Values());
00185     HIPS_ExitOnError(ierr);
00186     delete [] mapptr;
00187   }      
00188 
00189   if(Acrs)
00190     ierr = HIPS_MatrixDistrCSR(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),rowptr,gcolind,values,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_FOOL,0);  
00191   else{
00192     // Do things the hard way
00193      ierr = HIPS_AssemblyBegin(HIPS_id, nnz, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_FOOL,0);
00194      if(ierr!=HIPS_SUCCESS){HIPS_PrintError(ierr);IFPACK_CHK_ERR(-5);}
00195 
00196      // Matrix insert - RM Mode
00197      for(int i=0;i<N;i++){
00198        A_->ExtractMyRowCopy(i,maxnr,Nr,values,colind);
00199        for(int j=0;j<Nr;j++){      
00200          ierr = HIPS_AssemblySetValue(HIPS_id,RowMap0_->GID(i),ColMap0_->GID(colind[j]), values[j]);
00201          if(ierr!=HIPS_SUCCESS){HIPS_PrintError(ierr);IFPACK_CHK_ERR(-6);}
00202        }
00203      }
00204      ierr = HIPS_AssemblyEnd(HIPS_id);
00205   }
00206   if(ierr!=HIPS_SUCCESS){HIPS_PrintError(ierr);IFPACK_CHK_ERR(-7);}
00207 
00208   // Force factorization
00209   //NTS: This is odd.  There should be a better way to force this to happen.
00210   double *X=new double[A_->NumMyRows()];
00211   double *Y=new double[A_->NumMyRows()];
00212   for(int i=0;i<A_->NumMyRows();i++) X[i]=1.0;
00213 
00214   // Force HIPS to do it's own Import/Export
00215   ierr=HIPS_SetRHS(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),X,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_FOOL);
00216   if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-11);
00217   
00218   ierr=HIPS_GetSolution(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),Y,HIPS_ASSEMBLY_FOOL);
00219   if(ierr!=HIPS_SUCCESS) {
00220     HIPS_PrintError(ierr);
00221     IFPACK_CHK_ERR(-12);
00222   }
00223   
00224   // Reset output for iteration
00225   HIPS_SetOptionINT(HIPS_id,HIPS_VERBOSE,List_.get("hips: iteration output",0));
00226 
00227   // Set Label
00228   int nnzP=0;
00229   HIPS_GetInfoINT(HIPS_id,HIPS_INFO_NNZ,&nnzP);
00230   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());
00231   else sprintf(Label_,"Ifpack_HIPS [dt=%4.1e]",List_.get("hips: drop tolerance",1e-2));
00232   // NTS: fill requires a HIPS debug level of at least 2
00233   
00234   IsComputed_=true;
00235 
00236   // Cleanup
00237   if(!Acrs){
00238     delete [] colind;
00239     delete [] values;
00240   }
00241   else
00242     delete [] gcolind;
00243   
00244   delete [] X;
00245   delete [] Y;
00246   return 0;
00247 }
00248 
00249 
00250 
00251 
00252 int Ifpack_HIPS::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
00253   int rv;
00254   if (!IsComputed())
00255     IFPACK_CHK_ERR(-3);
00256 
00257   if (X.NumVectors() != Y.NumVectors())
00258     IFPACK_CHK_ERR(-2);
00259 
00260   // HAQ: For now
00261   if(X.NumVectors()!=1) IFPACK_CHK_ERR(-42);
00262 
00263   // Wrapping for X==Y
00264   Teuchos::RefCountPtr< Epetra_MultiVector > X2;
00265   if (X.Pointers()[0] == Y.Pointers()[0])
00266     X2 = Teuchos::rcp( new Epetra_MultiVector(X) );    
00267   else
00268     X2 = Teuchos::rcp( (Epetra_MultiVector*)&X, false );
00269 
00270   Time_.ResetStartTime();
00271 
00272   // Force HIPS to do it's own Import/Export
00273   rv=HIPS_SetRHS(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),(*X2)[0],HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_FOOL);
00274   if(rv!=HIPS_SUCCESS) IFPACK_CHK_ERR(-11);
00275   
00276   rv=HIPS_GetSolution(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),Y[0],HIPS_ASSEMBLY_FOOL);
00277   if(rv!=HIPS_SUCCESS) {
00278     HIPS_PrintError(rv);
00279     IFPACK_CHK_ERR(-12);
00280   }
00281   
00282   ++NumApplyInverse_;
00283   ApplyInverseTime_ += Time_.ElapsedTime();
00284 
00285   return(0);
00286 }
00287 
00288 
00289 ostream& Ifpack_HIPS::Print(ostream& os) const{
00290   os<<"Need to add meaningful output"<<endl;
00291   return os;
00292 }
00293 
00294 
00295 double Ifpack_HIPS::Condest(const Ifpack_CondestType CT, 
00296                              const int MaxIters,
00297                              const double Tol,
00298                              Epetra_RowMatrix* Matrix_in){
00299   return -1.0;
00300 }
00301 
00302 
00303 
00304 
00305 #endif
 All Classes Files Functions Variables Enumerations Friends