59 int N1,
int NRHS1,
double OneNorm1,
60 double * B1,
int LDB1,
61 double * X1,
int LDX1,
62 bool Upper,
bool verbose);
66bool Residual(
int N,
int NRHS,
double * A,
int LDA,
67 double * X,
int LDX,
double * B,
int LDB,
double * resid);
70int main(
int argc,
char *argv[])
72 int ierr = 0, i, j, k;
75 MPI_Init(&argc,&argv);
84 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
86 if(verbose && Comm.
MyPID()==0)
89 int rank = Comm.
MyPID();
96 if (verbose) std::cout << Comm << std::endl;
101 if (verbose && rank!=0) verbose =
false;
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];
111 double * B =
new double[N*NRHS];
112 double * B1 =
new double[N*NRHS];
123 for (
int kk=0; kk<2; kk++) {
124 for (i=1; i<=N; i++) {
127 for (j=1; j<=i; j++) OneNorm1 += 1.0/((
double) j);
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];
156 ierr =
check(solver, A1, LDA1, i, NRHS, OneNorm1, B1, LDB1, X1, LDX1, Upper, verbose);
160 if (verbose) std::cout <<
"Factorization failed due to bad conditioning. This is normal if SCOND is small."
179 double ScalarA = 2.0;
185 for (i=0; i<DM; i++) D[j][i] = (
double) (1+i+j*DM) ;
189 double NormInfD_ref = (double)(DM*(DN*(DN+1))/2);
190 double NormOneD_ref = NormInfD_ref;
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;
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;
229 if (verbose) std::cout <<
"\n\nComputing factor of an " << N <<
" x " << N <<
" SPD matrix...Please wait.\n\n" << std::endl;
233 A =
new double[LDA*N];
234 X =
new double[LDB*NRHS];
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));
259 ierr = BigSolver.
Factor();
260 if (ierr!=0 && verbose) std::cout <<
"Error in factorization = "<<ierr<< std::endl;
264 double FLOPS = counter.
Flops();
265 double MFLOPS = FLOPS/time/1000000.0;
266 if (verbose) std::cout <<
"MFLOPS for Factorization = " << MFLOPS << std::endl;
278 RHS.
Multiply(
'L', 1.0, OrigBigMatrix, LHS, 0.0);
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;
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;
297 FLOPS = BigSolver.
Flops();
298 MFLOPS = FLOPS/time/1000000.0;
299 if (verbose) std::cout <<
"MFLOPS for Solve (NRHS = " << NRHS <<
") = " << MFLOPS << std::endl;
301 double * resid =
new double[NRHS];
302 bool OK =
Residual(N, NRHS, A, LDA, BigSolver.
X(), BigSolver.
LDX(),
303 OrigRHS.
A(), OrigRHS.
LDA(), resid);
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;
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);}
324 B2.
Multiply(
'N',
'N', 1.0, OrigBigMatrix, X2, 0.0);
329 FLOPS = RHS_counter.
Flops();
330 MFLOPS = FLOPS/time/1000000.0;
331 if (verbose) std::cout <<
"MFLOPS to build single RHS = " << MFLOPS << std::endl;
337 ierr = BigSolver.
Solve();
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;
343 FLOPS = counter.
Flops();
344 MFLOPS = FLOPS/time/1000000.0;
345 if (verbose) std::cout <<
"MFLOPS to solve single RHS = " << MFLOPS << std::endl;
347 OK =
Residual(N, 1, A, LDA, BigSolver.
X(), BigSolver.
LDX(), OrigB2.
A(),
348 OrigB2.
LDA(), resid);
351 if (!OK) std::cout <<
"************* Residual do not meet tolerance *************" << std::endl;
352 std::cout <<
"Residual = "<< resid[0] << std::endl;
365 double * C1 =
new double[N*N];
377 for (j=0; j<N; j++) {
378 assert(C(i,j) == C1[i+j*N]);
379 assert(C(i,j) == C[j][i]);
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;
400 int N1,
int NRHS1,
double OneNorm1,
401 double * B1,
int LDB1,
402 double * X1,
int LDX1,
403 bool Upper,
bool verbose) {
410 if (verbose) std::cout <<
"\n\nNumber of Rows = " << M << std::endl<< std::endl;
414 if (verbose) std::cout <<
"\n\nNumber of Equations = " << N << std::endl<< std::endl;
417 int LDA = solver.
LDA();
418 if (verbose) std::cout <<
"\n\nLDA = " << LDA << std::endl<< std::endl;
421 int LDB = solver.
LDB();
422 if (verbose) std::cout <<
"\n\nLDB = " << LDB << std::endl<< std::endl;
425 int LDX = solver.
LDX();
426 if (verbose) std::cout <<
"\n\nLDX = " << LDX << std::endl<< std::endl;
429 int NRHS = solver.
NRHS();
430 if (verbose) std::cout <<
"\n\nNRHS = " << NRHS << std::endl<< std::endl;
433 assert(solver.
ANORM()==-1.0);
434 assert(solver.
RCOND()==-1.0);
436 assert(solver.
SCOND()==-1.0);
437 assert(solver.
AMAX()==-1.0);
454 int ierr = solver.
Factor();
458 if (ierr!=0)
return(ierr);
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;
470 ierr = solver.
Solve();
472 if (ierr!=0 && verbose) std::cout <<
"LAPACK rules suggest system should be equilibrated." << std::endl;
481 std::cout <<
"\n\nFERR[0] = "<< solver.
FERR()[0] << std::endl;
482 std::cout <<
"\n\nBERR[0] = "<< solver.
BERR()[0] << std::endl<< std::endl;
486 double * resid =
new double[NRHS];
487 OK =
Residual(N, NRHS, A1, LDA1, solver.
X(), solver.
LDX(), B1, LDB1, resid);
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;
508 assert(solver.
Solve()>-1);
512 OK =
Residual(N, NRHS, A1, LDA1, solver.
X(), solver.
LDX(), B1, LDB1, resid);
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;
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));
534bool Residual(
int N,
int NRHS,
double * A,
int LDA,
535 double * X,
int LDX,
double * B,
int LDB,
double * resid) {
539 Blas.
GEMM(Transa,
'N', N, NRHS, N, -1.0, A, LDA,
540 X, LDX, 1.0, B, LDB);
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;
std::string Epetra_Version()
Epetra_BLAS: The Epetra BLAS Wrapper Class.
float NRM2(const int N, const float *X, const int INCX=1) const
Epetra_BLAS norm function (SNRM2).
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.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
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.
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)