Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/SerialSpdDense/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_Map.h"
44 #include "Epetra_Time.h"
49 #ifdef EPETRA_MPI
50 #include "Epetra_MpiComm.h"
51 #include <mpi.h>
52 #endif
53 #include "Epetra_SerialComm.h"
54 #include "Epetra_Version.h"
55 
56 // prototypes
57 
58 int check(Epetra_SerialSpdDenseSolver & solver, double * A1, int LDA,
59  int N1, int NRHS1, double OneNorm1,
60  double * B1, int LDB1,
61  double * X1, int LDX1,
62  bool Upper, bool verbose);
63 
64 void GenerateHilbert(double *A, int LDA, int N);
65 
66 bool Residual( int N, int NRHS, double * A, int LDA,
67  double * X, int LDX, double * B, int LDB, double * resid);
68 
69 
70 int main(int argc, char *argv[])
71 {
72  int ierr = 0, i, j, k;
73 
74 #ifdef EPETRA_MPI
75  MPI_Init(&argc,&argv);
76  Epetra_MpiComm Comm( MPI_COMM_WORLD );
77 #else
78  Epetra_SerialComm Comm;
79 #endif
80 
81  bool verbose = false;
82 
83  // Check if we should print results to standard out
84  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
85 
86  if(verbose && Comm.MyPID()==0)
87  std::cout << Epetra_Version() << std::endl << std::endl;
88 
89  int rank = Comm.MyPID();
90  // char tmp;
91  // if (rank==0) std::cout << "Press any key to continue..."<< std::endl;
92  // if (rank==0) cin >> tmp;
93  // Comm.Barrier();
94 
95  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
96  if (verbose) std::cout << Comm << std::endl;
97 
98  // bool verbose1 = verbose;
99 
100  // Redefine verbose to only print on PE 0
101  if (verbose && rank!=0) verbose = false;
102 
103  int N = 20;
104  int NRHS = 4;
105  double * A = new double[N*N];
106  double * A1 = new double[N*N];
107  double * X = new double[(N+1)*NRHS];
108  double * X1 = new double[(N+1)*NRHS];
109  int LDX = N+1;
110  int LDX1 = N+1;
111  double * B = new double[N*NRHS];
112  double * B1 = new double[N*NRHS];
113  int LDB = N;
114  int LDB1 = N;
115 
116  int LDA = N;
117  int LDA1 = LDA;
118  double OneNorm1;
119  bool Upper = false;
120 
123  for (int kk=0; kk<2; kk++) {
124  for (i=1; i<=N; i++) {
125  GenerateHilbert(A, LDA, i);
126  OneNorm1 = 0.0;
127  for (j=1; j<=i; j++) OneNorm1 += 1.0/((double) j); // 1-Norm = 1 + 1/2 + ...+1/n
128 
129  if (kk==0) {
130  Matrix = new Epetra_SerialSymDenseMatrix(View, A, LDA, i);
131  LDA1 = LDA;
132  }
133  else {
134  Matrix = new Epetra_SerialSymDenseMatrix(Copy, A, LDA, i);
135  LDA1 = i;
136  }
137  GenerateHilbert(A1, LDA1, i);
138 
139  if (kk==1) {
140  solver.FactorWithEquilibration(true);
141  Matrix->SetUpper();
142  Upper = true;
143  solver.SolveToRefinedSolution(false);
144  }
145 
146  for (k=0; k<NRHS; k++)
147  for (j=0; j<i; j++) {
148  B[j+k*LDB] = 1.0/((double) (k+3)*(j+3));
149  B1[j+k*LDB1] = B[j+k*LDB1];
150  }
151  Epetra_SerialDenseMatrix Epetra_B(View, B, LDB, i, NRHS);
152  Epetra_SerialDenseMatrix Epetra_X(View, X, LDX, i, NRHS);
153  solver.SetMatrix(*Matrix);
154  solver.SetVectors(Epetra_X, Epetra_B);
155 
156  ierr = check(solver, A1, LDA1, i, NRHS, OneNorm1, B1, LDB1, X1, LDX1, Upper, verbose);
157  assert (ierr>-1);
158  delete Matrix;
159  if (ierr!=0) {
160  if (verbose) std::cout << "Factorization failed due to bad conditioning. This is normal if SCOND is small."
161  << std::endl;
162  break;
163  }
164  }
165  }
166 
167  delete [] A;
168  delete [] A1;
169  delete [] X;
170  delete [] X1;
171  delete [] B;
172  delete [] B1;
173 
175  // Now test norms and scaling functions
177 
179  double ScalarA = 2.0;
180 
181  int DM = 10;
182  int DN = 10;
183  D.Shape(DM);
184  for (j=0; j<DN; j++)
185  for (i=0; i<DM; i++) D[j][i] = (double) (1+i+j*DM) ;
186 
187  //std::cout << D << std::endl;
188 
189  double NormInfD_ref = (double)(DM*(DN*(DN+1))/2);
190  double NormOneD_ref = NormInfD_ref;
191 
192  double NormInfD = D.NormInf();
193  double NormOneD = D.NormOne();
194 
195  if (verbose) {
196  std::cout << " *** Before scaling *** " << std::endl
197  << " Computed one-norm of test matrix = " << NormOneD << std::endl
198  << " Expected one-norm = " << NormOneD_ref << std::endl
199  << " Computed inf-norm of test matrix = " << NormInfD << std::endl
200  << " Expected inf-norm = " << NormInfD_ref << std::endl;
201  }
202  D.Scale(ScalarA); // Scale entire D matrix by this value
203 
204  //std::cout << D << std::endl;
205 
206  NormInfD = D.NormInf();
207  NormOneD = D.NormOne();
208  if (verbose) {
209  std::cout << " *** After scaling *** " << std::endl
210  << " Computed one-norm of test matrix = " << NormOneD << std::endl
211  << " Expected one-norm = " << NormOneD_ref*ScalarA << std::endl
212  << " Computed inf-norm of test matrix = " << NormInfD << std::endl
213  << " Expected inf-norm = " << NormInfD_ref*ScalarA << std::endl;
214  }
215 
216 
217 
219  // Now test for larger system, both correctness and performance.
221 
222 
223  N = 2000;
224  NRHS = 5;
225  LDA = N;
226  LDB = N;
227  LDX = N;
228 
229  if (verbose) std::cout << "\n\nComputing factor of an " << N << " x " << N << " SPD matrix...Please wait.\n\n" << std::endl;
230 
231  // Define A and X
232 
233  A = new double[LDA*N];
234  X = new double[LDB*NRHS];
235 
236  for (j=0; j<N; j++) {
237  for (k=0; k<NRHS; k++) X[j+k*LDX] = 1.0/((double) (j+5+k));
238  for (i=0; i<N; i++) {
239  if (i==j) A[i+j*LDA] = 100.0 + i;
240  else A[i+j*LDA] = -1.0/((double) (i+10)*(j+10));
241  }
242  }
243 
244  // Define Epetra_SerialDenseMatrix object
245 
246  Epetra_SerialSymDenseMatrix BigMatrix(Copy, A, LDA, N);
247  Epetra_SerialSymDenseMatrix OrigBigMatrix(View, A, LDA, N);
248 
249  Epetra_SerialSpdDenseSolver BigSolver;
250  BigSolver.FactorWithEquilibration(true);
251  BigSolver.SetMatrix(BigMatrix);
252 
253  // Time factorization
254 
255  Epetra_Flops counter;
256  BigSolver.SetFlopCounter(counter);
257  Epetra_Time Timer(Comm);
258  double tstart = Timer.ElapsedTime();
259  ierr = BigSolver.Factor();
260  if (ierr!=0 && verbose) std::cout << "Error in factorization = "<<ierr<< std::endl;
261  assert(ierr==0);
262  double time = Timer.ElapsedTime() - tstart;
263 
264  double FLOPS = counter.Flops();
265  double MFLOPS = FLOPS/time/1000000.0;
266  if (verbose) std::cout << "MFLOPS for Factorization = " << MFLOPS << std::endl;
267 
268  // Define Left hand side and right hand side
269  Epetra_SerialDenseMatrix LHS(View, X, LDX, N, NRHS);
271  RHS.Shape(N,NRHS); // Allocate RHS
272 
273  // Compute RHS from A and X
274 
275  Epetra_Flops RHS_counter;
276  RHS.SetFlopCounter(RHS_counter);
277  tstart = Timer.ElapsedTime();
278  RHS.Multiply('L', 1.0, OrigBigMatrix, LHS, 0.0); // Symmetric Matrix-multiply
279  time = Timer.ElapsedTime() - tstart;
280 
281  Epetra_SerialDenseMatrix OrigRHS = RHS;
282 
283  FLOPS = RHS_counter.Flops();
284  MFLOPS = FLOPS/time/1000000.0;
285  if (verbose) std::cout << "MFLOPS to build RHS (NRHS = " << NRHS <<") = " << MFLOPS << std::endl;
286 
287  // Set LHS and RHS and solve
288  BigSolver.SetVectors(LHS, RHS);
289 
290  tstart = Timer.ElapsedTime();
291  ierr = BigSolver.Solve();
292  if (ierr==1 && verbose) std::cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << std::endl;
293  else if (ierr!=0 && verbose) std::cout << "Error in solve = "<<ierr<< std::endl;
294  assert(ierr>=0);
295  time = Timer.ElapsedTime() - tstart;
296 
297  FLOPS = BigSolver.Flops();
298  MFLOPS = FLOPS/time/1000000.0;
299  if (verbose) std::cout << "MFLOPS for Solve (NRHS = " << NRHS <<") = " << MFLOPS << std::endl;
300 
301  double * resid = new double[NRHS];
302  bool OK = Residual(N, NRHS, A, LDA, BigSolver.X(), BigSolver.LDX(),
303  OrigRHS.A(), OrigRHS.LDA(), resid);
304 
305  if (verbose) {
306  if (!OK) std::cout << "************* Residual do not meet tolerance *************" << std::endl;
307  for (i=0; i<NRHS; i++)
308  std::cout << "Residual[" << i <<"] = "<< resid[i] << std::endl;
309  std::cout << std::endl;
310  }
311 
312  // Solve again using the Epetra_SerialDenseVector class for LHS and RHS
313 
316  X2.Size(BigMatrix.N());
317  B2.Size(BigMatrix.M());
318  int length = BigMatrix.N();
319  {for (int kk=0; kk<length; kk++) X2[kk] = ((double ) kk)/ ((double) length);} // Define entries of X2
320 
321  RHS_counter.ResetFlops();
322  B2.SetFlopCounter(RHS_counter);
323  tstart = Timer.ElapsedTime();
324  B2.Multiply('N', 'N', 1.0, OrigBigMatrix, X2, 0.0); // Define B2 = A*X2
325  time = Timer.ElapsedTime() - tstart;
326 
327  Epetra_SerialDenseVector OrigB2 = B2;
328 
329  FLOPS = RHS_counter.Flops();
330  MFLOPS = FLOPS/time/1000000.0;
331  if (verbose) std::cout << "MFLOPS to build single RHS = " << MFLOPS << std::endl;
332 
333  // Set LHS and RHS and solve
334  BigSolver.SetVectors(X2, B2);
335 
336  tstart = Timer.ElapsedTime();
337  ierr = BigSolver.Solve();
338  time = Timer.ElapsedTime() - tstart;
339  if (ierr==1 && verbose) std::cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << std::endl;
340  else if (ierr!=0 && verbose) std::cout << "Error in solve = "<<ierr<< std::endl;
341  assert(ierr>=0);
342 
343  FLOPS = counter.Flops();
344  MFLOPS = FLOPS/time/1000000.0;
345  if (verbose) std::cout << "MFLOPS to solve single RHS = " << MFLOPS << std::endl;
346 
347  OK = Residual(N, 1, A, LDA, BigSolver.X(), BigSolver.LDX(), OrigB2.A(),
348  OrigB2.LDA(), resid);
349 
350  if (verbose) {
351  if (!OK) std::cout << "************* Residual do not meet tolerance *************" << std::endl;
352  std::cout << "Residual = "<< resid[0] << std::endl;
353  }
354  delete [] resid;
355  delete [] A;
356  delete [] X;
357 
359  // Now test default constructor and index operators
361 
362  N = 5;
363  Epetra_SerialSymDenseMatrix C; // Implicit call to default constructor, should not need to call destructor
364  C.Shape(5); // Make it 5 by 5
365  double * C1 = new double[N*N];
366  GenerateHilbert(C1, N, N); // Generate Hilber matrix
367 
368  C1[1+2*N] = 1000.0; // Make matrix nonsymmetric
369 
370  // Fill values of C with Hilbert values
371  for (i=0; i<N; i++)
372  for (j=0; j<N; j++)
373  C(i,j) = C1[i+j*N];
374 
375  // Test if values are correctly written and read
376  for (i=0; i<N; i++)
377  for (j=0; j<N; j++) {
378  assert(C(i,j) == C1[i+j*N]);
379  assert(C(i,j) == C[j][i]);
380  }
381 
382  if (verbose)
383  std::cout << "Default constructor and index operator check OK. Values of Hilbert matrix = "
384  << std::endl << C << std::endl
385  << "Values should be 1/(i+j+1), except value (1,2) should be 1000" << std::endl;
386 
387  delete [] C1;
388 
389 
390 #ifdef EPETRA_MPI
391  MPI_Finalize() ;
392 #endif
393 
394 /* end main
395 */
396 return ierr ;
397 }
398 
399 int check(Epetra_SerialSpdDenseSolver &solver, double * A1, int LDA1,
400  int N1, int NRHS1, double OneNorm1,
401  double * B1, int LDB1,
402  double * X1, int LDX1,
403  bool Upper, bool verbose) {
404  (void)OneNorm1;
405  int i;
406  bool OK;
407  // Test query functions
408 
409  int M= solver.M();
410  if (verbose) std::cout << "\n\nNumber of Rows = " << M << std::endl<< std::endl;
411  assert(M==N1);
412 
413  int N= solver.N();
414  if (verbose) std::cout << "\n\nNumber of Equations = " << N << std::endl<< std::endl;
415  assert(N==N1);
416 
417  int LDA = solver.LDA();
418  if (verbose) std::cout << "\n\nLDA = " << LDA << std::endl<< std::endl;
419  assert(LDA==LDA1);
420 
421  int LDB = solver.LDB();
422  if (verbose) std::cout << "\n\nLDB = " << LDB << std::endl<< std::endl;
423  assert(LDB==LDB1);
424 
425  int LDX = solver.LDX();
426  if (verbose) std::cout << "\n\nLDX = " << LDX << std::endl<< std::endl;
427  assert(LDX==LDX1);
428 
429  int NRHS = solver.NRHS();
430  if (verbose) std::cout << "\n\nNRHS = " << NRHS << std::endl<< std::endl;
431  assert(NRHS==NRHS1);
432 
433  assert(solver.ANORM()==-1.0);
434  assert(solver.RCOND()==-1.0);
435  if (!solver.A_Equilibrated() && !solver.B_Equilibrated()) {
436  assert(solver.SCOND()==-1.0);
437  assert(solver.AMAX()==-1.0);
438  }
439 
440 
441  // Other binary tests
442 
443  assert(!solver.Factored());
444  assert(solver.SymMatrix()->Upper()==Upper);
445  assert(!solver.SolutionErrorsEstimated());
446  assert(!solver.Inverted());
447  assert(!solver.ReciprocalConditionEstimated());
448  assert(!solver.Solved());
449 
450  assert(!solver.SolutionRefined());
451 
452  //std::cout << "Matrix before factorization " << std::endl << *solver.SymMatrix() << std::endl << std::endl;
453 
454  int ierr = solver.Factor();
455  //std::cout << "Matrix after factorization " << std::endl << *solver.SymMatrix() << std::endl << std::endl;
456  //std::cout << "Factor after factorization " << std::endl << *solver.SymFactoredMatrix() << std::endl << std::endl;
457  assert(ierr>-1);
458  if (ierr!=0) return(ierr); // Factorization failed due to poor conditioning.
459  double rcond;
460  ierr = solver.ReciprocalConditionEstimate(rcond);
461  assert(ierr==0);
462  if (verbose) {
463 
464  double rcond1 = 1.0/std::exp(3.5*((double)N));
465  if (N==1) rcond1 = 1.0;
466  std::cout << "\n\nSCOND = "<< rcond << " should be approx = "
467  << rcond1 << std::endl << std::endl;
468  }
469 
470  ierr = solver.Solve();
471  assert(ierr>-1);
472  if (ierr!=0 && verbose) std::cout << "LAPACK rules suggest system should be equilibrated." << std::endl;
473 
474  assert(solver.Factored());
475  assert(solver.SymMatrix()->Upper()==Upper);
476  assert(solver.ReciprocalConditionEstimated());
477  assert(solver.Solved());
478 
479  if (solver.SolutionErrorsEstimated()) {
480  if (verbose) {
481  std::cout << "\n\nFERR[0] = "<< solver.FERR()[0] << std::endl;
482  std::cout << "\n\nBERR[0] = "<< solver.BERR()[0] << std::endl<< std::endl;
483  }
484  }
485 
486  double * resid = new double[NRHS];
487  OK = Residual(N, NRHS, A1, LDA1, solver.X(), solver.LDX(), B1, LDB1, resid);
488  if (verbose) {
489  if (!OK) std::cout << "************* Residual do not meet tolerance *************" << std::endl;
490  std::cout << "\n\nResiduals using factorization to solve" << std::endl;
491  for (i=0; i<NRHS; i++)
492  std::cout << "Residual[" << i <<"] = "<< resid[i] << std::endl;
493  std::cout << std::endl;
494  }
495 
496 
497  ierr = solver.Invert();
498  assert(ierr>-1);
499 
500  assert(solver.Inverted());
501  assert(!solver.Factored());
502 
503  Epetra_SerialDenseMatrix RHS1(Copy, B1, LDB1, N, NRHS);
504  Epetra_SerialDenseMatrix LHS1(Copy, X1, LDX1, N, NRHS);
505  assert(solver.SetVectors(LHS1, RHS1)==0);
506  assert(!solver.Solved());
507 
508  assert(solver.Solve()>-1);
509 
510 
511 
512  OK = Residual(N, NRHS, A1, LDA1, solver.X(), solver.LDX(), B1, LDB1, resid);
513 
514  if (verbose) {
515  if (!OK) std::cout << "************* Residual do not meet tolerance *************" << std::endl;
516  std::cout << "Residuals using inverse to solve" << std::endl;
517  for (i=0; i<NRHS; i++)
518  std::cout << "Residual[" << i <<"] = "<< resid[i] << std::endl;
519  std::cout << std::endl;
520  }
521  delete [] resid;
522 
523 
524  return(0);
525 }
526 
527  void GenerateHilbert(double *A, int LDA, int N) {
528  for (int j=0; j<N; j++)
529  for (int i=0; i<N; i++)
530  A[i+j*LDA] = 1.0/((double)(i+j+1));
531  return;
532  }
533 
534 bool Residual( int N, int NRHS, double * A, int LDA,
535  double * X, int LDX, double * B, int LDB, double * resid) {
536 
537  Epetra_BLAS Blas;
538  char Transa = 'N';
539  Blas.GEMM(Transa, 'N', N, NRHS, N, -1.0, A, LDA,
540  X, LDX, 1.0, B, LDB);
541  bool OK = true;
542  for (int i=0; i<NRHS; i++) {
543  resid[i] = Blas.NRM2(N, B+i*LDB);
544  if (resid[i]>1.0E-7) OK = false;
545  }
546 
547  return(OK);
548 }
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS matrix-matrix multiply function (SGEMM)
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
int Invert(void)
Inverts the this matrix.
float NRM2(const int N, const float *X, const int INCX=1) const
Epetra_BLAS norm function (SNRM2).
Definition: Epetra_BLAS.cpp:81
bool B_Equilibrated()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
double ElapsedTime(void) const
Epetra_Time elapsed time function.
bool SolutionRefined()
Returns true if the current set of vectors has been refined.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
double * BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
int SetMatrix(Epetra_SerialSymDenseMatrix &A_in)
Sets the pointers for coefficient matrix; special version for symmetric matrices. ...
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
int M() const
Returns row dimension of system.
bool SolutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
std::string Epetra_Version()
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
double ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
Epetra_SerialSymDenseMatrix: A class for constructing and using symmetric positive definite dense mat...
double NormInf() const
Computes the Infinity-Norm of the this matrix.
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
double RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
double * A() const
Returns pointer to the this matrix.
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)
Epetra_SerialSpdDenseSolver: A class for constructing and using symmetric positive definite dense mat...
int N() const
Returns column dimension of system.
int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
int ReciprocalConditionEstimate(double &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
double SCOND()
Ratio of smallest to largest equilibration scale factors for the this matrix (returns -1 if not yet c...
Epetra_SerialSymDenseMatrix * SymMatrix() const
Returns pointer to current matrix.
int LDB() const
Returns the leading dimension of the RHS.
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
double * X() const
Returns pointer to current solution.
int LDA() const
Returns the leading dimension of the this matrix.
double NormOne() const
Computes the 1-Norm of the this matrix.
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
void SetUpper()
Specify that the upper triangle of the this matrix should be used.
bool Upper() const
Returns true if upper triangle of this matrix has and will be used.
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
void SolveToRefinedSolution(bool Flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
int LDA() const
Returns the leading dimension of the this matrix.
int LDX() const
Returns the leading dimension of the solution.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_SerialDenseMatrix &A, const Epetra_SerialDenseMatrix &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
bool A_Equilibrated()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
double Flops() const
Returns the number of floating point operations with this multi-vector.
void FactorWithEquilibration(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor...
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
double * FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
bool Solved()
Returns true if the current set of vectors has been solved.
int main(int argc, char *argv[])
int N() const
Returns column dimension of system.
int NRHS() const
Returns the number of current right hand sides and solution vectors.
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
int Factor(void)
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF...
int M() const
Returns row dimension of system.
void GenerateHilbert(double *A, int LDA, int N)
double AMAX()
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
int Shape(int NumRowsCols)
Set dimensions of a Epetra_SerialSymDenseMatrix object; init values to zero.