Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
49extern "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
59extern int ML_NODE_ID;
60#endif
61
62using Teuchos::RefCountPtr;
63using Teuchos::rcp;
64
65
66
67Ifpack_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
79void 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
87int Ifpack_HIPS::Initialize(){
88 if(Comm().NumProc() != 1) IsParallel_ = true;
89 else IsParallel_ = false;
90 IsInitialized_=true;
91 return 0;
92}
93
94int 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
141int 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
295int 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
332std::ostream& Ifpack_HIPS::Print(std::ostream& os) const{
333 os<<"Need to add meaningful output"<<endl;
334 return os;
335}
336
337
338double 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
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
#define IFPACK_CHK_ERR(ifpack_err)
int NumMyElements() const
virtual int NumProc() const=0
virtual int MyPID() const=0
virtual int ScanSum(double *MyVals, double *ScanSums, int Count) const=0
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
MPI_Comm GetMpiComm() const
int NumVectors() const
double ** Pointers() const
#define false
const int N