Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_HIPS.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // Ifpack: Object-Oriented Algebraic Preconditioner Package
6 // Copyright (2002) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #include "Ifpack_HIPS.h"
45 #if defined(HAVE_IFPACK_HIPS) && defined(HAVE_MPI)
46 
47 
48 #include "Ifpack_Utils.h"
49 extern "C" {
50 #include "hips.h"
51 }
52 //#include "io.h"
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;
60 #endif
61 
62 using Teuchos::RefCountPtr;
63 using Teuchos::rcp;
64 
65 
66 
67 Ifpack_HIPS::Ifpack_HIPS(Epetra_RowMatrix* A):
68  A_(rcp(A,false)),
69  HIPS_id(-1), // Assumes user will initialze HIPS outside
70  IsParallel_(false),
71  IsInitialized_(false),
72  IsComputed_(false),
73  Label_(),
74  Time_(A_->Comm())
75 {
76 
77 }
78 
79 void Ifpack_HIPS::Destroy(){
80  //NTS: Assume user will call HIPS_Finalize elsewhere - HIPS_Clean never needs
81  //to be called if HIPS_Finalize is called at the end, unless you want to reuse
82  //a slot.
83 }
84 
85 
86 
87 int Ifpack_HIPS::Initialize(){
88  if(Comm().NumProc() != 1) IsParallel_ = true;
89  else IsParallel_ = false;
90  IsInitialized_=true;
91  return 0;
92 }
93 
94 int Ifpack_HIPS::SetParameters(Teuchos::ParameterList& parameterlist){
95  List_=parameterlist;
96 
97  // Grab the hips ID
98  HIPS_id=List_.get("hips: id",-1);
99  if(HIPS_id==-1) IFPACK_CHK_ERR(-1);
100 
101  // Set Defaults
102  HIPS_SetDefaultOptions(HIPS_id,List_.get("hips: strategy",HIPS_ITERATIVE));
103 
104  // Set the communicator
105  const Epetra_MpiComm* MpiComm=dynamic_cast<const Epetra_MpiComm*>(&A_->Comm());
106  if(!MpiComm) IFPACK_CHK_ERR(-2);
107  HIPS_SetCommunicator(HIPS_id,MpiComm->GetMpiComm());
108 
109  // Options
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));
113  // Sadly, this fill option doesn't work for HIPS_ITERATIVE mode, meaning the
114  // only way to control fill-in is via the drop tolerance.
115 
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));
120 
121  // This will disable the GMRES wrap around HIPS.
122  // HIPS_SetOptionINT(HIPS_id,HIPS_ITMAX,1);
123  HIPS_SetOptionINT(HIPS_id,HIPS_ITMAX,-1);
124  // HIPS_SetOptionINT(HIPS_id,HIPS_KRYLOV_METHOD,List_.get("hips: krylov",0));
125  HIPS_SetOptionINT(HIPS_id,HIPS_KRYLOV_RESTART,-1);
126 
127  // Make sure the ILU always runs, by setting the internal tolerance really, really low.
128  HIPS_SetOptionREAL(HIPS_id,HIPS_PREC,1e-100);
129 
130  // Options for Iterative only
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));
135  }
136  // NTS: This is only a subset of the actual HIPS options.
137  return 0;
138 }
139 
140 
141 int Ifpack_HIPS::Compute(){
142  if(HIPS_id==-1) IFPACK_CHK_ERR(-1);
143  int N=A_->NumMyRows(), nnz=A_->NumMyNonzeros();
144  const Epetra_Comm &Comm=A_->Comm();
145  int mypid=Comm.MyPID();
146 
147  // Pull the column indices, if possible
148  int *rowptr,*colind,ierr,maxnr,Nr;
149  double *values;
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);
154  else{
155  maxnr=A_->MaxNumEntries();
156  colind=new int[maxnr];
157  values=new double[maxnr];
158  }
159 
160  // Create 0-to-N-1 consistent maps for rows and columns
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));
167 
168  int *gcolind=0;
169  if(Acrs){
170  // Global CIDs
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(),
174  rowptr,gcolind);
175  }
176  else{
177  // Do things the hard way
178  ierr=HIPS_GraphBegin(HIPS_id,A_->NumGlobalRows(),nnz);
179  if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-2);
180 
181  // Graph insert - RM mode
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);
187  }
188  }
189  ierr=HIPS_GraphEnd(HIPS_id);
190  }
191  if(ierr!=HIPS_SUCCESS) IFPACK_CHK_ERR(-4);
192 
193  /*Have processor 0 send in the partition*/
194  // NTS: This is really, really annoying. Look at all this import/export
195  // stuff. This is mind-numbingly unnecessary.
196 
197  Epetra_Map OnePerProcMap(-1,1,0,Comm);
198  Epetra_IntVector RowsPerProc(OnePerProcMap);
199  Epetra_IntVector RowGID(View,*RowMap0_,RowMap0_->MyGlobalElements());
200 
201  // Get the RPP partial sums
202  Comm.ScanSum(&N,&(RowsPerProc[0]),1);
203 
204  // Build the maps for xfer to proc 0
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_);
211 
212  // Pull the vectors to proc 0
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);
217 
218  // Setup the partition
219  if(!mypid){
220  int *mapptr=0;
221  mapptr=new int[Comm.NumProc()+1];
222  mapptr[0]=0;
223  for(int i=0;i<Comm.NumProc();i++)
224  mapptr[i+1]=OPP_0[i];
225 
226  // Call is only necessary on proc 0
227  ierr=HIPS_SetPartition(HIPS_id,A_->Comm().NumProc(),mapptr,RPP_0.Values());
228  HIPS_ExitOnError(ierr);
229  delete [] mapptr;
230  }
231 
232  if(Acrs)
233  ierr = HIPS_MatrixDistrCSR(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),rowptr,gcolind,values,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_OVW,HIPS_ASSEMBLY_FOOL,0);
234  else{
235  // Do things the hard way
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);}
238 
239  // Matrix insert - RM Mode
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);}
245  }
246  }
247  ierr = HIPS_AssemblyEnd(HIPS_id);
248  }
249  if(ierr!=HIPS_SUCCESS){HIPS_PrintError(ierr);IFPACK_CHK_ERR(-7);}
250 
251  // Force factorization
252  //NTS: This is odd. There should be a better way to force this to happen.
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;
256 
257  // Force HIPS to do it's own Import/Export
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);
260 
261  ierr=HIPS_GetSolution(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),Y,HIPS_ASSEMBLY_FOOL);
262  if(ierr!=HIPS_SUCCESS) {
263  HIPS_PrintError(ierr);
264  IFPACK_CHK_ERR(-12);
265  }
266 
267  // Reset output for iteration
268  HIPS_SetOptionINT(HIPS_id,HIPS_VERBOSE,List_.get("hips: iteration output",0));
269 
270  // Set Label
271  int nnzP=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));
275  // NTS: fill requires a HIPS debug level of at least 2
276 
277  IsComputed_=true;
278 
279  // Cleanup
280  if(!Acrs){
281  delete [] colind;
282  delete [] values;
283  }
284  else
285  delete [] gcolind;
286 
287  delete [] X;
288  delete [] Y;
289  return 0;
290 }
291 
292 
293 
294 
295 int Ifpack_HIPS::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
296  int rv;
297  if (!IsComputed())
298  IFPACK_CHK_ERR(-3);
299 
300  if (X.NumVectors() != Y.NumVectors())
301  IFPACK_CHK_ERR(-2);
302 
303  // HAQ: For now
304  if(X.NumVectors()!=1) IFPACK_CHK_ERR(-42);
305 
306  // Wrapping for X==Y
307  Teuchos::RefCountPtr< Epetra_MultiVector > X2;
308  if (X.Pointers()[0] == Y.Pointers()[0])
309  X2 = Teuchos::rcp( new Epetra_MultiVector(X) );
310  else
311  X2 = Teuchos::rcp( (Epetra_MultiVector*)&X, false );
312 
313  Time_.ResetStartTime();
314 
315  // Force HIPS to do it's own Import/Export
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);
318 
319  rv=HIPS_GetSolution(HIPS_id,A_->NumMyRows(),RowMap0_->MyGlobalElements(),Y[0],HIPS_ASSEMBLY_FOOL);
320  if(rv!=HIPS_SUCCESS) {
321  HIPS_PrintError(rv);
322  IFPACK_CHK_ERR(-12);
323  }
324 
325  ++NumApplyInverse_;
326  ApplyInverseTime_ += Time_.ElapsedTime();
327 
328  return(0);
329 }
330 
331 
332 std::ostream& Ifpack_HIPS::Print(std::ostream& os) const{
333  os<<"Need to add meaningful output"<<endl;
334  return os;
335 }
336 
337 
338 double Ifpack_HIPS::Condest(const Ifpack_CondestType CT,
339  const int MaxIters,
340  const double Tol,
341  Epetra_RowMatrix* Matrix_in){
342  return -1.0;
343 }
344 
345 
346 
347 
348 #endif
int MyGlobalElements(int *MyGlobalElementList) const
T & get(ParameterList &l, const std::string &name)
virtual int MyPID() const =0
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
int NumMyElements() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const int N
#define false
virtual int NumProc() const =0
#define IFPACK_CHK_ERR(ifpack_err)
virtual int ScanSum(double *MyVals, double *ScanSums, int Count) const =0
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const