Epetra Package Browser (Single Doxygen Collection) Development
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00043 #include "Epetra_Map.h"
00044 #include "Epetra_Time.h"
00045 #include "Epetra_SerialSymDenseMatrix.h" 
00046 #include "Epetra_SerialDenseMatrix.h" 
00047 #include "Epetra_SerialDenseVector.h"
00048 #include "Epetra_SerialSpdDenseSolver.h"
00049 #ifdef EPETRA_MPI
00050 #include "Epetra_MpiComm.h"
00051 #include <mpi.h>
00052 #endif
00053 #include "Epetra_SerialComm.h"
00054 #include "Epetra_Version.h"
00055 
00056 // prototypes
00057 
00058 int check(Epetra_SerialSpdDenseSolver & solver, double * A1, int LDA,
00059     int N1, int NRHS1, double OneNorm1, 
00060     double * B1, int LDB1, 
00061     double * X1, int LDX1,
00062     bool Upper, bool verbose);
00063 
00064 void GenerateHilbert(double *A, int LDA, int N);
00065 
00066 bool Residual( int N, int NRHS, double * A, int LDA,
00067          double * X, int LDX, double * B, int LDB, double * resid);
00068 
00069  
00070 int main(int argc, char *argv[])
00071 {
00072   int ierr = 0, i, j, k;
00073 
00074 #ifdef EPETRA_MPI
00075   MPI_Init(&argc,&argv);
00076   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00077 #else
00078   Epetra_SerialComm Comm;
00079 #endif
00080 
00081   bool verbose = false;
00082 
00083   // Check if we should print results to standard out
00084   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00085 
00086   if(verbose && Comm.MyPID()==0)
00087     cout << Epetra_Version() << endl << endl;
00088 
00089   int rank = Comm.MyPID();
00090   //  char tmp;
00091   //  if (rank==0) cout << "Press any key to continue..."<< endl;
00092   //  if (rank==0) cin >> tmp;
00093   //  Comm.Barrier();
00094 
00095   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00096   if (verbose) cout << Comm <<endl;
00097 
00098   //  bool verbose1 = verbose;
00099 
00100   // Redefine verbose to only print on PE 0
00101   if (verbose && rank!=0) verbose = false;
00102 
00103   int N = 20;
00104   int NRHS = 4;
00105   double * A = new double[N*N];
00106   double * A1 = new double[N*N];
00107   double * X = new double[(N+1)*NRHS];
00108   double * X1 = new double[(N+1)*NRHS];
00109   int LDX = N+1;
00110   int LDX1 = N+1;
00111   double * B = new double[N*NRHS];
00112   double * B1 = new double[N*NRHS];
00113   int LDB = N;
00114   int LDB1 = N;
00115 
00116   int LDA = N;
00117   int LDA1 = LDA;
00118   double OneNorm1;
00119   bool Upper = false;
00120   
00121   Epetra_SerialSpdDenseSolver solver;
00122   Epetra_SerialSymDenseMatrix * Matrix;
00123   for (int kk=0; kk<2; kk++) {
00124     for (i=1; i<=N; i++) {
00125       GenerateHilbert(A, LDA, i);
00126       OneNorm1 = 0.0;
00127       for (j=1; j<=i; j++) OneNorm1 += 1.0/((double) j); // 1-Norm = 1 + 1/2 + ...+1/n
00128       
00129       if (kk==0) {
00130   Matrix = new Epetra_SerialSymDenseMatrix(View, A, LDA, i);
00131   LDA1 = LDA;
00132       }
00133       else {
00134   Matrix = new Epetra_SerialSymDenseMatrix(Copy, A, LDA, i);
00135   LDA1 = i;
00136       }
00137       GenerateHilbert(A1, LDA1, i);
00138   
00139       if (kk==1) {
00140   solver.FactorWithEquilibration(true);
00141   Matrix->SetUpper();
00142   Upper = true;
00143   solver.SolveToRefinedSolution(false);
00144       } 
00145       
00146       for (k=0; k<NRHS; k++)
00147   for (j=0; j<i; j++) {
00148     B[j+k*LDB] = 1.0/((double) (k+3)*(j+3));
00149     B1[j+k*LDB1] = B[j+k*LDB1];
00150   }
00151       Epetra_SerialDenseMatrix Epetra_B(View, B, LDB, i, NRHS);
00152       Epetra_SerialDenseMatrix Epetra_X(View, X, LDX, i, NRHS);
00153       solver.SetMatrix(*Matrix);
00154       solver.SetVectors(Epetra_X, Epetra_B);
00155        
00156       ierr = check(solver, A1, LDA1,  i, NRHS, OneNorm1, B1, LDB1,  X1, LDX1, Upper, verbose);
00157       assert (ierr>-1);
00158       delete Matrix;
00159       if (ierr!=0) {
00160   if (verbose) cout << "Factorization failed due to bad conditioning.  This is normal if SCOND is small." 
00161         << endl;
00162   break;
00163       } 
00164     }
00165   }
00166 
00167   delete [] A;
00168   delete [] A1;
00169   delete [] X;
00170   delete [] X1;
00171   delete [] B;
00172   delete [] B1;
00173 
00175   // Now test norms and scaling functions
00177 
00178   Epetra_SerialSymDenseMatrix D;
00179   double ScalarA = 2.0;
00180 
00181   int DM = 10;
00182   int DN = 10;
00183   D.Shape(DM);
00184   for (j=0; j<DN; j++)
00185     for (i=0; i<DM; i++) D[j][i] = (double) (1+i+j*DM) ;
00186 
00187   //cout << D << endl;
00188 
00189   double NormInfD_ref = (double)(DM*(DN*(DN+1))/2);
00190   double NormOneD_ref = NormInfD_ref;
00191 
00192   double NormInfD = D.NormInf();
00193   double NormOneD = D.NormOne(); 
00194 
00195   if (verbose) {
00196     cout << " *** Before scaling *** " << endl
00197    << " Computed one-norm of test matrix = " << NormOneD << endl
00198    << " Expected one-norm                = " << NormOneD_ref << endl
00199    << " Computed inf-norm of test matrix = " << NormInfD << endl
00200    << " Expected inf-norm                = " << NormInfD_ref << endl;
00201   }
00202   D.Scale(ScalarA); // Scale entire D matrix by this value
00203 
00204   //cout << D << endl;
00205 
00206   NormInfD = D.NormInf();
00207   NormOneD = D.NormOne();
00208   if (verbose) {
00209     cout << " *** After scaling *** " << endl
00210    << " Computed one-norm of test matrix = " << NormOneD << endl
00211    << " Expected one-norm                = " << NormOneD_ref*ScalarA << endl
00212    << " Computed inf-norm of test matrix = " << NormInfD << endl
00213    << " Expected inf-norm                = " << NormInfD_ref*ScalarA << endl;
00214   }
00215 
00216   
00217 
00219   // Now test for larger system, both correctness and performance.
00221 
00222 
00223   N = 2000;
00224   NRHS = 5;
00225   LDA = N;
00226   LDB = N;
00227   LDX = N;
00228 
00229   if (verbose) cout << "\n\nComputing factor of an " << N << " x " << N << " SPD matrix...Please wait.\n\n" << endl;
00230 
00231   // Define A and X
00232 
00233   A = new double[LDA*N];
00234   X = new double[LDB*NRHS];
00235   
00236   for (j=0; j<N; j++) {
00237     for (k=0; k<NRHS; k++) X[j+k*LDX] = 1.0/((double) (j+5+k));
00238     for (i=0; i<N; i++) { 
00239       if (i==j) A[i+j*LDA] = 100.0 + i;
00240       else A[i+j*LDA] = -1.0/((double) (i+10)*(j+10));
00241     }
00242   }
00243 
00244   // Define Epetra_SerialDenseMatrix object
00245 
00246   Epetra_SerialSymDenseMatrix BigMatrix(Copy, A, LDA, N);
00247   Epetra_SerialSymDenseMatrix OrigBigMatrix(View, A, LDA, N);
00248 
00249   Epetra_SerialSpdDenseSolver BigSolver;
00250   BigSolver.FactorWithEquilibration(true);
00251   BigSolver.SetMatrix(BigMatrix);
00252 
00253   // Time factorization
00254 
00255   Epetra_Flops counter;
00256   BigSolver.SetFlopCounter(counter);
00257   Epetra_Time Timer(Comm);
00258   double tstart = Timer.ElapsedTime();
00259   ierr = BigSolver.Factor();
00260   if (ierr!=0 && verbose) cout << "Error in factorization = "<<ierr<< endl;
00261   assert(ierr==0);
00262   double time = Timer.ElapsedTime() - tstart;
00263 
00264   double FLOPS = counter.Flops();
00265   double MFLOPS = FLOPS/time/1000000.0; 
00266   if (verbose) cout << "MFLOPS for Factorization = " << MFLOPS << endl; 
00267 
00268   // Define Left hand side and right hand side 
00269   Epetra_SerialDenseMatrix LHS(View, X, LDX, N, NRHS);
00270   Epetra_SerialDenseMatrix RHS;
00271   RHS.Shape(N,NRHS); // Allocate RHS
00272 
00273   // Compute RHS from A and X
00274 
00275   Epetra_Flops RHS_counter;
00276   RHS.SetFlopCounter(RHS_counter);
00277   tstart = Timer.ElapsedTime();
00278   RHS.Multiply('L', 1.0, OrigBigMatrix, LHS, 0.0); // Symmetric Matrix-multiply
00279   time = Timer.ElapsedTime() - tstart;
00280 
00281   Epetra_SerialDenseMatrix OrigRHS = RHS;
00282 
00283   FLOPS = RHS_counter.Flops();
00284   MFLOPS = FLOPS/time/1000000.0;
00285   if (verbose) cout << "MFLOPS to build RHS (NRHS = " << NRHS <<") = " << MFLOPS << endl;
00286 
00287   // Set LHS and RHS and solve
00288   BigSolver.SetVectors(LHS, RHS);
00289 
00290   tstart = Timer.ElapsedTime();
00291   ierr = BigSolver.Solve();
00292   if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
00293   else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
00294   assert(ierr>=0);
00295   time = Timer.ElapsedTime() - tstart;
00296 
00297   FLOPS = BigSolver.Flops();
00298   MFLOPS = FLOPS/time/1000000.0;
00299   if (verbose) cout << "MFLOPS for Solve (NRHS = " << NRHS <<") = " << MFLOPS << endl;
00300 
00301   double * resid = new double[NRHS];
00302   bool OK = Residual(N, NRHS, A, LDA, BigSolver.X(), BigSolver.LDX(), 
00303          OrigRHS.A(), OrigRHS.LDA(), resid);
00304 
00305   if (verbose) {
00306     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00307     for (i=0; i<NRHS; i++)
00308       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00309     cout  << endl;
00310   }
00311 
00312   // Solve again using the Epetra_SerialDenseVector class for LHS and RHS
00313 
00314   Epetra_SerialDenseVector X2;
00315   Epetra_SerialDenseVector B2;
00316   X2.Size(BigMatrix.N());
00317   B2.Size(BigMatrix.M());
00318   int length = BigMatrix.N();
00319   {for (int kk=0; kk<length; kk++) X2[kk] = ((double ) kk)/ ((double) length);} // Define entries of X2
00320 
00321   RHS_counter.ResetFlops();
00322   B2.SetFlopCounter(RHS_counter);
00323   tstart = Timer.ElapsedTime();
00324   B2.Multiply('N', 'N', 1.0, OrigBigMatrix, X2, 0.0); // Define B2 = A*X2
00325   time = Timer.ElapsedTime() - tstart;
00326 
00327   Epetra_SerialDenseVector OrigB2 = B2;
00328 
00329   FLOPS = RHS_counter.Flops();
00330   MFLOPS = FLOPS/time/1000000.0;
00331   if (verbose) cout << "MFLOPS to build single RHS = " << MFLOPS << endl;
00332 
00333   // Set LHS and RHS and solve
00334   BigSolver.SetVectors(X2, B2);
00335 
00336   tstart = Timer.ElapsedTime();
00337   ierr = BigSolver.Solve();
00338   time = Timer.ElapsedTime() - tstart;
00339   if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
00340   else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
00341   assert(ierr>=0);
00342 
00343   FLOPS = counter.Flops();
00344   MFLOPS = FLOPS/time/1000000.0;
00345   if (verbose) cout << "MFLOPS to solve single RHS = " << MFLOPS << endl;
00346 
00347   OK = Residual(N, 1, A, LDA, BigSolver.X(), BigSolver.LDX(), OrigB2.A(), 
00348     OrigB2.LDA(), resid);
00349 
00350   if (verbose) {
00351     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00352       cout << "Residual = "<< resid[0] << endl;
00353   }
00354   delete [] resid;
00355   delete [] A;
00356   delete [] X;
00357 
00359   // Now test default constructor and index operators
00361 
00362   N = 5;
00363   Epetra_SerialSymDenseMatrix C; // Implicit call to default constructor, should not need to call destructor
00364   C.Shape(5); // Make it 5 by 5
00365   double * C1 = new double[N*N];
00366   GenerateHilbert(C1, N, N); // Generate Hilber matrix
00367 
00368   C1[1+2*N] = 1000.0;  // Make matrix nonsymmetric
00369 
00370   // Fill values of C with Hilbert values
00371   for (i=0; i<N; i++) 
00372     for (j=0; j<N; j++)
00373       C(i,j) = C1[i+j*N];
00374 
00375   // Test if values are correctly written and read
00376   for (i=0; i<N; i++) 
00377     for (j=0; j<N; j++) {
00378       assert(C(i,j) == C1[i+j*N]);
00379       assert(C(i,j) == C[j][i]);
00380     }
00381 
00382   if (verbose)
00383     cout << "Default constructor and index operator check OK.  Values of Hilbert matrix = " 
00384    << endl << C << endl
00385    << "Values should be 1/(i+j+1), except value (1,2) should be 1000" << endl; 
00386   
00387   delete [] C1;
00388 
00389 
00390 #ifdef EPETRA_MPI
00391   MPI_Finalize() ;
00392 #endif
00393 
00394 /* end main
00395 */
00396 return ierr ;
00397 }
00398 
00399 int check(Epetra_SerialSpdDenseSolver &solver, double * A1, int LDA1, 
00400     int N1, int NRHS1, double OneNorm1, 
00401     double * B1, int LDB1, 
00402     double * X1, int LDX1,
00403     bool Upper, bool verbose) {  
00404   (void)OneNorm1;
00405   int i;
00406   bool OK;
00407   // Test query functions
00408 
00409   int M= solver.M();
00410   if (verbose) cout << "\n\nNumber of Rows = " << M << endl<< endl;
00411   assert(M==N1);
00412 
00413   int N= solver.N();
00414   if (verbose) cout << "\n\nNumber of Equations = " << N << endl<< endl;
00415   assert(N==N1);
00416 
00417   int LDA = solver.LDA();
00418   if (verbose) cout << "\n\nLDA = " << LDA << endl<< endl;
00419   assert(LDA==LDA1);
00420 
00421   int LDB = solver.LDB();
00422   if (verbose) cout << "\n\nLDB = " << LDB << endl<< endl;
00423   assert(LDB==LDB1);
00424 
00425   int LDX = solver.LDX();
00426   if (verbose) cout << "\n\nLDX = " << LDX << endl<< endl;
00427   assert(LDX==LDX1);
00428 
00429   int NRHS = solver.NRHS();
00430   if (verbose) cout << "\n\nNRHS = " << NRHS << endl<< endl;
00431   assert(NRHS==NRHS1);
00432 
00433   assert(solver.ANORM()==-1.0);
00434   assert(solver.RCOND()==-1.0);
00435   if (!solver.A_Equilibrated() && !solver.B_Equilibrated()) {
00436     assert(solver.SCOND()==-1.0);
00437     assert(solver.AMAX()==-1.0);
00438   }
00439 
00440 
00441   // Other binary tests
00442 
00443   assert(!solver.Factored()); 
00444   assert(solver.SymMatrix()->Upper()==Upper);
00445   assert(!solver.SolutionErrorsEstimated());
00446   assert(!solver.Inverted());
00447   assert(!solver.ReciprocalConditionEstimated());
00448   assert(!solver.Solved());
00449 
00450   assert(!solver.SolutionRefined());
00451       
00452   //cout << "Matrix before factorization " << endl << *solver.SymMatrix() << endl << endl;
00453   
00454   int ierr = solver.Factor();
00455   //cout << "Matrix after factorization " << endl << *solver.SymMatrix() << endl << endl;
00456   //cout << "Factor after factorization " << endl << *solver.SymFactoredMatrix() << endl << endl;
00457   assert(ierr>-1);
00458   if (ierr!=0) return(ierr); // Factorization failed due to poor conditioning.
00459   double rcond;
00460   ierr = solver.ReciprocalConditionEstimate(rcond);
00461   assert(ierr==0);
00462   if (verbose) {
00463     
00464     double rcond1 = 1.0/std::exp(3.5*((double)N));
00465     if (N==1) rcond1 = 1.0;
00466     cout << "\n\nSCOND = "<< rcond << " should be approx = " 
00467         << rcond1 << endl << endl;
00468   }
00469   
00470   ierr = solver.Solve();
00471   assert(ierr>-1);
00472   if (ierr!=0 && verbose) cout << "LAPACK rules suggest system should be equilibrated." << endl;
00473 
00474   assert(solver.Factored());
00475   assert(solver.SymMatrix()->Upper()==Upper);
00476   assert(solver.ReciprocalConditionEstimated());
00477   assert(solver.Solved());
00478 
00479   if (solver.SolutionErrorsEstimated()) {
00480     if (verbose) {
00481       cout << "\n\nFERR[0] = "<< solver.FERR()[0] << endl;
00482       cout << "\n\nBERR[0] = "<< solver.BERR()[0] << endl<< endl;
00483     }
00484   }
00485   
00486   double * resid = new double[NRHS];
00487   OK = Residual(N, NRHS, A1, LDA1, solver.X(), solver.LDX(), B1, LDB1, resid);
00488   if (verbose) {
00489     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00490     cout << "\n\nResiduals using factorization to solve" << endl;
00491     for (i=0; i<NRHS; i++)
00492       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00493     cout  << endl;
00494   }
00495 
00496 
00497   ierr = solver.Invert();
00498   assert(ierr>-1);
00499 
00500   assert(solver.Inverted());
00501   assert(!solver.Factored());
00502   
00503   Epetra_SerialDenseMatrix RHS1(Copy, B1, LDB1, N, NRHS);
00504   Epetra_SerialDenseMatrix LHS1(Copy, X1, LDX1, N, NRHS);
00505   assert(solver.SetVectors(LHS1, RHS1)==0);
00506   assert(!solver.Solved());
00507 
00508   assert(solver.Solve()>-1);
00509    
00510   
00511 
00512   OK = Residual(N, NRHS, A1, LDA1, solver.X(), solver.LDX(), B1, LDB1, resid);
00513 
00514   if (verbose) {
00515     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00516     cout << "Residuals using inverse to solve" << endl;
00517     for (i=0; i<NRHS; i++)
00518       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00519     cout  << endl;
00520   }
00521   delete [] resid;
00522   
00523 
00524   return(0);
00525 }
00526 
00527  void GenerateHilbert(double *A, int LDA, int N) {
00528    for (int j=0; j<N; j++)
00529      for (int i=0; i<N; i++)
00530        A[i+j*LDA] = 1.0/((double)(i+j+1));
00531    return;
00532  }
00533 
00534 bool Residual( int N, int NRHS, double * A, int LDA, 
00535          double * X, int LDX, double * B, int LDB, double * resid) {
00536 
00537   Epetra_BLAS Blas;
00538   char Transa = 'N';
00539   Blas.GEMM(Transa, 'N', N, NRHS, N, -1.0, A, LDA,
00540       X, LDX, 1.0, B, LDB);
00541   bool OK = true;
00542   for (int i=0; i<NRHS; i++) {
00543     resid[i] = Blas.NRM2(N, B+i*LDB);
00544     if (resid[i]>1.0E-7) OK = false;
00545   }
00546 
00547   return(OK);
00548 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines