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