test/SerialSpdDense/cxx_main.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #include "Epetra_Map.h"
00030 #include "Epetra_Time.h"
00031 #include "Epetra_SerialSymDenseMatrix.h" 
00032 #include "Epetra_SerialDenseMatrix.h" 
00033 #include "Epetra_SerialDenseVector.h"
00034 #include "Epetra_SerialSpdDenseSolver.h"
00035 #ifdef EPETRA_MPI
00036 #include "Epetra_MpiComm.h"
00037 #include <mpi.h>
00038 #endif
00039 #include "Epetra_SerialComm.h"
00040 #include "Epetra_Version.h"
00041 
00042 // prototypes
00043 
00044 int check(Epetra_SerialSpdDenseSolver & solver, double * A1, int LDA,
00045     int N1, int NRHS1, double OneNorm1, 
00046     double * B1, int LDB1, 
00047     double * X1, int LDX1,
00048     bool Upper, bool verbose);
00049 
00050 void GenerateHilbert(double *A, int LDA, int N);
00051 
00052 bool Residual( int N, int NRHS, double * A, int LDA,
00053          double * X, int LDX, double * B, int LDB, double * resid);
00054 
00055  
00056 int main(int argc, char *argv[])
00057 {
00058   int ierr = 0, i, j, k;
00059 
00060 #ifdef EPETRA_MPI
00061   MPI_Init(&argc,&argv);
00062   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00063 #else
00064   Epetra_SerialComm Comm;
00065 #endif
00066 
00067   bool verbose = false;
00068 
00069   // Check if we should print results to standard out
00070   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00071 
00072   if(verbose && Comm.MyPID()==0)
00073     cout << Epetra_Version() << endl << endl;
00074 
00075   int rank = Comm.MyPID();
00076   //  char tmp;
00077   //  if (rank==0) cout << "Press any key to continue..."<< endl;
00078   //  if (rank==0) cin >> tmp;
00079   //  Comm.Barrier();
00080 
00081   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00082   if (verbose) cout << Comm <<endl;
00083 
00084   //  bool verbose1 = verbose;
00085 
00086   // Redefine verbose to only print on PE 0
00087   if (verbose && rank!=0) verbose = false;
00088 
00089   int N = 20;
00090   int NRHS = 4;
00091   double * A = new double[N*N];
00092   double * A1 = new double[N*N];
00093   double * X = new double[(N+1)*NRHS];
00094   double * X1 = new double[(N+1)*NRHS];
00095   int LDX = N+1;
00096   int LDX1 = N+1;
00097   double * B = new double[N*NRHS];
00098   double * B1 = new double[N*NRHS];
00099   int LDB = N;
00100   int LDB1 = N;
00101 
00102   int LDA = N;
00103   int LDA1 = LDA;
00104   double OneNorm1;
00105   bool Upper = false;
00106   
00107   Epetra_SerialSpdDenseSolver solver;
00108   Epetra_SerialSymDenseMatrix * Matrix;
00109   for (int kk=0; kk<2; kk++) {
00110     for (i=1; i<=N; i++) {
00111       GenerateHilbert(A, LDA, i);
00112       OneNorm1 = 0.0;
00113       for (j=1; j<=i; j++) OneNorm1 += 1.0/((double) j); // 1-Norm = 1 + 1/2 + ...+1/n
00114       
00115       if (kk==0) {
00116   Matrix = new Epetra_SerialSymDenseMatrix(View, A, LDA, i);
00117   LDA1 = LDA;
00118       }
00119       else {
00120   Matrix = new Epetra_SerialSymDenseMatrix(Copy, A, LDA, i);
00121   LDA1 = i;
00122       }
00123       GenerateHilbert(A1, LDA1, i);
00124   
00125       if (kk==1) {
00126   solver.FactorWithEquilibration(true);
00127   Matrix->SetUpper();
00128   Upper = true;
00129   solver.SolveToRefinedSolution(false);
00130       } 
00131       
00132       for (k=0; k<NRHS; k++)
00133   for (j=0; j<i; j++) {
00134     B[j+k*LDB] = 1.0/((double) (k+3)*(j+3));
00135     B1[j+k*LDB1] = B[j+k*LDB1];
00136   }
00137       Epetra_SerialDenseMatrix Epetra_B(View, B, LDB, i, NRHS);
00138       Epetra_SerialDenseMatrix Epetra_X(View, X, LDX, i, NRHS);
00139       solver.SetMatrix(*Matrix);
00140       solver.SetVectors(Epetra_X, Epetra_B);
00141        
00142       ierr = check(solver, A1, LDA1,  i, NRHS, OneNorm1, B1, LDB1,  X1, LDX1, Upper, verbose);
00143       assert (ierr>-1);
00144       delete Matrix;
00145       if (ierr!=0) {
00146   if (verbose) cout << "Factorization failed due to bad conditioning.  This is normal if SCOND is small." 
00147         << endl;
00148   break;
00149       } 
00150     }
00151   }
00152 
00153   delete [] A;
00154   delete [] A1;
00155   delete [] X;
00156   delete [] X1;
00157   delete [] B;
00158   delete [] B1;
00159 
00161   // Now test norms and scaling functions
00163 
00164   Epetra_SerialSymDenseMatrix D;
00165   double ScalarA = 2.0;
00166 
00167   int DM = 10;
00168   int DN = 10;
00169   D.Shape(DM);
00170   for (j=0; j<DN; j++)
00171     for (i=0; i<DM; i++) D[j][i] = (double) (1+i+j*DM) ;
00172 
00173   //cout << D << endl;
00174 
00175   double NormInfD_ref = (double)(DM*(DN*(DN+1))/2);
00176   double NormOneD_ref = NormInfD_ref;
00177 
00178   double NormInfD = D.NormInf();
00179   double NormOneD = D.NormOne(); 
00180 
00181   if (verbose) {
00182     cout << " *** Before scaling *** " << endl
00183    << " Computed one-norm of test matrix = " << NormOneD << endl
00184    << " Expected one-norm                = " << NormOneD_ref << endl
00185    << " Computed inf-norm of test matrix = " << NormInfD << endl
00186    << " Expected inf-norm                = " << NormInfD_ref << endl;
00187   }
00188   D.Scale(ScalarA); // Scale entire D matrix by this value
00189 
00190   //cout << D << endl;
00191 
00192   NormInfD = D.NormInf();
00193   NormOneD = D.NormOne();
00194   if (verbose) {
00195     cout << " *** After scaling *** " << endl
00196    << " Computed one-norm of test matrix = " << NormOneD << endl
00197    << " Expected one-norm                = " << NormOneD_ref*ScalarA << endl
00198    << " Computed inf-norm of test matrix = " << NormInfD << endl
00199    << " Expected inf-norm                = " << NormInfD_ref*ScalarA << endl;
00200   }
00201 
00202   
00203 
00205   // Now test for larger system, both correctness and performance.
00207 
00208 
00209   N = 2000;
00210   NRHS = 5;
00211   LDA = N;
00212   LDB = N;
00213   LDX = N;
00214 
00215   if (verbose) cout << "\n\nComputing factor of an " << N << " x " << N << " SPD matrix...Please wait.\n\n" << endl;
00216 
00217   // Define A and X
00218 
00219   A = new double[LDA*N];
00220   X = new double[LDB*NRHS];
00221   
00222   for (j=0; j<N; j++) {
00223     for (k=0; k<NRHS; k++) X[j+k*LDX] = 1.0/((double) (j+5+k));
00224     for (i=0; i<N; i++) { 
00225       if (i==j) A[i+j*LDA] = 100.0 + i;
00226       else A[i+j*LDA] = -1.0/((double) (i+10)*(j+10));
00227     }
00228   }
00229 
00230   // Define Epetra_SerialDenseMatrix object
00231 
00232   Epetra_SerialSymDenseMatrix BigMatrix(Copy, A, LDA, N);
00233   Epetra_SerialSymDenseMatrix OrigBigMatrix(View, A, LDA, N);
00234 
00235   Epetra_SerialSpdDenseSolver BigSolver;
00236   BigSolver.FactorWithEquilibration(true);
00237   BigSolver.SetMatrix(BigMatrix);
00238 
00239   // Time factorization
00240 
00241   Epetra_Flops counter;
00242   BigSolver.SetFlopCounter(counter);
00243   Epetra_Time Timer(Comm);
00244   double tstart = Timer.ElapsedTime();
00245   ierr = BigSolver.Factor();
00246   if (ierr!=0 && verbose) cout << "Error in factorization = "<<ierr<< endl;
00247   assert(ierr==0);
00248   double time = Timer.ElapsedTime() - tstart;
00249 
00250   double FLOPS = counter.Flops();
00251   double MFLOPS = FLOPS/time/1000000.0; 
00252   if (verbose) cout << "MFLOPS for Factorization = " << MFLOPS << endl; 
00253 
00254   // Define Left hand side and right hand side 
00255   Epetra_SerialDenseMatrix LHS(View, X, LDX, N, NRHS);
00256   Epetra_SerialDenseMatrix RHS;
00257   RHS.Shape(N,NRHS); // Allocate RHS
00258 
00259   // Compute RHS from A and X
00260 
00261   Epetra_Flops RHS_counter;
00262   RHS.SetFlopCounter(RHS_counter);
00263   tstart = Timer.ElapsedTime();
00264   RHS.Multiply('L', 1.0, OrigBigMatrix, LHS, 0.0); // Symmetric Matrix-multiply
00265   time = Timer.ElapsedTime() - tstart;
00266 
00267   Epetra_SerialDenseMatrix OrigRHS = RHS;
00268 
00269   FLOPS = RHS_counter.Flops();
00270   MFLOPS = FLOPS/time/1000000.0;
00271   if (verbose) cout << "MFLOPS to build RHS (NRHS = " << NRHS <<") = " << MFLOPS << endl;
00272 
00273   // Set LHS and RHS and solve
00274   BigSolver.SetVectors(LHS, RHS);
00275 
00276   tstart = Timer.ElapsedTime();
00277   ierr = BigSolver.Solve();
00278   if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
00279   else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
00280   assert(ierr>=0);
00281   time = Timer.ElapsedTime() - tstart;
00282 
00283   FLOPS = BigSolver.Flops();
00284   MFLOPS = FLOPS/time/1000000.0;
00285   if (verbose) cout << "MFLOPS for Solve (NRHS = " << NRHS <<") = " << MFLOPS << endl;
00286 
00287   double * resid = new double[NRHS];
00288   bool OK = Residual(N, NRHS, A, LDA, BigSolver.X(), BigSolver.LDX(), 
00289          OrigRHS.A(), OrigRHS.LDA(), resid);
00290 
00291   if (verbose) {
00292     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00293     for (i=0; i<NRHS; i++)
00294       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00295     cout  << endl;
00296   }
00297 
00298   // Solve again using the Epetra_SerialDenseVector class for LHS and RHS
00299 
00300   Epetra_SerialDenseVector X2;
00301   Epetra_SerialDenseVector B2;
00302   X2.Size(BigMatrix.N());
00303   B2.Size(BigMatrix.M());
00304   int length = BigMatrix.N();
00305   {for (int kk=0; kk<length; kk++) X2[kk] = ((double ) kk)/ ((double) length);} // Define entries of X2
00306 
00307   RHS_counter.ResetFlops();
00308   B2.SetFlopCounter(RHS_counter);
00309   tstart = Timer.ElapsedTime();
00310   B2.Multiply('N', 'N', 1.0, OrigBigMatrix, X2, 0.0); // Define B2 = A*X2
00311   time = Timer.ElapsedTime() - tstart;
00312 
00313   Epetra_SerialDenseVector OrigB2 = B2;
00314 
00315   FLOPS = RHS_counter.Flops();
00316   MFLOPS = FLOPS/time/1000000.0;
00317   if (verbose) cout << "MFLOPS to build single RHS = " << MFLOPS << endl;
00318 
00319   // Set LHS and RHS and solve
00320   BigSolver.SetVectors(X2, B2);
00321 
00322   tstart = Timer.ElapsedTime();
00323   ierr = BigSolver.Solve();
00324   time = Timer.ElapsedTime() - tstart;
00325   if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
00326   else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
00327   assert(ierr>=0);
00328 
00329   FLOPS = counter.Flops();
00330   MFLOPS = FLOPS/time/1000000.0;
00331   if (verbose) cout << "MFLOPS to solve single RHS = " << MFLOPS << endl;
00332 
00333   OK = Residual(N, 1, A, LDA, BigSolver.X(), BigSolver.LDX(), OrigB2.A(), 
00334     OrigB2.LDA(), resid);
00335 
00336   if (verbose) {
00337     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00338       cout << "Residual = "<< resid[0] << endl;
00339   }
00340   delete [] resid;
00341   delete [] A;
00342   delete [] X;
00343 
00345   // Now test default constructor and index operators
00347 
00348   N = 5;
00349   Epetra_SerialSymDenseMatrix C; // Implicit call to default constructor, should not need to call destructor
00350   C.Shape(5); // Make it 5 by 5
00351   double * C1 = new double[N*N];
00352   GenerateHilbert(C1, N, N); // Generate Hilber matrix
00353 
00354   C1[1+2*N] = 1000.0;  // Make matrix nonsymmetric
00355 
00356   // Fill values of C with Hilbert values
00357   for (i=0; i<N; i++) 
00358     for (j=0; j<N; j++)
00359       C(i,j) = C1[i+j*N];
00360 
00361   // Test if values are correctly written and read
00362   for (i=0; i<N; i++) 
00363     for (j=0; j<N; j++) {
00364       assert(C(i,j) == C1[i+j*N]);
00365       assert(C(i,j) == C[j][i]);
00366     }
00367 
00368   if (verbose)
00369     cout << "Default constructor and index operator check OK.  Values of Hilbert matrix = " 
00370    << endl << C << endl
00371    << "Values should be 1/(i+j+1), except value (1,2) should be 1000" << endl; 
00372   
00373   delete [] C1;
00374 
00375 
00376 #ifdef EPETRA_MPI
00377   MPI_Finalize() ;
00378 #endif
00379 
00380 /* end main
00381 */
00382 return ierr ;
00383 }
00384 
00385 int check(Epetra_SerialSpdDenseSolver &solver, double * A1, int LDA1, 
00386     int N1, int NRHS1, double OneNorm1, 
00387     double * B1, int LDB1, 
00388     double * X1, int LDX1,
00389     bool Upper, bool verbose) {  
00390   (void)OneNorm1;
00391   int i;
00392   bool OK;
00393   // Test query functions
00394 
00395   int M= solver.M();
00396   if (verbose) cout << "\n\nNumber of Rows = " << M << endl<< endl;
00397   assert(M==N1);
00398 
00399   int N= solver.N();
00400   if (verbose) cout << "\n\nNumber of Equations = " << N << endl<< endl;
00401   assert(N==N1);
00402 
00403   int LDA = solver.LDA();
00404   if (verbose) cout << "\n\nLDA = " << LDA << endl<< endl;
00405   assert(LDA==LDA1);
00406 
00407   int LDB = solver.LDB();
00408   if (verbose) cout << "\n\nLDB = " << LDB << endl<< endl;
00409   assert(LDB==LDB1);
00410 
00411   int LDX = solver.LDX();
00412   if (verbose) cout << "\n\nLDX = " << LDX << endl<< endl;
00413   assert(LDX==LDX1);
00414 
00415   int NRHS = solver.NRHS();
00416   if (verbose) cout << "\n\nNRHS = " << NRHS << endl<< endl;
00417   assert(NRHS==NRHS1);
00418 
00419   assert(solver.ANORM()==-1.0);
00420   assert(solver.RCOND()==-1.0);
00421   if (!solver.A_Equilibrated() && !solver.B_Equilibrated()) {
00422     assert(solver.SCOND()==-1.0);
00423     assert(solver.AMAX()==-1.0);
00424   }
00425 
00426 
00427   // Other binary tests
00428 
00429   assert(!solver.Factored()); 
00430   assert(solver.SymMatrix()->Upper()==Upper);
00431   assert(!solver.SolutionErrorsEstimated());
00432   assert(!solver.Inverted());
00433   assert(!solver.ReciprocalConditionEstimated());
00434   assert(!solver.Solved());
00435 
00436   assert(!solver.SolutionRefined());
00437       
00438   //cout << "Matrix before factorization " << endl << *solver.SymMatrix() << endl << endl;
00439   
00440   int ierr = solver.Factor();
00441   //cout << "Matrix after factorization " << endl << *solver.SymMatrix() << endl << endl;
00442   //cout << "Factor after factorization " << endl << *solver.SymFactoredMatrix() << endl << endl;
00443   assert(ierr>-1);
00444   if (ierr!=0) return(ierr); // Factorization failed due to poor conditioning.
00445   double rcond;
00446   assert(solver.ReciprocalConditionEstimate(rcond)==0);
00447   if (verbose) {
00448     
00449     double rcond1 = 1.0/std::exp(3.5*((double)N));
00450     if (N==1) rcond1 = 1.0;
00451     cout << "\n\nSCOND = "<< rcond << " should be approx = " 
00452         << rcond1 << endl << endl;
00453   }
00454   
00455   ierr = solver.Solve();
00456   assert(ierr>-1);
00457   if (ierr!=0 && verbose) cout << "LAPACK rules suggest system should be equilibrated." << endl;
00458 
00459   assert(solver.Factored());
00460   assert(solver.SymMatrix()->Upper()==Upper);
00461   assert(solver.ReciprocalConditionEstimated());
00462   assert(solver.Solved());
00463 
00464   if (solver.SolutionErrorsEstimated()) {
00465     if (verbose) {
00466       cout << "\n\nFERR[0] = "<< solver.FERR()[0] << endl;
00467       cout << "\n\nBERR[0] = "<< solver.BERR()[0] << endl<< endl;
00468     }
00469   }
00470   
00471   double * resid = new double[NRHS];
00472   OK = Residual(N, NRHS, A1, LDA1, solver.X(), solver.LDX(), B1, LDB1, resid);
00473   if (verbose) {
00474     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00475     cout << "\n\nResiduals using factorization to solve" << endl;
00476     for (i=0; i<NRHS; i++)
00477       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00478     cout  << endl;
00479   }
00480 
00481 
00482   ierr = solver.Invert();
00483   assert(ierr>-1);
00484 
00485   assert(solver.Inverted());
00486   assert(!solver.Factored());
00487   
00488   Epetra_SerialDenseMatrix RHS1(Copy, B1, LDB1, N, NRHS);
00489   Epetra_SerialDenseMatrix LHS1(Copy, X1, LDX1, N, NRHS);
00490   assert(solver.SetVectors(LHS1, RHS1)==0);
00491   assert(!solver.Solved());
00492 
00493   assert(solver.Solve()>-1);
00494    
00495   
00496 
00497   OK = Residual(N, NRHS, A1, LDA1, solver.X(), solver.LDX(), B1, LDB1, resid);
00498 
00499   if (verbose) {
00500     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00501     cout << "Residuals using inverse to solve" << endl;
00502     for (i=0; i<NRHS; i++)
00503       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00504     cout  << endl;
00505   }
00506   delete [] resid;
00507   
00508 
00509   return(0);
00510 }
00511 
00512  void GenerateHilbert(double *A, int LDA, int N) {
00513    for (int j=0; j<N; j++)
00514      for (int i=0; i<N; i++)
00515        A[i+j*LDA] = 1.0/((double)(i+j+1));
00516    return;
00517  }
00518 
00519 bool Residual( int N, int NRHS, double * A, int LDA, 
00520          double * X, int LDX, double * B, int LDB, double * resid) {
00521 
00522   Epetra_BLAS Blas;
00523   char Transa = 'N';
00524   Blas.GEMM(Transa, 'N', N, NRHS, N, -1.0, A, LDA,
00525       X, LDX, 1.0, B, LDB);
00526   bool OK = true;
00527   for (int i=0; i<NRHS; i++) {
00528     resid[i] = Blas.NRM2(N, B+i*LDB);
00529     if (resid[i]>1.0E-7) OK = false;
00530   }
00531 
00532   return(OK);
00533 }

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1