61#include "Trilinos_Util.h"
62#include "Ifpack_IlukGraph.h"
63#include "Ifpack_CrsRiluk.h"
67 int Maxiter,
double Tolerance,
68 double *residual,
bool verbose);
70int main(
int argc,
char *argv[]) {
73 MPI_Init(&argc,&argv);
81 int MyPID = Comm.
MyPID();
85 if (MyPID==0) verbose =
true;
87 if(argc < 2 && verbose) {
88 cerr <<
"Usage: " << argv[0]
89 <<
" HPC_filename [level_fill [level_overlap [absolute_threshold [ relative_threshold]]]]" << endl
91 <<
"HB_filename - filename and path of a Harwell-Boeing data set" << endl
92 <<
"level_fill - The amount of fill to use for ILU(k) preconditioner (default 0)" << endl
93 <<
"level_overlap - The amount of overlap used for overlapping Schwarz subdomains (default 0)" << endl
94 <<
"absolute_threshold - The minimum value to place on the diagonal prior to factorization (default 0.0)" << endl
95 <<
"relative_threshold - The relative amount to perturb the diagonal prior to factorization (default 1.0)" << endl << endl
96 <<
"To specify a non-default value for one of these parameters, you must specify all" << endl
97 <<
" preceding values but not any subsequent parameters. Example:" << endl
98 <<
"ifpackHpcSerialMsr.exe mymatrix.hpc 1 - loads mymatrix.hpc, uses level fill of one, all other values are defaults" << endl
106 if (MyPID==0) cin >> tmp;
115 Trilinos_Util_ReadHpc2Epetra(argv[1], Comm, map, A, x, b, xexact);
117 bool smallProblem =
false;
121 cout <<
"Original Matrix = " << endl << *A << endl;
130 cout <<
"Inf norm of Original Matrix = " << A->
NormInf() << endl
131 <<
"Inf norm of Original RHS = " << normb << endl;
136 double reduceInitTime = ReductionTimer.
ElapsedTime();
139 double reduceAnalyzeTime = ReductionTimer.
ElapsedTime() - reduceInitTime;
142 cout <<
"Singletons found" << endl;
144 cout <<
"Singletons not found" << endl;
149 double reduceConstructTime = ReductionTimer.
ElapsedTime() - reduceInitTime;
151 double totalReduceTime = ReductionTimer.
ElapsedTime();
154 cout <<
"\n\n****************************************************" << endl
155 <<
" Reduction init time (sec) = " << reduceInitTime<< endl
156 <<
" Reduction Analyze time (sec) = " << reduceAnalyzeTime << endl
157 <<
" Construct Reduced Problem time (sec) = " << reduceConstructTime << endl
158 <<
" Reduction Total time (sec) = " << totalReduceTime << endl<< endl;
170 cout <<
" Reduced Matrix = " << endl << *Ap << endl
171 <<
" LHS before sol = " << endl << *xp << endl
172 <<
" RHS = " << endl << *bp << endl;
176 double elapsed_time, total_flops, MFLOPs;
180 if (argc > 2) LevelFill = atoi(argv[2]);
181 if (verbose) cout <<
"Using Level Fill = " << LevelFill << endl;
183 if (argc > 3) Overlap = atoi(argv[3]);
184 if (verbose) cout <<
"Using Level Overlap = " << Overlap << endl;
185 double Athresh = 0.0;
186 if (argc > 4) Athresh = atof(argv[4]);
187 if (verbose) cout <<
"Using Absolute Threshold Value of = " << Athresh << endl;
189 double Rthresh = 1.0;
190 if (argc > 5) Rthresh = atof(argv[5]);
191 if (verbose) cout <<
"Using Relative Threshold Value of = " << Rthresh << endl;
193 Ifpack_IlukGraph * IlukGraph = 0;
194 Ifpack_CrsRiluk * ILUK = 0;
198 IlukGraph =
new Ifpack_IlukGraph(Ap->
Graph(), LevelFill, Overlap);
199 assert(IlukGraph->ConstructFilledGraph()==0);
201 if (verbose) cout <<
"Time to construct ILUK graph = " << elapsed_time << endl;
207 ILUK =
new Ifpack_CrsRiluk(*IlukGraph);
208 ILUK->SetFlopCounter(fact_counter);
209 ILUK->SetAbsoluteThreshold(Athresh);
210 ILUK->SetRelativeThreshold(Rthresh);
212 int initerr = ILUK->InitValues(*Ap);
214 cout << endl << Comm << endl <<
" InitValues error = " << initerr;
215 if (initerr==1) cout <<
" Zero diagonal found, warning error only";
216 cout << endl << endl;
218 assert(ILUK->Factor()==0);
220 total_flops = ILUK->Flops();
221 MFLOPs = total_flops/elapsed_time/1000000.0;
222 if (verbose) cout <<
"Time to compute preconditioner values = "
223 << elapsed_time << endl
224 <<
"MFLOPS for Factorization = " << MFLOPs << endl;
227 ILUK->Condest(
false, Condest);
229 if (verbose) cout <<
"Condition number estimate for this preconditioner = " << Condest << endl;
232 double Tolerance = 1.0E-8;
238 if (ILUK!=0) ILUK->SetFlopCounter(*Ap);
242 double normreducedb, normreduceda;
246 cout <<
"Inf norm of Reduced Matrix = " << normreduceda << endl
247 <<
"Inf norm of Reduced RHS = " << normreducedb << endl;
250 BiCGSTAB(*Ap, *xp, *bp, ILUK, Maxiter, Tolerance, &residual, verbose);
253 total_flops = counter.
Flops();
254 MFLOPs = total_flops/elapsed_time/1000000.0;
255 if (verbose) cout <<
"Time to compute solution = "
256 << elapsed_time << endl
257 <<
"Number of operations in solve = " << total_flops << endl
258 <<
"MFLOPS for Solve = " << MFLOPs<< endl << endl;
263 cout <<
" Reduced LHS after sol = " << endl << *xp << endl
264 <<
" Full LHS after sol = " << endl << *x << endl
265 <<
" Full Exact LHS = " << endl << *xexact << endl;
269 resid.
Update(1.0, *x, -1.0, *xexact, 0.0);
271 resid.
Norm2(&residual);
272 double normx, normxexact;
274 xexact->
Norm2(&normxexact);
277 cout <<
"2-norm of computed solution = " << normx << endl
278 <<
"2-norm of exact solution = " << normxexact << endl
279 <<
"2-norm of difference between computed and exact solution = " << residual << endl;
281 if (verbose1 && residual>1.0e-5) {
283 cout <<
"Difference between computed and exact solution appears large..." << endl
284 <<
"Computing norm of A times this difference. If this norm is small, then matrix is singular"
287 assert(A->
Multiply(
false, resid, bdiff)==0);
288 assert(bdiff.
Norm2(&residual)==0);
290 cout <<
"2-norm of A times difference between computed and exact solution = " << residual << endl;
294 cout <<
"********************************************************" << endl
295 <<
" Solving again with 2*Ax=2*b" << endl
296 <<
"********************************************************" << endl;
304 cout <<
"Inf norm of Original Matrix = " << norma << endl
305 <<
"Inf norm of Original RHS = " << normb << endl;
306 double updateReducedProblemTime = ReductionTimer.
ElapsedTime();
309 updateReducedProblemTime = ReductionTimer.
ElapsedTime() - updateReducedProblemTime;
311 cout <<
"\n\n****************************************************" << endl
312 <<
" Update Reduced Problem time (sec) = " << updateReducedProblemTime<< endl
313 <<
"****************************************************" << endl;
322 int initerr = ILUK->InitValues(*Ap);
324 cout << endl << Comm << endl <<
" InitValues error = " << initerr;
325 if (initerr==1) cout <<
" Zero diagonal found, warning error only";
326 cout << endl << endl;
328 assert(ILUK->Factor()==0);
330 total_flops = ILUK->Flops();
331 MFLOPs = total_flops/elapsed_time/1000000.0;
332 if (verbose) cout <<
"Time to compute preconditioner values = "
333 << elapsed_time << endl
334 <<
"MFLOPS for Factorization = " << MFLOPs << endl;
336 ILUK->Condest(
false, Condest);
338 if (verbose) cout <<
"Condition number estimate for this preconditioner = " << Condest << endl;
343 cout <<
"Inf norm of Reduced Matrix = " << normreduceda << endl
344 <<
"Inf norm of Reduced RHS = " << normreducedb << endl;
346 BiCGSTAB(*Ap, *xp, *bp, ILUK, Maxiter, Tolerance, &residual, verbose);
349 total_flops = counter.
Flops();
350 MFLOPs = total_flops/elapsed_time/1000000.0;
351 if (verbose) cout <<
"Time to compute solution = "
352 << elapsed_time << endl
353 <<
"Number of operations in solve = " << total_flops << endl
354 <<
"MFLOPS for Solve = " << MFLOPs<< endl << endl;
359 cout <<
" Reduced LHS after sol = " << endl << *xp << endl
360 <<
" Full LHS after sol = " << endl << *x << endl
361 <<
" Full Exact LHS = " << endl << *xexact << endl;
363 resid.
Update(1.0, *x, -1.0, *xexact, 0.0);
365 resid.
Norm2(&residual);
367 xexact->
Norm2(&normxexact);
370 cout <<
"2-norm of computed solution = " << normx << endl
371 <<
"2-norm of exact solution = " << normxexact << endl
372 <<
"2-norm of difference between computed and exact solution = " << residual << endl;
374 if (verbose1 && residual>1.0e-5) {
376 cout <<
"Difference between computed and exact solution appears large..." << endl
377 <<
"Computing norm of A times this difference. If this norm is small, then matrix is singular"
380 assert(A->
Multiply(
false, resid, bdiff)==0);
381 assert(bdiff.
Norm2(&residual)==0);
383 cout <<
"2-norm of A times difference between computed and exact solution = " << residual << endl;
389 if (ILUK!=0)
delete ILUK;
390 if (IlukGraph!=0)
delete IlukGraph;
404 double *residual,
bool verbose) {
420 double r_norm, b_norm, scaled_r_norm, rhon, rhonm1 = 1.0;
421 double alpha = 1.0, omega = 1.0, sigma;
422 double omega_num, omega_den;
425 scaled_r_norm = r_norm/b_norm;
428 if (verbose) cout <<
"Initial residual = " << r_norm
429 <<
" Scaled residual = " << scaled_r_norm << endl;
432 for (
int i=0; i<Maxiter; i++) {
434 double beta = (rhon/rhonm1) * (alpha/omega);
441 double dtemp = - beta*omega;
443 p.
Update(1.0, r, dtemp, v, beta);
447 M->Solve(
false, p, phat);
451 rtilde.
Dot(v,&sigma);
458 s.
Update(-alpha, v, 1.0, r, 0.0);
462 M->Solve(
false, s, shat);
465 r.
Dot(s, &omega_num);
466 r.
Dot(r, &omega_den);
467 omega = omega_num / omega_den;
472 x.
Update(alpha, phat, omega, shat, 1.0);
476 scaled_r_norm = r_norm/b_norm;
479 if (verbose) cout <<
"Iter "<< i <<
" residual = " << r_norm
480 <<
" Scaled residual = " << scaled_r_norm << endl;
482 if (r_norm < Tolerance)
break;
498 cout <<
"Full System characteristics:" << endl << endl
499 <<
" Dimension = " << fnrows << endl
500 <<
" Number of nonzeros = " <<fnnz << endl
501 <<
" Maximum Number of Row Entries = " << fmaxentries << endl << endl
502 <<
"Reduced System characteristics:" << endl << endl
506 <<
"Singleton information: " << endl
507 <<
" Total number of singletons = " << filter.
NumSingletons() << endl
508 <<
" Number of rows with single entries = " << filter.
NumRowSingletons() << endl
509 <<
" Number of columns with single entries " << endl
510 <<
" (that were not already counted as " << endl
512 <<
"Ratios: " << endl
513 <<
" Percent reduction in dimension = " << filter.
RatioOfDimensions()*100.0 << endl
514 <<
" Percent reduction in nonzero count = " << filter.
RatioOfNonzeros()*100.0 << endl << endl;
int NumGlobalElements() const
Number of elements across all processors.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
virtual int MyPID() const =0
Return my process ID.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int GlobalMaxNumEntries() const
Returns the maximum number of nonzero entries across all rows on all processors.
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A <- ScalarConstant * A).
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
int NumGlobalRows() const
Returns the number of global matrix rows.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
double NormInf() const
Returns the infinity norm of the global matrix.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
Epetra_Flops: The Epetra Floating Point Operations Class.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Epetra_LinearProblem: The Epetra Linear Problem Class.
Epetra_MultiVector * GetRHS() const
Get a pointer to the right-hand-side B.
Epetra_RowMatrix * GetMatrix() const
Get a pointer to the matrix A.
Epetra_MultiVector * GetLHS() const
Get a pointer to the left-hand-side X.
Epetra_Map: A class for partitioning vectors and matrices.
Epetra_MpiComm: The Epetra MPI Communication Class.
void Barrier() const
Epetra_MpiComm Barrier function.
int MyPID() const
Return my process ID.
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
virtual int NumGlobalNonzeros() const =0
Returns the number of nonzero entries in the global matrix.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
virtual int NumGlobalRows() const =0
Returns the number of global matrix rows.
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Time: The Epetra Timing Class.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char *argv[])
void BiCGSTAB(Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, Ifpack_CrsRiluk *M, int Maxiter, double Tolerance, double *residual, bool verbose)
int Statistics(const Epetra_CrsSingletonFilter &filter)