Epetra Package Browser (Single Doxygen Collection)  Development
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"
50 #include "Epetra_LinearProblem.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);
69 int Statistics(const Epetra_CrsSingletonFilter & filter);
70 int main(int argc, char *argv[]) {
71 
72 #ifdef EPETRA_MPI
73  MPI_Init(&argc,&argv);
74  Epetra_MpiComm Comm (MPI_COMM_WORLD);
75 #else
76  Epetra_SerialComm Comm;
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;
110  Epetra_CrsMatrix * A;
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 
396 return 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 }
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
int GlobalMaxNumEntries() const
Returns the maximum number of nonzero entries across all rows on all processors.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
void Barrier() const
Epetra_SerialComm Barrier function.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
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.
void BiCGSTAB(Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, Ifpack_CrsRiluk *M, int Maxiter, double Tolerance, double *residual, bool verbose)
Epetra_LinearProblem * ReducedProblem() const
Returns pointer to the derived reduced Epetra_LinearProblem.
int UpdateReducedProblem(Epetra_LinearProblem *Problem)
Update a reduced linear problem using new values.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
int Analyze(Epetra_RowMatrix *FullMatrix)
Analyze the input matrix, removing row/column pairs that have singletons.
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int ComputeFullSolution()
Compute a solution for the full problem using the solution of the reduced problem, put in LHS of FullProblem().
double RatioOfDimensions() const
Returns ratio of reduced system to full system dimensions, returns -1.0 if reduced problem not constr...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
Epetra_MpiComm: The Epetra MPI Communication Class.
virtual int NumGlobalNonzeros() const =0
Returns the number of nonzero entries in the global matrix.
virtual int MyPID() const =0
Return my process ID.
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
int NumColSingletons() const
Return number of columns that contain a single entry that are not associated with singleton row...
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int Statistics(const Epetra_CrsSingletonFilter &filter)
Epetra_CrsMatrix * ReducedMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
Epetra_RowMatrix * FullMatrix() const
Returns pointer to Epetra_CrsMatrix from full problem.
double RatioOfNonzeros() const
Returns ratio of reduced system to full system nonzero count, returns -1.0 if reduced problem not con...
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
Epetra_CrsSingletonFilter: A class for explicitly eliminating matrix rows and columns.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
int NumSingletons() const
Return total number of singletons detected, returns -1 if Analysis not performed yet.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
Epetra_MultiVector * GetRHS() const
Get a pointer to the right-hand-side B.
int ConstructReducedProblem(Epetra_LinearProblem *Problem)
Return a reduced linear problem based on results of Analyze().
int MyPID() const
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int NumGlobalElements() const
Number of elements across all processors.
double NormInf() const
Returns the infinity norm of the global matrix.
Epetra_MultiVector * GetLHS() const
Get a pointer to the left-hand-side X.
Epetra_RowMatrix * GetMatrix() const
Get a pointer to the matrix A.
Epetra_LinearProblem: The Epetra Linear Problem Class.
int NumGlobalRows() const
Returns the number of global matrix rows.
int main(int argc, char *argv[])
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A <- ScalarConstant * A)...
bool SingletonsDetected() const
Returns true if singletons were detected in this matrix (must be called after Analyze() to be effecti...
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
virtual int NumGlobalRows() const =0
Returns the number of global matrix rows.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
int NumRowSingletons() const
Return number of rows that contain a single entry, returns -1 if Analysis not performed yet...