Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
example/ifpack_hb/cxx_main.cpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_ConfigDefs.h"
44 #include <stdio.h>
45 #include <stdlib.h>
46 #include <ctype.h>
47 #include <assert.h>
48 #include <string.h>
49 #include <math.h>
50 #ifdef EPETRA_MPI
51 #include "mpi.h"
52 #include "Epetra_MpiComm.h"
53 #else
54 #include "Epetra_SerialComm.h"
55 #endif
56 #include "Trilinos_Util.h"
57 #include "Epetra_Comm.h"
58 #include "Epetra_Map.h"
59 #include "Epetra_Time.h"
60 #include "Epetra_BlockMap.h"
61 #include "Epetra_MultiVector.h"
62 #include "Epetra_Vector.h"
63 #include "Epetra_Export.h"
64 
65 #include "Epetra_VbrMatrix.h"
66 #include "Epetra_CrsMatrix.h"
67 #include "Ifpack_IlukGraph.h"
68 #include "Ifpack_CrsRiluk.h"
69 
71  Ifpack_CrsRiluk *M,
72  int Maxiter, double Tolerance,
73  double *residual, bool verbose);
74 
75 int main(int argc, char *argv[]) {
76  using std::cout;
77  using std::endl;
78 
79 #ifdef EPETRA_MPI
80  MPI_Init(&argc,&argv);
81  Epetra_MpiComm Comm (MPI_COMM_WORLD);
82 #else
83  Epetra_SerialComm Comm;
84 #endif
85 
86  cout << Comm << endl;
87 
88  int MyPID = Comm.MyPID();
89 
90  bool verbose = false;
91  if (MyPID==0) verbose = true;
92 
93  if(argc < 2 && verbose) {
94  cerr << "Usage: " << argv[0]
95  << " HB_filename [level_fill [level_overlap [absolute_threshold [ relative_threshold]]]]" << endl
96  << "where:" << endl
97  << "HB_filename - filename and path of a Harwell-Boeing data set" << endl
98  << "level_fill - The amount of fill to use for ILU(k) preconditioner (default 0)" << endl
99  << "level_overlap - The amount of overlap used for overlapping Schwarz subdomains (default 0)" << endl
100  << "absolute_threshold - The minimum value to place on the diagonal prior to factorization (default 0.0)" << endl
101  << "relative_threshold - The relative amount to perturb the diagonal prior to factorization (default 1.0)" << endl << endl
102  << "To specify a non-default value for one of these parameters, you must specify all" << endl
103  << " preceding values but not any subsequent parameters. Example:" << endl
104  << "ifpackHbSerialMsr.exe mymatrix.hb 1 - loads mymatrix.hb, uses level fill of one, all other values are defaults" << endl
105  << endl;
106  return(1);
107 
108  }
109 
110  // Uncomment the next three lines to debug in mpi mode
111  //int tmp;
112  //if (MyPID==0) cin >> tmp;
113  //Comm.Barrier();
114 
115  Epetra_Map * readMap;
116  Epetra_CrsMatrix * readA;
117  Epetra_Vector * readx;
118  Epetra_Vector * readb;
119  Epetra_Vector * readxexact;
120 
121  // Call routine to read in HB problem
122  Trilinos_Util_ReadHb2Epetra(argv[1], Comm, readMap, readA, readx, readb, readxexact);
123 
124  // Create uniform distributed map
125  Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);
126 
127  // Create Exporter to distribute read-in matrix and vectors
128 
129  Epetra_Export exporter(*readMap, map);
130  Epetra_CrsMatrix A(Copy, map, 0);
131  Epetra_Vector x(map);
132  Epetra_Vector b(map);
133  Epetra_Vector xexact(map);
134 
135  Epetra_Time FillTimer(Comm);
136  x.Export(*readx, exporter, Add);
137  b.Export(*readb, exporter, Add);
138  xexact.Export(*readxexact, exporter, Add);
139  Comm.Barrier();
140  double vectorRedistributeTime = FillTimer.ElapsedTime();
141  A.Export(*readA, exporter, Add);
142  Comm.Barrier();
143  double matrixRedistributeTime = FillTimer.ElapsedTime() - vectorRedistributeTime;
144  assert(A.FillComplete()==0);
145  Comm.Barrier();
146  double fillCompleteTime = FillTimer.ElapsedTime() - matrixRedistributeTime;
147  if (Comm.MyPID()==0) {
148  cout << "\n\n****************************************************" << endl;
149  cout << "\n Vector redistribute time (sec) = " << vectorRedistributeTime<< endl;
150  cout << " Matrix redistribute time (sec) = " << matrixRedistributeTime << endl;
151  cout << " Transform to Local time (sec) = " << fillCompleteTime << endl<< endl;
152  }
153  Epetra_Vector tmp1(*readMap);
154  Epetra_Vector tmp2(map);
155  readA->Multiply(false, *readxexact, tmp1);
156 
157  A.Multiply(false, xexact, tmp2);
158  double residual;
159  tmp1.Norm2(&residual);
160  if (verbose) cout << "Norm of Ax from file = " << residual << endl;
161  tmp2.Norm2(&residual);
162  if (verbose) cout << "Norm of Ax after redistribution = " << residual << endl << endl << endl;
163 
164  //cout << "A from file = " << *readA << endl << endl << endl;
165 
166  //cout << "A after dist = " << A << endl << endl << endl;
167 
168  delete readA;
169  delete readx;
170  delete readb;
171  delete readxexact;
172  delete readMap;
173 
174  Comm.Barrier();
175 
176  // Construct ILU preconditioner
177 
178  double elapsed_time, total_flops, MFLOPs;
179  Epetra_Time timer(Comm);
180 
181  int LevelFill = 0;
182  if (argc > 2) LevelFill = atoi(argv[2]);
183  if (verbose) cout << "Using Level Fill = " << LevelFill << endl;
184  int Overlap = 0;
185  if (argc > 3) Overlap = atoi(argv[3]);
186  if (verbose) cout << "Using Level Overlap = " << Overlap << endl;
187  double Athresh = 0.0;
188  if (argc > 4) Athresh = atof(argv[4]);
189  if (verbose) cout << "Using Absolute Threshold Value of = " << Athresh << endl;
190 
191  double Rthresh = 1.0;
192  if (argc > 5) Rthresh = atof(argv[5]);
193  if (verbose) cout << "Using Relative Threshold Value of = " << Rthresh << endl;
194 
195  Ifpack_IlukGraph * IlukGraph = 0;
196  Ifpack_CrsRiluk * ILUK = 0;
197 
198  if (LevelFill>-1) {
199  elapsed_time = timer.ElapsedTime();
200  IlukGraph = new Ifpack_IlukGraph(A.Graph(), LevelFill, Overlap);
201  assert(IlukGraph->ConstructFilledGraph()==0);
202  elapsed_time = timer.ElapsedTime() - elapsed_time;
203  if (verbose) cout << "Time to construct ILUK graph = " << elapsed_time << endl;
204 
205 
206  Epetra_Flops fact_counter;
207 
208  elapsed_time = timer.ElapsedTime();
209  ILUK = new Ifpack_CrsRiluk(*IlukGraph);
210  ILUK->SetFlopCounter(fact_counter);
211  ILUK->SetAbsoluteThreshold(Athresh);
212  ILUK->SetRelativeThreshold(Rthresh);
213  //assert(ILUK->InitValues()==0);
214  int initerr = ILUK->InitValues(A);
215  if (initerr!=0) cout << Comm << "InitValues error = " << initerr;
216  assert(ILUK->Factor()==0);
217  elapsed_time = timer.ElapsedTime() - elapsed_time;
218  total_flops = ILUK->Flops();
219  MFLOPs = total_flops/elapsed_time/1000000.0;
220  if (verbose) cout << "Time to compute preconditioner values = "
221  << elapsed_time << endl
222  << "MFLOPS for Factorization = " << MFLOPs << endl;
223  //cout << *ILUK << endl;
224  }
225  double Condest;
226  ILUK->Condest(false, Condest);
227 
228  if (verbose) cout << "Condition number estimate for this preconditioner = " << Condest << endl;
229  int Maxiter = 500;
230  double Tolerance = 1.0E-14;
231 
232  Epetra_Vector xcomp(map);
233  Epetra_Vector resid(map);
234 
235  Epetra_Flops counter;
236  A.SetFlopCounter(counter);
237  xcomp.SetFlopCounter(A);
238  b.SetFlopCounter(A);
239  resid.SetFlopCounter(A);
240  ILUK->SetFlopCounter(A);
241 
242  elapsed_time = timer.ElapsedTime();
243 
244  BiCGSTAB(A, xcomp, b, ILUK, Maxiter, Tolerance, &residual, verbose);
245 
246  elapsed_time = timer.ElapsedTime() - elapsed_time;
247  total_flops = counter.Flops();
248  MFLOPs = total_flops/elapsed_time/1000000.0;
249  if (verbose) cout << "Time to compute solution = "
250  << elapsed_time << endl
251  << "Number of operations in solve = " << total_flops << endl
252  << "MFLOPS for Solve = " << MFLOPs<< endl << endl;
253 
254  resid.Update(1.0, xcomp, -1.0, xexact, 0.0); // resid = xcomp - xexact
255 
256  resid.Norm2(&residual);
257 
258  if (verbose) cout << "Norm of the difference between exact and computed solutions = " << residual << endl;
259 
260 
261 
262 
263  if (ILUK!=0) delete ILUK;
264  if (IlukGraph!=0) delete IlukGraph;
265 
266 #ifdef EPETRA_MPI
267  MPI_Finalize() ;
268 #endif
269 
270 return 0 ;
271 }
273  Epetra_Vector &x,
274  Epetra_Vector &b,
275  Ifpack_CrsRiluk *M,
276  int Maxiter,
277  double Tolerance,
278  double *residual, bool verbose) {
279 
280  // Allocate vectors needed for iterations
281  Epetra_Vector phat(x.Map()); phat.SetFlopCounter(x);
282  Epetra_Vector p(x.Map()); p.SetFlopCounter(x);
283  Epetra_Vector shat(x.Map()); shat.SetFlopCounter(x);
284  Epetra_Vector s(x.Map()); s.SetFlopCounter(x);
285  Epetra_Vector r(x.Map()); r.SetFlopCounter(x);
286  Epetra_Vector rtilde(x.Map()); rtilde.Random(); rtilde.SetFlopCounter(x);
287  Epetra_Vector v(x.Map()); v.SetFlopCounter(x);
288 
289 
290  A.Multiply(false, x, r); // r = A*x
291 
292  r.Update(1.0, b, -1.0); // r = b - A*x
293 
294  double r_norm, b_norm, scaled_r_norm, rhon, rhonm1 = 1.0;
295  double alpha = 1.0, omega = 1.0, sigma;
296  double omega_num, omega_den;
297  r.Norm2(&r_norm);
298  b.Norm2(&b_norm);
299  scaled_r_norm = r_norm/b_norm;
300  r.Dot(rtilde,&rhon);
301 
302  if (verbose) cout << "Initial residual = " << r_norm
303  << " Scaled residual = " << scaled_r_norm << endl;
304 
305 
306  for (int i=0; i<Maxiter; i++) { // Main iteration loop
307 
308  double beta = (rhon/rhonm1) * (alpha/omega);
309  rhonm1 = rhon;
310 
311  /* p = r + beta*(p - omega*v) */
312  /* phat = M^-1 p */
313  /* v = A phat */
314 
315  double dtemp = - beta*omega;
316 
317  p.Update(1.0, r, dtemp, v, beta);
318  if (M==0)
319  phat.Scale(1.0, p);
320  else
321  M->Solve(false, p, phat);
322  A.Multiply(false, phat, v);
323 
324 
325  rtilde.Dot(v,&sigma);
326  alpha = rhon/sigma;
327 
328  /* s = r - alpha*v */
329  /* shat = M^-1 s */
330  /* r = A shat (r is a tmp here for t ) */
331 
332  s.Update(-alpha, v, 1.0, r, 0.0);
333  if (M==0)
334  shat.Scale(1.0, s);
335  else
336  M->Solve(false, s, shat);
337  A.Multiply(false, shat, r);
338 
339  r.Dot(s, &omega_num);
340  r.Dot(r, &omega_den);
341  omega = omega_num / omega_den;
342 
343  /* x = x + alpha*phat + omega*shat */
344  /* r = s - omega*r */
345 
346  x.Update(alpha, phat, omega, shat, 1.0);
347  r.Update(1.0, s, -omega);
348 
349  r.Norm2(&r_norm);
350  scaled_r_norm = r_norm/b_norm;
351  r.Dot(rtilde,&rhon);
352 
353  if (verbose) cout << "Iter "<< i << " residual = " << r_norm
354  << " Scaled residual = " << scaled_r_norm << endl;
355 
356  if (scaled_r_norm < Tolerance) break;
357  }
358  return;
359 }
int NumGlobalElements() const
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
double ElapsedTime(void) const
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int MyPID() const
void SetAbsoluteThreshold(double Athresh)
Set absolute threshold value.
int FillComplete(bool OptimizeDataStorage=true)
void SetRelativeThreshold(double Rthresh)
Set relative threshold value.
virtual const Epetra_BlockMap & Map() const =0
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
int main(int argc, char *argv[])
int Solve(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_CrsRiluk forward/back solve on a Epetra_MultiVector X in Y (works for ...
void BiCGSTAB(Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, Ifpack_CrsRiluk *M, int Maxiter, double Tolerance, double *residual, bool verbose)
double Flops() const
void Barrier() const
const Epetra_CrsGraph & Graph() const
virtual int ConstructFilledGraph()
Does the actual construction of the graph.
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...