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

Generated on Wed May 12 21:41:04 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7