Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
58int 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
64void GenerateHilbert(double *A, int LDA, int N);
65
66bool Residual( int N, int NRHS, double * A, int LDA,
67 double * X, int LDX, double * B, int LDB, double * resid);
68
69
70int 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
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
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*/
396return ierr ;
397}
398
399int 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
534bool 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}
@ View
@ Copy
std::string Epetra_Version()
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
float NRM2(const int N, const float *X, const int INCX=1) const
Epetra_BLAS norm function (SNRM2).
Definition: Epetra_BLAS.cpp:81
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)
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
double Flops() const
Returns the number of floating point operations with this multi-vector.
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
Epetra_MpiComm: The Epetra MPI Communication Class.
int MyPID() const
Return my process ID.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
double * A() const
Returns pointer to the this matrix.
int LDA() const
Returns the leading dimension of the this matrix.
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
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.
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
int M() const
Returns row dimension of system.
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
int LDB() const
Returns the leading dimension of the RHS.
int LDX() const
Returns the leading dimension of the solution.
bool SolutionRefined()
Returns true if the current set of vectors has been refined.
void SolveToRefinedSolution(bool Flag)
Causes all solves to compute solution to best ability using iterative refinement.
double * X() const
Returns pointer to current solution.
bool SolutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
int N() const
Returns column dimension of system.
bool A_Equilibrated()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
double * FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
int NRHS() const
Returns the number of current right hand sides and solution vectors.
double ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
double RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
int LDA() const
Returns the leading dimension of the this matrix.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
bool B_Equilibrated()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
double * BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
bool Solved()
Returns true if the current set of vectors has been solved.
void FactorWithEquilibration(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor.
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
Epetra_SerialSpdDenseSolver: A class for constructing and using symmetric positive definite dense mat...
int Invert(void)
Inverts the this matrix.
int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
int Factor(void)
Computes the in-place Cholesky factorization of the matrix using the LAPACK routine DPOTRF.
int SetMatrix(Epetra_SerialSymDenseMatrix &A_in)
Sets the pointers for coefficient matrix; special version for symmetric matrices.
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 ReciprocalConditionEstimate(double &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
double AMAX()
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
Epetra_SerialSymDenseMatrix: A class for constructing and using symmetric positive definite dense mat...
int Shape(int NumRowsCols)
Set dimensions of a Epetra_SerialSymDenseMatrix object; init values to zero.
bool Upper() const
Returns true if upper triangle of this matrix has and will be used.
double NormInf() const
Computes the Infinity-Norm of the this matrix.
void SetUpper()
Specify that the upper triangle of the this matrix should be used.
double NormOne() const
Computes the 1-Norm of the this matrix.
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
double ElapsedTime(void) const
Epetra_Time elapsed time function.
bool Residual(int N, int NRHS, double *A, int LDA, double *X, int LDX, double *B, int LDB, double *resid)
int main(int argc, char *argv[])
void GenerateHilbert(double *A, int LDA, int N)
int check(Epetra_SerialSpdDenseSolver &solver, double *A1, int LDA, int N1, int NRHS1, double OneNorm1, double *B1, int LDB1, double *X1, int LDX1, bool Upper, bool verbose)