Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
72 int Maxiter, double Tolerance,
73 double *residual, bool verbose);
74
75int 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
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
270return 0 ;
271}
273 Epetra_Vector &x,
274 Epetra_Vector &b,
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}
Add
Copy
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
int FillComplete(bool OptimizeDataStorage=true)
const Epetra_CrsGraph & Graph() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
const Epetra_BlockMap & Map() const
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
void Barrier() const
int MyPID() const
int Dot(const Epetra_MultiVector &A, double *Result) const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Norm2(double *Result) const
double ElapsedTime(void) const
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
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 ...
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int main(int argc, char *argv[])