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);
75 int main(
int argc,
char *argv[]) {
80 MPI_Init(&argc,&argv);
88 int MyPID = Comm.
MyPID();
91 if (MyPID==0) verbose =
true;
93 if(argc < 2 && verbose) {
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);
136 x.Export(*readx, exporter,
Add);
137 b.Export(*readb, exporter,
Add);
138 xexact.Export(*readxexact, exporter,
Add);
140 double vectorRedistributeTime = FillTimer.
ElapsedTime();
141 A.Export(*readA, exporter,
Add);
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;
210 ILUK->SetFlopCounter(fact_counter);
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;
236 A.SetFlopCounter(counter);
237 xcomp.SetFlopCounter(A);
239 resid.SetFlopCounter(A);
240 ILUK->SetFlopCounter(A);
244 BiCGSTAB(A, xcomp, b, ILUK, Maxiter, Tolerance, &residual, verbose);
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) {
292 r.Update(1.0, b, -1.0);
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);
347 r.Update(1.0, s, -omega);
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
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...
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)
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...