52#include "Epetra_MpiComm.h"
54#include "Epetra_SerialComm.h"
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"
65#include "Epetra_VbrMatrix.h"
66#include "Epetra_CrsMatrix.h"
72 int Maxiter,
double Tolerance,
73 double *residual,
bool verbose);
75int main(
int argc,
char *argv[]) {
80 MPI_Init(&argc,&argv);
88 int MyPID = Comm.
MyPID();
94 cerr <<
"Usage: " << argv[0]
95 <<
" HB_filename [level_fill [level_overlap [absolute_threshold [ relative_threshold]]]]" << 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
122 Trilinos_Util_ReadHb2Epetra(argv[1], Comm, readMap, readA, readx, readb, readxexact);
138 xexact.
Export(*readxexact, exporter,
Add);
140 double vectorRedistributeTime = FillTimer.
ElapsedTime();
143 double matrixRedistributeTime = FillTimer.
ElapsedTime() - vectorRedistributeTime;
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;
155 readA->
Multiply(
false, *readxexact, tmp1);
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;
178 double elapsed_time, total_flops, MFLOPs;
182 if (argc > 2) LevelFill = atoi(argv[2]);
183 if (
verbose) cout <<
"Using Level Fill = " << LevelFill << endl;
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;
191 double Rthresh = 1.0;
192 if (argc > 5) Rthresh = atof(argv[5]);
193 if (
verbose) cout <<
"Using Relative Threshold Value of = " << Rthresh << endl;
203 if (
verbose) cout <<
"Time to construct ILUK graph = " << elapsed_time << endl;
215 if (initerr!=0) cout << Comm <<
"InitValues error = " << initerr;
216 assert(ILUK->
Factor()==0);
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;
228 if (
verbose) cout <<
"Condition number estimate for this preconditioner = " << Condest << endl;
230 double Tolerance = 1.0E-14;
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;
254 resid.
Update(1.0, xcomp, -1.0, xexact, 0.0);
256 resid.
Norm2(&residual);
258 if (
verbose) cout <<
"Norm of the difference between exact and computed solutions = " << residual << endl;
263 if (ILUK!=0)
delete ILUK;
264 if (IlukGraph!=0)
delete IlukGraph;
278 double *residual,
bool verbose) {
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;
299 scaled_r_norm = r_norm/b_norm;
302 if (
verbose) cout <<
"Initial residual = " << r_norm
303 <<
" Scaled residual = " << scaled_r_norm << endl;
306 for (
int i=0; i<Maxiter; i++) {
308 double beta = (rhon/rhonm1) * (alpha/omega);
315 double dtemp = - beta*omega;
317 p.
Update(1.0, r, dtemp, v, beta);
321 M->
Solve(
false, p, phat);
325 rtilde.
Dot(v,&sigma);
332 s.
Update(-alpha, v, 1.0, r, 0.0);
336 M->
Solve(
false, s, shat);
339 r.
Dot(s, &omega_num);
340 r.
Dot(r, &omega_den);
341 omega = omega_num / omega_den;
346 x.
Update(alpha, phat, omega, shat, 1.0);
350 scaled_r_norm = r_norm/b_norm;
353 if (
verbose) cout <<
"Iter "<< i <<
" residual = " << r_norm
354 <<
" Scaled residual = " << scaled_r_norm << endl;
356 if (scaled_r_norm < Tolerance)
break;
int NumGlobalElements() const
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)
int Scale(double ScalarValue)
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 Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
void SetAbsoluteThreshold(double Athresh)
Set absolute threshold value.
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
void SetRelativeThreshold(double Rthresh)
Set relative threshold value.
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors.
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...
virtual int ConstructFilledGraph()
Does the actual construction of the graph.
int main(int argc, char *argv[])