test/SerialDense/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_SerialDenseMatrix.h" 
00032 #include "Epetra_SerialDenseVector.h"
00033 #include "Epetra_SerialDenseSolver.h"
00034 #ifdef EPETRA_MPI
00035 #include "Epetra_MpiComm.h"
00036 #include <mpi.h>
00037 #endif
00038 #include "Epetra_SerialComm.h"
00039 #include "../epetra_test_err.h"
00040 #include "Epetra_Version.h"
00041 
00042 // prototypes
00043 
00044 int check(Epetra_SerialDenseSolver & solver, double * A1, int LDA,
00045     int N1, int NRHS1, double OneNorm1, 
00046     double * B1, int LDB1, 
00047     double * X1, int LDX1,
00048     bool Transpose, bool verbose);
00049 
00050 void GenerateHilbert(double *A, int LDA, int N);
00051 
00052 bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
00053          double * X, int LDX, double * B, int LDB, double * resid);
00054 
00055 int matrixCpyCtr(bool verbose, bool debug);
00056 int matrixAssignment(bool verbose, bool debug);
00057 void printHeading(const char* heading);
00058 double* getRandArray(int length);
00059 void printMat(const char* name, Epetra_SerialDenseMatrix& matrix);
00060 void printArray(double* array, int length);
00061 bool identicalSignatures(Epetra_SerialDenseMatrix& a, Epetra_SerialDenseMatrix& b, bool testLDA = true);
00062 bool seperateData(Epetra_SerialDenseMatrix& a, Epetra_SerialDenseMatrix& b);
00063 
00064  
00065 int main(int argc, char *argv[])
00066 {
00067   int ierr = 0, i, j, k;
00068   bool debug = false;
00069 
00070 #ifdef EPETRA_MPI
00071   MPI_Init(&argc,&argv);
00072   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00073 #else
00074   Epetra_SerialComm Comm;
00075 #endif
00076 
00077   bool verbose = false;
00078 
00079   // Check if we should print results to standard out
00080   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00081 
00082   if (verbose && Comm.MyPID()==0)
00083     cout << Epetra_Version() << endl << endl;
00084 
00085   int rank = Comm.MyPID();
00086   //  char tmp;
00087   //  if (rank==0) cout << "Press any key to continue..."<< endl;
00088   //  if (rank==0) cin >> tmp;
00089   //  Comm.Barrier();
00090 
00091   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00092   if (verbose) cout << Comm <<endl;
00093 
00094   //  bool verbose1 = verbose;
00095 
00096   // Redefine verbose to only print on PE 0
00097   if (verbose && rank!=0) verbose = false;
00098   
00099   int N = 20;
00100   int NRHS = 4;
00101   double * A = new double[N*N];
00102   double * A1 = new double[N*N];
00103   double * X = new double[(N+1)*NRHS];
00104   double * X1 = new double[(N+1)*NRHS];
00105   int LDX = N+1;
00106   int LDX1 = N+1;
00107   double * B = new double[N*NRHS];
00108   double * B1 = new double[N*NRHS];
00109   int LDB = N;
00110   int LDB1 = N;
00111 
00112   int LDA = N;
00113   int LDA1 = LDA;
00114   double OneNorm1;
00115   bool Transpose = false;
00116   
00117   Epetra_SerialDenseSolver solver;
00118   Epetra_SerialDenseMatrix * Matrix;
00119   for (int kk=0; kk<2; kk++) {
00120     for (i=1; i<=N; i++) {
00121       GenerateHilbert(A, LDA, i);
00122       OneNorm1 = 0.0;
00123       for (j=1; j<=i; j++) OneNorm1 += 1.0/((double) j); // 1-Norm = 1 + 1/2 + ...+1/n
00124       
00125       if (kk==0) {
00126   Matrix = new Epetra_SerialDenseMatrix(View, A, LDA, i, i);
00127   LDA1 = LDA;
00128       }
00129       else {
00130   Matrix = new Epetra_SerialDenseMatrix(Copy, A, LDA, i, i);
00131   LDA1 = i;
00132       }
00133       GenerateHilbert(A1, LDA1, i);
00134   
00135       if (kk==1) {
00136   solver.FactorWithEquilibration(true);
00137   solver.SolveWithTranspose(true);
00138   Transpose = true;
00139   solver.SolveToRefinedSolution(true);
00140       }
00141       
00142       for (k=0; k<NRHS; k++)
00143   for (j=0; j<i; j++) {
00144     B[j+k*LDB] = 1.0/((double) (k+3)*(j+3));
00145     B1[j+k*LDB1] = B[j+k*LDB1];
00146   }
00147       Epetra_SerialDenseMatrix Epetra_B(View, B, LDB, i, NRHS);
00148       Epetra_SerialDenseMatrix Epetra_X(View, X, LDX, i, NRHS);
00149       solver.SetMatrix(*Matrix);
00150       solver.SetVectors(Epetra_X, Epetra_B);
00151       
00152       ierr = check(solver, A1, LDA1,  i, NRHS, OneNorm1, B1, LDB1,  X1, LDX1, Transpose, verbose);
00153       assert (ierr>-1);
00154       delete Matrix;
00155       if (ierr!=0) {
00156   if (verbose) cout << "Factorization failed due to bad conditioning.  This is normal if RCOND is small." 
00157         << endl;
00158   break;
00159       } 
00160     }
00161   }
00162 
00163   delete [] A;
00164   delete [] A1;
00165   delete [] X;
00166   delete [] X1;
00167   delete [] B;
00168   delete [] B1;
00169 
00171   // Now test norms and scaling functions
00173 
00174   Epetra_SerialDenseMatrix D;
00175   double ScalarA = 2.0;
00176 
00177   int DM = 10;
00178   int DN = 8;
00179   D.Shape(DM, DN);
00180   for (j=0; j<DN; j++)
00181     for (i=0; i<DM; i++) D[j][i] = (double) (1+i+j*DM) ;
00182 
00183   //cout << D << endl;
00184 
00185   double NormInfD_ref = (double)(DM*(DN*(DN+1))/2);
00186   double NormOneD_ref = (double)((DM*DN*(DM*DN+1))/2 - (DM*(DN-1)*(DM*(DN-1)+1))/2 );
00187 
00188   double NormInfD = D.NormInf();
00189   double NormOneD = D.NormOne();
00190 
00191   if (verbose) {
00192     cout << " *** Before scaling *** " << endl
00193    << " Computed one-norm of test matrix = " << NormOneD << endl
00194    << " Expected one-norm                = " << NormOneD_ref << endl
00195    << " Computed inf-norm of test matrix = " << NormInfD << endl
00196    << " Expected inf-norm                = " << NormInfD_ref << endl;
00197   }
00198   D.Scale(ScalarA); // Scale entire D matrix by this value
00199   NormInfD = D.NormInf();
00200   NormOneD = D.NormOne();
00201   if (verbose) {
00202     cout << " *** After scaling *** " << endl
00203    << " Computed one-norm of test matrix = " << NormOneD << endl
00204    << " Expected one-norm                = " << NormOneD_ref*ScalarA << endl
00205    << " Computed inf-norm of test matrix = " << NormInfD << endl
00206    << " Expected inf-norm                = " << NormInfD_ref*ScalarA << endl;
00207   }
00208 
00209   
00211   // Now test that A.Multiply(false, x, y) produces the same result
00212   // as y.Multiply('N','N', 1.0, A, x, 0.0).
00214 
00215   N = 10;
00216   int M = 10;
00217   LDA = N;
00218   Epetra_SerialDenseMatrix smallA(N, M, false);
00219   Epetra_SerialDenseMatrix x(N, 1, false);
00220   Epetra_SerialDenseMatrix y1(N, 1, false);
00221   Epetra_SerialDenseMatrix y2(N, 1, false);
00222 
00223   for(i=0; i<N; ++i) {
00224     for(j=0; j<M; ++j) {
00225       smallA(i,j) = 1.0*i+2.0*j+1.0;
00226     }
00227     x(i,0) = 1.0;
00228     y1(i,0) = 0.0;
00229     y2(i,0) = 0.0;
00230   }
00231 
00232   int err1 = smallA.Multiply(false, x, y1);
00233   int err2 = y2.Multiply('N','N', 1.0, smallA, x, 0.0);
00234   if (err1 != 0 || err2 != 0) {
00235     if (verbose) cout << "err in Epetra_SerialDenseMatrix::Multiply"<<endl;
00236     return(err1+err2);
00237   }
00238 
00239   for(i=0; i<N; ++i) {
00240     if (y1(i,0) != y2(i,0)) {
00241       if (verbose) cout << "different versions of Multiply don't match."<<endl;
00242       return(-99);
00243     }
00244   }
00245  
00247   // Now test for larger system, both correctness and performance.
00249 
00250 
00251   N = 2000;
00252   NRHS = 5;
00253   LDA = N;
00254   LDB = N;
00255   LDX = N;
00256 
00257   if (verbose) cout << "\n\nComputing factor of an " << N << " x " << N << " general matrix...Please wait.\n\n" << endl;
00258 
00259   // Define A and X
00260 
00261   A = new double[LDA*N];
00262   X = new double[LDB*NRHS];
00263   
00264   for (j=0; j<N; j++) {
00265     for (k=0; k<NRHS; k++) X[j+k*LDX] = 1.0/((double) (j+5+k));
00266     for (i=0; i<N; i++) {
00267       if (i==((j+2)%N)) A[i+j*LDA] = 100.0 + i;
00268       else A[i+j*LDA] = -11.0/((double) (i+5)*(j+2));
00269     }
00270   }
00271 
00272   // Define Epetra_SerialDenseMatrix object
00273 
00274   Epetra_SerialDenseMatrix BigMatrix(Copy, A, LDA, N, N);
00275   Epetra_SerialDenseMatrix OrigBigMatrix(View, A, LDA, N, N);
00276 
00277   Epetra_SerialDenseSolver BigSolver;
00278   BigSolver.FactorWithEquilibration(true);
00279   BigSolver.SetMatrix(BigMatrix);
00280 
00281   // Time factorization
00282 
00283   Epetra_Flops counter;
00284   BigSolver.SetFlopCounter(counter);
00285   Epetra_Time Timer(Comm);
00286   double tstart = Timer.ElapsedTime();
00287   ierr = BigSolver.Factor();
00288   if (ierr!=0 && verbose) cout << "Error in factorization = "<<ierr<< endl;
00289   assert(ierr==0);
00290   double time = Timer.ElapsedTime() - tstart;
00291 
00292   double FLOPS = counter.Flops();
00293   double MFLOPS = FLOPS/time/1000000.0;
00294   if (verbose) cout << "MFLOPS for Factorization = " << MFLOPS << endl;
00295 
00296   // Define Left hand side and right hand side 
00297   Epetra_SerialDenseMatrix LHS(View, X, LDX, N, NRHS);
00298   Epetra_SerialDenseMatrix RHS;
00299   RHS.Shape(N,NRHS); // Allocate RHS
00300 
00301   // Compute RHS from A and X
00302 
00303   Epetra_Flops RHS_counter;
00304   RHS.SetFlopCounter(RHS_counter);
00305   tstart = Timer.ElapsedTime();
00306   RHS.Multiply('N', 'N', 1.0, OrigBigMatrix, LHS, 0.0);
00307   time = Timer.ElapsedTime() - tstart;
00308 
00309   Epetra_SerialDenseMatrix OrigRHS = RHS;
00310 
00311   FLOPS = RHS_counter.Flops();
00312   MFLOPS = FLOPS/time/1000000.0;
00313   if (verbose) cout << "MFLOPS to build RHS (NRHS = " << NRHS <<") = " << MFLOPS << endl;
00314 
00315   // Set LHS and RHS and solve
00316   BigSolver.SetVectors(LHS, RHS);
00317 
00318   tstart = Timer.ElapsedTime();
00319   ierr = BigSolver.Solve();
00320   if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
00321   else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
00322   assert(ierr>=0);
00323   time = Timer.ElapsedTime() - tstart;
00324 
00325   FLOPS = BigSolver.Flops();
00326   MFLOPS = FLOPS/time/1000000.0;
00327   if (verbose) cout << "MFLOPS for Solve (NRHS = " << NRHS <<") = " << MFLOPS << endl;
00328 
00329   double * resid = new double[NRHS];
00330   bool OK = Residual(N, NRHS, A, LDA, BigSolver.Transpose(), BigSolver.X(), BigSolver.LDX(), 
00331          OrigRHS.A(), OrigRHS.LDA(), resid);
00332 
00333   if (verbose) {
00334     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00335     for (i=0; i<NRHS; i++)
00336       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00337     cout  << endl;
00338   }
00339 
00340   // Solve again using the Epetra_SerialDenseVector class for LHS and RHS
00341 
00342   Epetra_SerialDenseVector X2;
00343   Epetra_SerialDenseVector B2;
00344   X2.Size(BigMatrix.N());
00345   B2.Size(BigMatrix.M());
00346   int length = BigMatrix.N();
00347   {for (int kk=0; kk<length; kk++) X2[kk] = ((double ) kk)/ ((double) length);} // Define entries of X2
00348 
00349   RHS_counter.ResetFlops();
00350   B2.SetFlopCounter(RHS_counter);
00351   tstart = Timer.ElapsedTime();
00352   B2.Multiply('N', 'N', 1.0, OrigBigMatrix, X2, 0.0); // Define B2 = A*X2
00353   time = Timer.ElapsedTime() - tstart;
00354 
00355   Epetra_SerialDenseVector OrigB2 = B2;
00356 
00357   FLOPS = RHS_counter.Flops();
00358   MFLOPS = FLOPS/time/1000000.0;
00359   if (verbose) cout << "MFLOPS to build single RHS = " << MFLOPS << endl;
00360 
00361   // Set LHS and RHS and solve
00362   BigSolver.SetVectors(X2, B2);
00363 
00364   tstart = Timer.ElapsedTime();
00365   ierr = BigSolver.Solve();
00366   time = Timer.ElapsedTime() - tstart;
00367   if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
00368   else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
00369   assert(ierr>=0);
00370 
00371   FLOPS = counter.Flops();
00372   MFLOPS = FLOPS/time/1000000.0;
00373   if (verbose) cout << "MFLOPS to solve single RHS = " << MFLOPS << endl;
00374 
00375   OK = Residual(N, 1, A, LDA, BigSolver.Transpose(), BigSolver.X(), BigSolver.LDX(), OrigB2.A(), 
00376     OrigB2.LDA(), resid);
00377 
00378   if (verbose) {
00379     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00380       cout << "Residual = "<< resid[0] << endl;
00381   }
00382   delete [] resid;
00383   delete [] A;
00384   delete [] X;
00385 
00387   // Now test default constructor and index operators
00389 
00390   N = 5;
00391   Epetra_SerialDenseMatrix C; // Implicit call to default constructor, should not need to call destructor
00392   C.Shape(5,5); // Make it 5 by 5
00393   double * C1 = new double[N*N];
00394   GenerateHilbert(C1, N, N); // Generate Hilber matrix
00395 
00396   C1[1+2*N] = 1000.0;  // Make matrix nonsymmetric
00397 
00398   // Fill values of C with Hilbert values
00399   for (i=0; i<N; i++) 
00400     for (j=0; j<N; j++)
00401       C(i,j) = C1[i+j*N];
00402 
00403   // Test if values are correctly written and read
00404   for (i=0; i<N; i++) 
00405     for (j=0; j<N; j++) {
00406       assert(C(i,j) == C1[i+j*N]);
00407       assert(C(i,j) == C[j][i]);
00408     }
00409 
00410   if (verbose)
00411     cout << "Default constructor and index operator check OK.  Values of Hilbert matrix = " 
00412    << endl << C << endl
00413    << "Values should be 1/(i+j+1), except value (1,2) should be 1000" << endl; 
00414   
00415   delete [] C1;
00416 
00417   // now test sized/shaped constructor
00418   Epetra_SerialDenseMatrix shapedMatrix(10, 12);
00419   assert(shapedMatrix.M() == 10);
00420   assert(shapedMatrix.N() == 12);
00421   for(i = 0; i < 10; i++)
00422     for(j = 0; j < 12; j++)
00423       assert(shapedMatrix(i, j) == 0.0);
00424   Epetra_SerialDenseVector sizedVector(20);
00425   assert(sizedVector.Length() == 20);
00426   for(i = 0; i < 20; i++)
00427     assert(sizedVector(i) == 0.0);
00428   if (verbose)
00429     cout << "Shaped/sized constructors check OK." << endl;
00430   
00431   // test Copy/View mode in op= and cpy ctr
00432   int temperr = 0;
00433   temperr = matrixAssignment(verbose, debug);
00434   if(verbose && temperr == 0)
00435     cout << "Operator = checked OK." << endl;
00436   EPETRA_TEST_ERR(temperr, ierr);
00437   temperr = matrixCpyCtr(verbose, debug);
00438   if(verbose && temperr == 0)
00439     cout << "Copy ctr checked OK." << endl;
00440   EPETRA_TEST_ERR(temperr, ierr);
00441   
00442   // Test some vector methods
00443 
00444   Epetra_SerialDenseVector v1(3);
00445   v1[0] = 1.0;
00446   v1[1] = 3.0;
00447   v1[2] = 2.0;
00448 
00449   Epetra_SerialDenseVector v2(3);
00450   v2[0] = 2.0;
00451   v2[1] = 1.0;
00452   v2[2] = -2.0;
00453 
00454   temperr = 0;
00455   if (v1.Norm1()!=6.0) temperr++;
00456   if (fabs(sqrt(14.0)-v1.Norm2())>1.0e-6) temperr++;
00457   if (v1.NormInf()!=3.0) temperr++;
00458   if(verbose && temperr == 0)
00459     cout << "Vector Norms checked OK." << endl;
00460   temperr = 0;
00461   if (v1.Dot(v2)!=1.0) temperr++;
00462   if(verbose && temperr == 0)
00463     cout << "Vector Dot product checked OK." << endl;
00464 
00465 #ifdef EPETRA_MPI
00466   MPI_Finalize() ;
00467 #endif
00468 
00469 /* end main
00470 */
00471 return ierr ;
00472 }
00473 
00474 int check(Epetra_SerialDenseSolver &solver, double * A1, int LDA1, 
00475     int N1, int NRHS1, double OneNorm1, 
00476     double * B1, int LDB1, 
00477     double * X1, int LDX1,
00478     bool Transpose, bool verbose) {  
00479 
00480   int i;
00481   bool OK;
00482   // Test query functions
00483 
00484   int M= solver.M();
00485   if (verbose) cout << "\n\nNumber of Rows = " << M << endl<< endl;
00486   assert(M==N1);
00487 
00488   int N= solver.N();
00489   if (verbose) cout << "\n\nNumber of Equations = " << N << endl<< endl;
00490   assert(N==N1);
00491 
00492   int LDA = solver.LDA();
00493   if (verbose) cout << "\n\nLDA = " << LDA << endl<< endl;
00494   assert(LDA==LDA1);
00495 
00496   int LDB = solver.LDB();
00497   if (verbose) cout << "\n\nLDB = " << LDB << endl<< endl;
00498   assert(LDB==LDB1);
00499 
00500   int LDX = solver.LDX();
00501   if (verbose) cout << "\n\nLDX = " << LDX << endl<< endl;
00502   assert(LDX==LDX1);
00503 
00504   int NRHS = solver.NRHS();
00505   if (verbose) cout << "\n\nNRHS = " << NRHS << endl<< endl;
00506   assert(NRHS==NRHS1);
00507 
00508   assert(solver.ANORM()==-1.0);
00509   assert(solver.RCOND()==-1.0);
00510   if (!solver.A_Equilibrated() && !solver.B_Equilibrated()) {
00511     assert(solver.ROWCND()==-1.0);
00512     assert(solver.COLCND()==-1.0);
00513     assert(solver.AMAX()==-1.0);
00514   }
00515 
00516 
00517   // Other binary tests
00518 
00519   assert(!solver.Factored()); 
00520   assert(solver.Transpose()==Transpose);
00521   assert(!solver.SolutionErrorsEstimated());
00522   assert(!solver.Inverted());
00523   assert(!solver.ReciprocalConditionEstimated());
00524   assert(!solver.Solved());
00525 
00526   assert(!solver.SolutionRefined());
00527       
00528   
00529   int ierr = solver.Factor();
00530   assert(ierr>-1);
00531   if (ierr!=0) return(ierr); // Factorization failed due to poor conditioning.
00532   double rcond;
00533   assert(solver.ReciprocalConditionEstimate(rcond)==0);
00534   if (verbose) {
00535     
00536     double rcond1 = 1.0/exp(3.5*((double)N));
00537     if (N==1) rcond1 = 1.0;
00538     cout << "\n\nRCOND = "<< rcond << " should be approx = " 
00539         << rcond1 << endl << endl;
00540   }
00541   
00542   ierr = solver.Solve();
00543   assert(ierr>-1);
00544   if (ierr!=0 && verbose) cout << "LAPACK rules suggest system should be equilibrated." << endl;
00545 
00546   assert(solver.Factored());
00547   assert(solver.Transpose()==Transpose);
00548   assert(solver.ReciprocalConditionEstimated());
00549   assert(solver.Solved());
00550 
00551   if (solver.SolutionErrorsEstimated()) {
00552     if (verbose) {
00553       cout << "\n\nFERR[0] = "<< solver.FERR()[0] << endl;
00554       cout << "\n\nBERR[0] = "<< solver.BERR()[0] << endl<< endl;
00555     }
00556   }
00557   
00558   double * resid = new double[NRHS];
00559   OK = Residual(N, NRHS, A1, LDA1, solver.Transpose(), solver.X(), solver.LDX(), B1, LDB1, resid);
00560   if (verbose) {
00561     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00562  /*
00563     if (solver.A_Equilibrated()) {
00564       double * R = solver.R();
00565       double * C = solver.C();
00566       for (i=0; i<solver.M(); i++) 
00567       cout << "R[" << i <<"] = "<< R[i] << endl;
00568       for (i=0; i<solver.N(); i++) 
00569       cout << "C[" << i <<"] = "<< C[i] << endl;
00570     }
00571  */
00572     cout << "\n\nResiduals using factorization to solve" << endl;
00573     for (i=0; i<NRHS; i++)
00574       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00575     cout  << endl;
00576   }
00577 
00578 
00579   ierr = solver.Invert();
00580   assert(ierr>-1);
00581 
00582   assert(solver.Inverted());
00583   assert(!solver.Factored());
00584   assert(solver.Transpose()==Transpose);
00585 
00586   
00587   Epetra_SerialDenseMatrix RHS1(Copy, B1, LDB1, N, NRHS);
00588   Epetra_SerialDenseMatrix LHS1(Copy, X1, LDX1, N, NRHS);
00589   assert(solver.SetVectors(LHS1, RHS1)==0);
00590   assert(!solver.Solved());
00591 
00592   assert(solver.Solve()>-1);
00593    
00594   
00595 
00596   OK = Residual(N, NRHS, A1, LDA1, solver.Transpose(), solver.X(), solver.LDX(), B1, LDB1, resid);
00597 
00598   if (verbose) {
00599     if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
00600     cout << "Residuals using inverse to solve" << endl;
00601     for (i=0; i<NRHS; i++)
00602       cout << "Residual[" << i <<"] = "<< resid[i] << endl;
00603     cout  << endl;
00604   }
00605   delete [] resid;
00606   
00607 
00608   return(0);
00609 }
00610 
00611  void GenerateHilbert(double *A, int LDA, int N) {
00612    for (int j=0; j<N; j++)
00613      for (int i=0; i<N; i++)
00614        A[i+j*LDA] = 1.0/((double)(i+j+1));
00615    return;
00616  }
00617 
00618 bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
00619          double * X, int LDX, double * B, int LDB, double * resid) {
00620 
00621   Epetra_BLAS Blas;
00622   char Transa = 'N';
00623   if (Transpose) Transa = 'T';
00624   Blas.GEMM(Transa, 'N', N, NRHS, N, -1.0, A, LDA,
00625       X, LDX, 1.0, B, LDB);
00626   bool OK = true;
00627   for (int i=0; i<NRHS; i++) {
00628     resid[i] = Blas.NRM2(N, B+i*LDB);
00629     if (resid[i]>1.0E-7) OK = false;
00630   }
00631 
00632   return(OK);
00633 }
00634 
00635 
00636 //=========================================================================
00637 // test matrix operator= (copy & view)
00638 int matrixAssignment(bool verbose, bool debug) {
00639   int ierr = 0;
00640   int returnierr = 0;
00641   if(verbose) printHeading("Testing matrix operator=");
00642 
00643   // each section is in its own block so we can reuse variable names
00644   // lhs = left hand side, rhs = right hand side
00645   
00646   {
00647     // copy->copy (more space needed)
00648     // orig and dup should have same signature
00649     // modifying orig or dup should have no effect on the other
00650     if(verbose) cout << "Checking copy->copy (new alloc)" << endl;
00651     Epetra_SerialDenseMatrix lhs(2,2);
00652     double* rand1 = getRandArray(25);
00653     Epetra_SerialDenseMatrix rhs(Copy, rand1, 5, 5, 5);
00654     if(debug) {
00655       cout << "before assignment:" << endl;
00656       printMat("rhs",rhs);
00657       printMat("lhs",lhs);
00658     }
00659     lhs = rhs;
00660     if(debug) {
00661       cout << "after assignment:" << endl;
00662       printMat("rhs",rhs);
00663       printMat("lhs",lhs);
00664     }
00665     EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
00666     EPETRA_TEST_ERR(!seperateData(rhs,lhs), ierr);
00667     delete[] rand1;
00668   }
00669   returnierr += ierr;
00670   if(ierr == 0)
00671     if(verbose) cout << "Checked OK." << endl;
00672   ierr = 0;
00673   {
00674     // copy->copy (have enough space)
00675     // orig and dup should have same signature
00676     // modifying orig or dup should have no effect on the other
00677     if(verbose) cout << "\nChecking copy->copy (no alloc)" << endl;
00678     double* rand1 = getRandArray(25);
00679     double* rand2 = getRandArray(20);
00680     Epetra_SerialDenseMatrix lhs(Copy, rand1, 5, 5, 5);
00681     Epetra_SerialDenseMatrix rhs(Copy, rand2, 4, 4, 5);
00682     double* origA = lhs.A();
00683     int origLDA = lhs.LDA();
00684     if(debug) {
00685       cout << "before assignment:" << endl;
00686       printMat("rhs",rhs);
00687       printMat("lhs",lhs);
00688     }
00689     lhs = rhs;
00690     if(debug) {
00691       cout << "after assignment:" << endl;
00692       printMat("rhs",rhs);
00693       printMat("lhs",lhs);
00694     }
00695     // in this case, instead of doing a "normal" LDA test in identSig,
00696     // we do our own. Since we had enough space already, A and LDA should
00697     // not have been changed by the assignment. (The extra parameter to
00698     // identicalSignatures tells it not to test LDA).
00699     EPETRA_TEST_ERR((lhs.A() != origA) || (lhs.LDA() != origLDA), ierr);
00700     EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs,false), ierr);
00701     EPETRA_TEST_ERR(!seperateData(rhs,lhs), ierr);
00702   }
00703   returnierr += ierr;
00704   if(ierr == 0)
00705     if(verbose) cout << "Checked OK." << endl;
00706   ierr = 0;
00707   {
00708     // view->copy
00709     // orig and dup should have same signature
00710     // modifying orig or dup should have no effect on the other
00711     if(verbose) cout << "\nChecking view->copy" << endl;
00712     double* rand1 = getRandArray(25);
00713     double* rand2 = getRandArray(64);
00714     Epetra_SerialDenseMatrix lhs(View, rand1, 5, 5, 5);
00715     Epetra_SerialDenseMatrix rhs(Copy, rand2, 8, 8, 8);
00716     if(debug) {
00717       cout << "before assignment:" << endl;
00718       printMat("rhs",rhs);
00719       printMat("lhs",lhs);
00720     }
00721     lhs = rhs;
00722     if(debug) {
00723       cout << "after assignment:" << endl;
00724       printMat("rhs",rhs);
00725       printMat("lhs",lhs);
00726     }
00727     EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
00728     EPETRA_TEST_ERR(!seperateData(rhs,lhs), ierr);
00729     delete[] rand1;
00730     delete[] rand2;
00731   }
00732   returnierr += ierr;
00733   if(ierr == 0)
00734     if(verbose) cout << "Checked OK." << endl;
00735   ierr = 0;
00736   {
00737     // copy->view
00738     // orig and dup should have same signature
00739     // modifying orig or dup should change the other
00740     if(verbose) cout << "\nChecking copy->view" << endl;
00741     double* rand1 = getRandArray(10);
00742     Epetra_SerialDenseMatrix lhs(4,4);
00743     Epetra_SerialDenseMatrix rhs(View, rand1, 2, 2, 5);
00744     if(debug) {
00745       cout << "before assignment:" << endl;
00746       printMat("rhs",rhs);
00747       printMat("lhs",lhs);
00748     }
00749     lhs = rhs;
00750     if(debug) {
00751       cout << "after assignment:" << endl;
00752       printMat("rhs",rhs);
00753       printMat("lhs",lhs);
00754     }
00755     EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
00756     EPETRA_TEST_ERR(seperateData(rhs,lhs), ierr);
00757     delete[] rand1;
00758   }
00759   returnierr += ierr;
00760   if(ierr == 0)
00761     if(verbose) cout << "Checked OK." << endl;
00762   ierr = 0; 
00763   {
00764     // view->view
00765     // orig and dup should have same signature
00766     // modifying orig or dup should change the other
00767     if(verbose) cout << "\nChecking view->view" << endl;
00768     double* rand1 = getRandArray(9);
00769     double* rand2 = getRandArray(18);
00770     Epetra_SerialDenseMatrix lhs(View, rand1, 3, 3, 3);
00771     Epetra_SerialDenseMatrix rhs(View, rand2, 3, 3, 6);
00772     if(debug) {
00773       cout << "before assignment:" << endl;
00774       printMat("rhs",rhs);
00775       printMat("lhs",lhs);
00776     }
00777     lhs = rhs;
00778     if(debug) {
00779       cout << "after assignment:" << endl;
00780       printMat("rhs",rhs);
00781       printMat("lhs",lhs);
00782     }
00783     EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
00784     EPETRA_TEST_ERR(seperateData(rhs,lhs), ierr);
00785     delete[] rand1;
00786     delete[] rand2;
00787   }
00788   returnierr += ierr;
00789   if(ierr == 0)
00790     if(verbose) cout << "Checked OK." << endl;
00791   ierr = 0;
00792 
00793   return(returnierr);
00794 }
00795 
00796 //=========================================================================
00797 // test matrix copy constructor (copy & view)
00798 int matrixCpyCtr(bool verbose, bool debug) {
00799   const int m1rows = 5;
00800   const int m1cols = 4;
00801   const int m2rows = 2;
00802   const int m2cols = 6;
00803 
00804   int ierr = 0;
00805   int returnierr = 0;
00806   if(verbose) printHeading("Testing matrix copy constructors");
00807 
00808   if(verbose) cout << "checking copy constructor (view)" << endl;
00809   double* m1rand = getRandArray(m1rows * m1cols);
00810   if(debug) printArray(m1rand, m1rows * m1cols);
00811   Epetra_SerialDenseMatrix m1(View, m1rand, m1rows, m1rows, m1cols);
00812   if(debug) {
00813     cout << "original matrix:" << endl;
00814     printMat("m1",m1);
00815   }
00816   Epetra_SerialDenseMatrix m1clone(m1);
00817   if(debug) {
00818     cout << "clone matrix:" << endl;
00819     printMat("m1clone",m1clone);
00820   }
00821   if(verbose) cout << "making sure signatures match" << endl;
00822   EPETRA_TEST_ERR(!identicalSignatures(m1, m1clone), ierr);
00823   delete[] m1rand;
00824   returnierr += ierr;
00825   if(ierr == 0)
00826     if(verbose) cout << "Checked OK." << endl;
00827   ierr = 0;
00828   
00829   if(verbose) cout << "\nchecking copy constructor (copy)" << endl;
00830   double* m2rand = getRandArray(m2rows * m2cols);
00831   if(debug) printArray(m2rand, m2rows * m2cols);
00832   Epetra_SerialDenseMatrix m2(Copy, m2rand, m2rows, m2rows, m2cols);
00833   if(debug) {
00834     cout << "original matrix:" << endl;
00835     printMat("m2",m2);
00836   }
00837   Epetra_SerialDenseMatrix m2clone(m2);
00838   if(debug) {
00839     cout << "clone matrix:" << endl;
00840     printMat("m2clone",m2clone);
00841   }
00842   if(verbose) cout << "checking that signatures match" << endl;
00843   EPETRA_TEST_ERR(!identicalSignatures(m2, m2clone), ierr);
00844   returnierr += ierr;
00845   if(ierr == 0)
00846     if(verbose) cout << "Checked OK." << endl;
00847   ierr = 0;
00848 
00849   if(verbose) cout << "\nmodifying entry in m2, m2clone should be unchanged" << endl;
00850   EPETRA_TEST_ERR(!seperateData(m2, m2clone), ierr);
00851   if(debug) {
00852     printArray(m2rand, m2rows * m2cols);
00853     cout << "orig:" << endl;
00854     printMat("m2",m2);
00855     cout << "clone:" << endl;
00856     printMat("m2clone",m2clone);
00857   }
00858   delete[] m2rand;
00859   returnierr += ierr;
00860   if(ierr == 0)
00861     if(verbose) cout << "Checked OK." << endl;
00862   ierr = 0;
00863   
00864   return(returnierr);
00865 }
00866 
00867 //=========================================================================
00868 // prints section heading with spacers/formatting
00869 void printHeading(const char* heading) {
00870   cout << "\n==================================================================\n";
00871   cout << heading << endl;
00872   cout << "==================================================================\n";
00873 }
00874 
00875 //=========================================================================
00876 // prints SerialDenseMatrix/Vector with formatting
00877 void printMat(const char* name, Epetra_SerialDenseMatrix& matrix) {
00878   //cout << "--------------------" << endl;
00879   cout << "*** " << name << " ***" << endl;
00880   cout << matrix;
00881   //cout << "--------------------" << endl;
00882 }
00883 
00884 //=========================================================================
00885 // prints double* array with formatting
00886 void printArray(double* array, int length) {
00887   cout << "user array (size " << length << "): ";
00888   for(int i = 0; i < length; i++)
00889     cout << array[i] << "  ";
00890   cout << endl;
00891 }
00892 
00893 //=========================================================================
00894 // returns a double* array of a given length, with random values on interval (-1,1).
00895 // this is the same generator used in SerialDenseMatrix
00896 double* getRandArray(int length) {
00897   const double a = 16807.0;
00898   const double BigInt = 2147483647.0;
00899   const double DbleOne = 1.0;
00900   const double DbleTwo = 2.0;
00901   double seed = rand();
00902 
00903   double* array = new double[length];
00904 
00905   for(int i = 0; i < length; i++) {
00906     seed = fmod(a * seed, BigInt);
00907     array[i] = DbleTwo * (seed / BigInt) - DbleOne;
00908   }
00909 
00910   return(array);
00911 }
00912 
00913 //=========================================================================
00914 // checks the signatures of two matrices
00915 bool identicalSignatures(Epetra_SerialDenseMatrix& a, Epetra_SerialDenseMatrix& b, bool testLDA) {
00916 
00917   if((a.M()  != b.M()  )|| // check properties first
00918      (a.N()  != b.N()  )||
00919      (a.CV() != b.CV() ))
00920     return(false);
00921 
00922   if(testLDA == true)      // if we are coming from op= c->c #2 (have enough space)
00923     if(a.LDA() != b.LDA()) // then we don't check LDA (but we do check it in the test function)
00924       return(false);
00925 
00926   if(a.CV() == View) { // if we're still here, we need to check the data
00927     if(a.A() != b.A()) // for a view, this just means checking the pointers
00928       return(false);   // for a copy, this means checking each element
00929   }
00930   else { // CV == Copy
00931     const int m = a.M();
00932     const int n = a.N();
00933     for(int i = 0; i < m; i++)
00934       for(int j = 0; j < n; j++) {
00935         if(a(i,j) != b(i,j))
00936           return(false);
00937       }
00938   }
00939 
00940   return(true); // if we're still here, signatures are identical
00941 }
00942 
00943 //=========================================================================
00944 // checks if two matrices are independent or not
00945 bool seperateData(Epetra_SerialDenseMatrix& a, Epetra_SerialDenseMatrix& b) {
00946   bool seperate;
00947 
00948   int r = EPETRA_MIN(a.M(),b.M()) / 2; // ensures (r,c) is valid
00949   int c = EPETRA_MIN(a.N(),b.N()) / 2; // in both matrices
00950 
00951   double orig_a = a(r,c);
00952   double new_value = a(r,c) + 1;
00953   if(b(r,c) == new_value) // there's a chance b could be independent, but
00954     new_value++;          // already have new_value in (r,c).
00955   
00956   a(r,c) = new_value;
00957   if(b(r,c) == new_value)
00958     seperate = false;
00959   else
00960     seperate = true;
00961 
00962   a(r,c) = orig_a; // undo change we made to a
00963 
00964   return(seperate);
00965 }

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