Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
example/ReducedLinearProblem/cxx_main.cpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Epetra: Linear Algebra Services Package
5// Copyright 2011 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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 "Epetra_Comm.h"
44#include "Epetra_Map.h"
45#include "Epetra_Time.h"
46#include "Epetra_BlockMap.h"
47#include "Epetra_MultiVector.h"
48#include "Epetra_Vector.h"
49#include "Epetra_Export.h"
52
53#include "Epetra_VbrMatrix.h"
54#include "Epetra_CrsMatrix.h"
55#ifdef EPETRA_MPI
56#include "mpi.h"
57#include "Epetra_MpiComm.h"
58#else
59#include "Epetra_SerialComm.h"
60#endif
61#include "Trilinos_Util.h"
62#include "Ifpack_IlukGraph.h"
63#include "Ifpack_CrsRiluk.h"
64
66 Ifpack_CrsRiluk *M,
67 int Maxiter, double Tolerance,
68 double *residual, bool verbose);
69int Statistics(const Epetra_CrsSingletonFilter & filter);
70int main(int argc, char *argv[]) {
71
72#ifdef EPETRA_MPI
73 MPI_Init(&argc,&argv);
74 Epetra_MpiComm Comm (MPI_COMM_WORLD);
75#else
77#endif
78
79 cout << Comm << endl;
80
81 int MyPID = Comm.MyPID();
82
83 bool verbose = false;
84 bool verbose1 = true;
85 if (MyPID==0) verbose = true;
86
87 if(argc < 2 && verbose) {
88 cerr << "Usage: " << argv[0]
89 << " HPC_filename [level_fill [level_overlap [absolute_threshold [ relative_threshold]]]]" << endl
90 << "where:" << 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
99 << endl;
100 return(1);
101
102 }
103
104 // Uncomment the next three lines to debug in mpi mode
105 int tmp;
106 if (MyPID==0) cin >> tmp;
107 Comm.Barrier();
108
109 Epetra_Map * map;
111 Epetra_Vector * x;
112 Epetra_Vector * b;
113 Epetra_Vector * xexact;
114
115 Trilinos_Util_ReadHpc2Epetra(argv[1], Comm, map, A, x, b, xexact);
116
117 bool smallProblem = false;
118 if (A->RowMap().NumGlobalElements()<100) smallProblem = true;
119
120 if (smallProblem)
121 cout << "Original Matrix = " << endl << *A << endl;
122
123 x->PutScalar(0.0);
124
125 Epetra_LinearProblem FullProblem(A, x, b);
126 double normb, norma;
127 b->NormInf(&normb);
128 norma = A->NormInf();
129 if (verbose)
130 cout << "Inf norm of Original Matrix = " << A->NormInf() << endl
131 << "Inf norm of Original RHS = " << normb << endl;
132
133 Epetra_Time ReductionTimer(Comm);
134 Epetra_CrsSingletonFilter SingletonFilter;
135 Comm.Barrier();
136 double reduceInitTime = ReductionTimer.ElapsedTime();
137 SingletonFilter.Analyze(A);
138 Comm.Barrier();
139 double reduceAnalyzeTime = ReductionTimer.ElapsedTime() - reduceInitTime;
140
141 if (SingletonFilter.SingletonsDetected())
142 cout << "Singletons found" << endl;
143 else {
144 cout << "Singletons not found" << endl;
145 exit(1);
146 }
147 SingletonFilter.ConstructReducedProblem(&FullProblem);
148 Comm.Barrier();
149 double reduceConstructTime = ReductionTimer.ElapsedTime() - reduceInitTime;
150
151 double totalReduceTime = ReductionTimer.ElapsedTime();
152
153 if (verbose)
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;
159
160 Statistics(SingletonFilter);
161
162 Epetra_LinearProblem * ReducedProblem = SingletonFilter.ReducedProblem();
163
164 Epetra_CrsMatrix * Ap = dynamic_cast<Epetra_CrsMatrix *>(ReducedProblem->GetMatrix());
165 Epetra_Vector * bp = (*ReducedProblem->GetRHS())(0);
166 Epetra_Vector * xp = (*ReducedProblem->GetLHS())(0);
167
168
169 if (smallProblem)
170 cout << " Reduced Matrix = " << endl << *Ap << endl
171 << " LHS before sol = " << endl << *xp << endl
172 << " RHS = " << endl << *bp << endl;
173
174 // Construct ILU preconditioner
175
176 double elapsed_time, total_flops, MFLOPs;
177 Epetra_Time timer(Comm);
178
179 int LevelFill = 0;
180 if (argc > 2) LevelFill = atoi(argv[2]);
181 if (verbose) cout << "Using Level Fill = " << LevelFill << endl;
182 int Overlap = 0;
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;
188
189 double Rthresh = 1.0;
190 if (argc > 5) Rthresh = atof(argv[5]);
191 if (verbose) cout << "Using Relative Threshold Value of = " << Rthresh << endl;
192
193 Ifpack_IlukGraph * IlukGraph = 0;
194 Ifpack_CrsRiluk * ILUK = 0;
195
196 if (LevelFill>-1) {
197 elapsed_time = timer.ElapsedTime();
198 IlukGraph = new Ifpack_IlukGraph(Ap->Graph(), LevelFill, Overlap);
199 assert(IlukGraph->ConstructFilledGraph()==0);
200 elapsed_time = timer.ElapsedTime() - elapsed_time;
201 if (verbose) cout << "Time to construct ILUK graph = " << elapsed_time << endl;
202
203
204 Epetra_Flops fact_counter;
205
206 elapsed_time = timer.ElapsedTime();
207 ILUK = new Ifpack_CrsRiluk(*IlukGraph);
208 ILUK->SetFlopCounter(fact_counter);
209 ILUK->SetAbsoluteThreshold(Athresh);
210 ILUK->SetRelativeThreshold(Rthresh);
211 //assert(ILUK->InitValues()==0);
212 int initerr = ILUK->InitValues(*Ap);
213 if (initerr!=0) {
214 cout << endl << Comm << endl << " InitValues error = " << initerr;
215 if (initerr==1) cout << " Zero diagonal found, warning error only";
216 cout << endl << endl;
217 }
218 assert(ILUK->Factor()==0);
219 elapsed_time = timer.ElapsedTime() - elapsed_time;
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;
225 //cout << *ILUK << endl;
226 double Condest;
227 ILUK->Condest(false, Condest);
228
229 if (verbose) cout << "Condition number estimate for this preconditioner = " << Condest << endl;
230 }
231 int Maxiter = 100;
232 double Tolerance = 1.0E-8;
233
234 Epetra_Flops counter;
235 Ap->SetFlopCounter(counter);
236 xp->SetFlopCounter(*Ap);
237 bp->SetFlopCounter(*Ap);
238 if (ILUK!=0) ILUK->SetFlopCounter(*Ap);
239
240 elapsed_time = timer.ElapsedTime();
241
242 double normreducedb, normreduceda;
243 bp->NormInf(&normreducedb);
244 normreduceda = Ap->NormInf();
245 if (verbose)
246 cout << "Inf norm of Reduced Matrix = " << normreduceda << endl
247 << "Inf norm of Reduced RHS = " << normreducedb << endl;
248
249 double residual;
250 BiCGSTAB(*Ap, *xp, *bp, ILUK, Maxiter, Tolerance, &residual, verbose);
251
252 elapsed_time = timer.ElapsedTime() - elapsed_time;
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;
259
260 SingletonFilter.ComputeFullSolution();
261
262 if (smallProblem)
263 cout << " Reduced LHS after sol = " << endl << *xp << endl
264 << " Full LHS after sol = " << endl << *x << endl
265 << " Full Exact LHS = " << endl << *xexact << endl;
266
267 Epetra_Vector resid(*x);
268
269 resid.Update(1.0, *x, -1.0, *xexact, 0.0); // resid = xcomp - xexact
270
271 resid.Norm2(&residual);
272 double normx, normxexact;
273 x->Norm2(&normx);
274 xexact->Norm2(&normxexact);
275
276 if (verbose)
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;
280
281 if (verbose1 && residual>1.0e-5) {
282 if (verbose)
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"
285 << endl;
286 Epetra_Vector bdiff(*b);
287 assert(A->Multiply(false, resid, bdiff)==0);
288 assert(bdiff.Norm2(&residual)==0);
289 if (verbose)
290 cout << "2-norm of A times difference between computed and exact solution = " << residual << endl;
291
292 }
293 if (verbose)
294 cout << "********************************************************" << endl
295 << " Solving again with 2*Ax=2*b" << endl
296 << "********************************************************" << endl;
297
298 A->Scale(1.0); // A = 2*A
299 b->Scale(1.0); // b = 2*b
300 x->PutScalar(0.0);
301 b->NormInf(&normb);
302 norma = A->NormInf();
303 if (verbose)
304 cout << "Inf norm of Original Matrix = " << norma << endl
305 << "Inf norm of Original RHS = " << normb << endl;
306 double updateReducedProblemTime = ReductionTimer.ElapsedTime();
307 SingletonFilter.UpdateReducedProblem(&FullProblem);
308 Comm.Barrier();
309 updateReducedProblemTime = ReductionTimer.ElapsedTime() - updateReducedProblemTime;
310 if (verbose)
311 cout << "\n\n****************************************************" << endl
312 << " Update Reduced Problem time (sec) = " << updateReducedProblemTime<< endl
313 << "****************************************************" << endl;
314 Statistics(SingletonFilter);
315
316 if (LevelFill>-1) {
317
318 Epetra_Flops fact_counter;
319
320 elapsed_time = timer.ElapsedTime();
321
322 int initerr = ILUK->InitValues(*Ap);
323 if (initerr!=0) {
324 cout << endl << Comm << endl << " InitValues error = " << initerr;
325 if (initerr==1) cout << " Zero diagonal found, warning error only";
326 cout << endl << endl;
327 }
328 assert(ILUK->Factor()==0);
329 elapsed_time = timer.ElapsedTime() - elapsed_time;
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;
335 double Condest;
336 ILUK->Condest(false, Condest);
337
338 if (verbose) cout << "Condition number estimate for this preconditioner = " << Condest << endl;
339 }
340 bp->NormInf(&normreducedb);
341 normreduceda = Ap->NormInf();
342 if (verbose)
343 cout << "Inf norm of Reduced Matrix = " << normreduceda << endl
344 << "Inf norm of Reduced RHS = " << normreducedb << endl;
345
346 BiCGSTAB(*Ap, *xp, *bp, ILUK, Maxiter, Tolerance, &residual, verbose);
347
348 elapsed_time = timer.ElapsedTime() - elapsed_time;
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;
355
356 SingletonFilter.ComputeFullSolution();
357
358 if (smallProblem)
359 cout << " Reduced LHS after sol = " << endl << *xp << endl
360 << " Full LHS after sol = " << endl << *x << endl
361 << " Full Exact LHS = " << endl << *xexact << endl;
362
363 resid.Update(1.0, *x, -1.0, *xexact, 0.0); // resid = xcomp - xexact
364
365 resid.Norm2(&residual);
366 x->Norm2(&normx);
367 xexact->Norm2(&normxexact);
368
369 if (verbose)
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;
373
374 if (verbose1 && residual>1.0e-5) {
375 if (verbose)
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"
378 << endl;
379 Epetra_Vector bdiff(*b);
380 assert(A->Multiply(false, resid, bdiff)==0);
381 assert(bdiff.Norm2(&residual)==0);
382 if (verbose)
383 cout << "2-norm of A times difference between computed and exact solution = " << residual << endl;
384
385 }
386
387
388
389 if (ILUK!=0) delete ILUK;
390 if (IlukGraph!=0) delete IlukGraph;
391
392#ifdef EPETRA_MPI
393 MPI_Finalize() ;
394#endif
395
396return 0 ;
397}
399 Epetra_Vector &x,
400 Epetra_Vector &b,
401 Ifpack_CrsRiluk *M,
402 int Maxiter,
403 double Tolerance,
404 double *residual, bool verbose) {
405
406 // Allocate vectors needed for iterations
407 Epetra_Vector phat(x.Map()); phat.SetFlopCounter(x);
408 Epetra_Vector p(x.Map()); p.SetFlopCounter(x);
409 Epetra_Vector shat(x.Map()); shat.SetFlopCounter(x);
410 Epetra_Vector s(x.Map()); s.SetFlopCounter(x);
411 Epetra_Vector r(b.Map()); r.SetFlopCounter(x);
412 Epetra_Vector rtilde(x.Map()); rtilde.PutScalar(1.0); rtilde.SetFlopCounter(x);
413 Epetra_Vector v(x.Map()); v.SetFlopCounter(x);
414
415
416 A.Multiply(false, x, r); // r = A*x
417
418 r.Update(1.0, b, -1.0); // r = b - A*x
419
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;
423 r.Norm2(&r_norm);
424 b.Norm2(&b_norm);
425 scaled_r_norm = r_norm/b_norm;
426 r.Dot(rtilde,&rhon);
427
428 if (verbose) cout << "Initial residual = " << r_norm
429 << " Scaled residual = " << scaled_r_norm << endl;
430
431
432 for (int i=0; i<Maxiter; i++) { // Main iteration loop
433
434 double beta = (rhon/rhonm1) * (alpha/omega);
435 rhonm1 = rhon;
436
437 /* p = r + beta*(p - omega*v) */
438 /* phat = M^-1 p */
439 /* v = A phat */
440
441 double dtemp = - beta*omega;
442
443 p.Update(1.0, r, dtemp, v, beta);
444 if (M==0)
445 phat.Scale(1.0, p);
446 else
447 M->Solve(false, p, phat);
448 A.Multiply(false, phat, v);
449
450
451 rtilde.Dot(v,&sigma);
452 alpha = rhon/sigma;
453
454 /* s = r - alpha*v */
455 /* shat = M^-1 s */
456 /* r = A shat (r is a tmp here for t ) */
457
458 s.Update(-alpha, v, 1.0, r, 0.0);
459 if (M==0)
460 shat.Scale(1.0, s);
461 else
462 M->Solve(false, s, shat);
463 A.Multiply(false, shat, r);
464
465 r.Dot(s, &omega_num);
466 r.Dot(r, &omega_den);
467 omega = omega_num / omega_den;
468
469 /* x = x + alpha*phat + omega*shat */
470 /* r = s - omega*r */
471
472 x.Update(alpha, phat, omega, shat, 1.0);
473 r.Update(1.0, s, -omega);
474
475 r.Norm2(&r_norm);
476 scaled_r_norm = r_norm/b_norm;
477 r.Dot(rtilde,&rhon);
478
479 if (verbose) cout << "Iter "<< i << " residual = " << r_norm
480 << " Scaled residual = " << scaled_r_norm << endl;
481
482 if (r_norm < Tolerance) break;
483 }
484 return;
485}
486//==============================================================================
488
489 // Create local variables of some of the stats we need because we don't know if the
490 // method calls are collective or not for the particular implementation of the Row Matrix
491 int fmaxentries;
492 int maxentries = filter.FullMatrix()->MaxNumEntries();
493 filter.FullMatrix()->RowMatrixRowMap().Comm().SumAll(&maxentries, &fmaxentries, 1);
494 int fnrows = filter.FullMatrix()->NumGlobalRows();
495 int fnnz = filter.FullMatrix()->NumGlobalNonzeros();
496 if (filter.FullMatrix()->RowMatrixRowMap().Comm().MyPID()!=0) return(0);
497
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
503 << " Dimension = " << filter.ReducedMatrix()->NumGlobalRows() << endl
504 << " Number of nonzeros = " << filter.ReducedMatrix()->NumGlobalNonzeros() << endl
505 << " Maximum Number of Row Entries = " << filter.ReducedMatrix()->GlobalMaxNumEntries() << 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
511 << " row singletons) = " << filter.NumColSingletons() << endl << 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;
515
516 return(0);
517
518}
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.
Epetra_CrsSingletonFilter: A class for explicitly eliminating matrix rows and columns.
int ComputeFullSolution()
Compute a solution for the full problem using the solution of the reduced problem,...
int Analyze(Epetra_RowMatrix *FullMatrix)
Analyze the input matrix, removing row/column pairs that have singletons.
int ConstructReducedProblem(Epetra_LinearProblem *Problem)
Return a reduced linear problem based on results of Analyze().
Epetra_LinearProblem * ReducedProblem() const
Returns pointer to the derived reduced Epetra_LinearProblem.
Epetra_CrsMatrix * ReducedMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
int NumRowSingletons() const
Return number of rows that contain a single entry, returns -1 if Analysis not performed yet.
bool SingletonsDetected() const
Returns true if singletons were detected in this matrix (must be called after Analyze() to be effecti...
double RatioOfDimensions() const
Returns ratio of reduced system to full system dimensions, returns -1.0 if reduced problem not constr...
int NumSingletons() const
Return total number of singletons detected, returns -1 if Analysis not performed yet.
int UpdateReducedProblem(Epetra_LinearProblem *Problem)
Update a reduced linear problem using new values.
Epetra_RowMatrix * FullMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
int NumColSingletons() const
Return number of columns that contain a single entry that are not associated with singleton row,...
double RatioOfNonzeros() const
Returns ratio of reduced system to full system nonzero count, returns -1.0 if reduced problem not con...
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
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.
Definition: Epetra_Map.h:119
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.
Definition: Epetra_Time.h:75
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)