Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/SerialDense/cxx_main.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Epetra: Linear Algebra Services Package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 
43 #include "Epetra_Map.h"
44 #include "Epetra_Time.h"
48 #ifdef EPETRA_MPI
49 #include "Epetra_MpiComm.h"
50 #include <mpi.h>
51 #endif
52 #include "Epetra_SerialComm.h"
53 #include "../epetra_test_err.h"
54 #include "Epetra_Version.h"
55 
56 // prototypes
57 
58 int check(Epetra_SerialDenseSolver & solver, double * A1, int LDA,
59  int N1, int NRHS1, double OneNorm1,
60  double * B1, int LDB1,
61  double * X1, int LDX1,
62  bool Transpose, bool verbose);
63 
64 void GenerateHilbert(double *A, int LDA, int N);
65 
66 bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
67  double * X, int LDX, double * B, int LDB, double * resid);
68 
69 int matrixCpyCtr(bool verbose, bool debug);
70 int matrixAssignment(bool verbose, bool debug);
71 void printHeading(const char* heading);
72 double* getRandArray(int length);
73 void printMat(const char* name, Epetra_SerialDenseMatrix& matrix);
74 void printArray(double* array, int length);
77 
78 
79 int main(int argc, char *argv[])
80 {
81  int ierr = 0, i, j, k;
82  bool debug = false;
83 
84 #ifdef EPETRA_MPI
85  MPI_Init(&argc,&argv);
86  Epetra_MpiComm Comm( MPI_COMM_WORLD );
87 #else
88  Epetra_SerialComm Comm;
89 #endif
90 
91  bool verbose = false;
92 
93  // Check if we should print results to standard out
94  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
95 
96  if (verbose && Comm.MyPID()==0)
97  cout << Epetra_Version() << endl << endl;
98 
99  int rank = Comm.MyPID();
100  // char tmp;
101  // if (rank==0) cout << "Press any key to continue..."<< endl;
102  // if (rank==0) cin >> tmp;
103  // Comm.Barrier();
104 
105  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
106  if (verbose) cout << Comm <<endl;
107 
108  // bool verbose1 = verbose;
109 
110  // Redefine verbose to only print on PE 0
111  if (verbose && rank!=0) verbose = false;
112 
113  int N = 20;
114  int NRHS = 4;
115  double * A = new double[N*N];
116  double * A1 = new double[N*N];
117  double * X = new double[(N+1)*NRHS];
118  double * X1 = new double[(N+1)*NRHS];
119  int LDX = N+1;
120  int LDX1 = N+1;
121  double * B = new double[N*NRHS];
122  double * B1 = new double[N*NRHS];
123  int LDB = N;
124  int LDB1 = N;
125 
126  int LDA = N;
127  int LDA1 = LDA;
128  double OneNorm1;
129  bool Transpose = false;
130 
132  Epetra_SerialDenseMatrix * Matrix;
133  for (int kk=0; kk<2; kk++) {
134  for (i=1; i<=N; i++) {
135  GenerateHilbert(A, LDA, i);
136  OneNorm1 = 0.0;
137  for (j=1; j<=i; j++) OneNorm1 += 1.0/((double) j); // 1-Norm = 1 + 1/2 + ...+1/n
138 
139  if (kk==0) {
140  Matrix = new Epetra_SerialDenseMatrix(View, A, LDA, i, i);
141  LDA1 = LDA;
142  }
143  else {
144  Matrix = new Epetra_SerialDenseMatrix(Copy, A, LDA, i, i);
145  LDA1 = i;
146  }
147 
148  GenerateHilbert(A1, LDA1, i);
149 
150  if (kk==1) {
151  solver.FactorWithEquilibration(true);
152  solver.SolveWithTranspose(true);
153  Transpose = true;
154  solver.SolveToRefinedSolution(true);
155  }
156 
157  for (k=0; k<NRHS; k++)
158  for (j=0; j<i; j++) {
159  B[j+k*LDB] = 1.0/((double) (k+3)*(j+3));
160  B1[j+k*LDB1] = B[j+k*LDB1];
161  }
162  Epetra_SerialDenseMatrix Epetra_B(View, B, LDB, i, NRHS);
163  Epetra_SerialDenseMatrix Epetra_X(View, X, LDX, i, NRHS);
164 
165  solver.SetMatrix(*Matrix);
166  solver.SetVectors(Epetra_X, Epetra_B);
167 
168  ierr = check(solver, A1, LDA1, i, NRHS, OneNorm1, B1, LDB1, X1, LDX1, Transpose, verbose);
169  assert (ierr>-1);
170  delete Matrix;
171  if (ierr!=0) {
172  if (verbose) cout << "Factorization failed due to bad conditioning. This is normal if RCOND is small."
173  << endl;
174  break;
175  }
176  }
177  }
178 
179  delete [] A;
180  delete [] A1;
181  delete [] X;
182  delete [] X1;
183  delete [] B;
184  delete [] B1;
185 
187  // Now test norms and scaling functions
189 
191  double ScalarA = 2.0;
192 
193  int DM = 10;
194  int DN = 8;
195  D.Shape(DM, DN);
196  for (j=0; j<DN; j++)
197  for (i=0; i<DM; i++) D[j][i] = (double) (1+i+j*DM) ;
198 
199  //cout << D << endl;
200 
201  double NormInfD_ref = (double)(DM*(DN*(DN+1))/2);
202  double NormOneD_ref = (double)((DM*DN*(DM*DN+1))/2 - (DM*(DN-1)*(DM*(DN-1)+1))/2 );
203 
204  double NormInfD = D.NormInf();
205  double NormOneD = D.NormOne();
206 
207  if (verbose) {
208  cout << " *** Before scaling *** " << endl
209  << " Computed one-norm of test matrix = " << NormOneD << endl
210  << " Expected one-norm = " << NormOneD_ref << endl
211  << " Computed inf-norm of test matrix = " << NormInfD << endl
212  << " Expected inf-norm = " << NormInfD_ref << endl;
213  }
214  D.Scale(ScalarA); // Scale entire D matrix by this value
215  NormInfD = D.NormInf();
216  NormOneD = D.NormOne();
217  if (verbose) {
218  cout << " *** After scaling *** " << endl
219  << " Computed one-norm of test matrix = " << NormOneD << endl
220  << " Expected one-norm = " << NormOneD_ref*ScalarA << endl
221  << " Computed inf-norm of test matrix = " << NormInfD << endl
222  << " Expected inf-norm = " << NormInfD_ref*ScalarA << endl;
223  }
224 
225 
227  // Now test that A.Multiply(false, x, y) produces the same result
228  // as y.Multiply('N','N', 1.0, A, x, 0.0).
230 
231  N = 10;
232  int M = 10;
233  LDA = N;
234  Epetra_SerialDenseMatrix smallA(N, M, false);
235  Epetra_SerialDenseMatrix x(N, 1, false);
236  Epetra_SerialDenseMatrix y1(N, 1, false);
237  Epetra_SerialDenseMatrix y2(N, 1, false);
238 
239  for(i=0; i<N; ++i) {
240  for(j=0; j<M; ++j) {
241  smallA(i,j) = 1.0*i+2.0*j+1.0;
242  }
243  x(i,0) = 1.0;
244  y1(i,0) = 0.0;
245  y2(i,0) = 0.0;
246  }
247 
248  //quick check of operator==
249  if (x == y1) {
250  if (verbose) cout << "err in Epetra_SerialDenseMatrix::operator==, "
251  << "erroneously returned true." << std::endl;
252  return(-1);
253  }
254 
255  //quick check of operator!=
256  if (x != x) {
257  if (verbose) cout << "err in Epetra_SerialDenseMatrix::operator==, "
258  << "erroneously returned true." << std::endl;
259  return(-1);
260  }
261 
262  int err1 = smallA.Multiply(false, x, y1);
263  int err2 = y2.Multiply('N','N', 1.0, smallA, x, 0.0);
264  if (err1 != 0 || err2 != 0) {
265  if (verbose) cout << "err in Epetra_SerialDenseMatrix::Multiply"<<endl;
266  return(err1+err2);
267  }
268 
269  for(i=0; i<N; ++i) {
270  if (y1(i,0) != y2(i,0)) {
271  if (verbose) cout << "different versions of Multiply don't match."<<endl;
272  return(-99);
273  }
274  }
275 
277  // Now test for larger system, both correctness and performance.
279 
280 
281  N = 2000;
282  NRHS = 5;
283  LDA = N;
284  LDB = N;
285  LDX = N;
286 
287  if (verbose) cout << "\n\nComputing factor of an " << N << " x " << N << " general matrix...Please wait.\n\n" << endl;
288 
289  // Define A and X
290 
291  A = new double[LDA*N];
292  X = new double[LDB*NRHS];
293 
294  for (j=0; j<N; j++) {
295  for (k=0; k<NRHS; k++) X[j+k*LDX] = 1.0/((double) (j+5+k));
296  for (i=0; i<N; i++) {
297  if (i==((j+2)%N)) A[i+j*LDA] = 100.0 + i;
298  else A[i+j*LDA] = -11.0/((double) (i+5)*(j+2));
299  }
300  }
301 
302  // Define Epetra_SerialDenseMatrix object
303 
304  Epetra_SerialDenseMatrix BigMatrix(Copy, A, LDA, N, N);
305  Epetra_SerialDenseMatrix OrigBigMatrix(View, A, LDA, N, N);
306 
307  Epetra_SerialDenseSolver BigSolver;
308  BigSolver.FactorWithEquilibration(true);
309  BigSolver.SetMatrix(BigMatrix);
310 
311  // Time factorization
312 
313  Epetra_Flops counter;
314  BigSolver.SetFlopCounter(counter);
315  Epetra_Time Timer(Comm);
316  double tstart = Timer.ElapsedTime();
317  ierr = BigSolver.Factor();
318  if (ierr!=0 && verbose) cout << "Error in factorization = "<<ierr<< endl;
319  assert(ierr==0);
320  double time = Timer.ElapsedTime() - tstart;
321 
322  double FLOPS = counter.Flops();
323  double MFLOPS = FLOPS/time/1000000.0;
324  if (verbose) cout << "MFLOPS for Factorization = " << MFLOPS << endl;
325 
326  // Define Left hand side and right hand side
327  Epetra_SerialDenseMatrix LHS(View, X, LDX, N, NRHS);
329  RHS.Shape(N,NRHS); // Allocate RHS
330 
331  // Compute RHS from A and X
332 
333  Epetra_Flops RHS_counter;
334  RHS.SetFlopCounter(RHS_counter);
335  tstart = Timer.ElapsedTime();
336  RHS.Multiply('N', 'N', 1.0, OrigBigMatrix, LHS, 0.0);
337  time = Timer.ElapsedTime() - tstart;
338 
339  Epetra_SerialDenseMatrix OrigRHS = RHS;
340 
341  FLOPS = RHS_counter.Flops();
342  MFLOPS = FLOPS/time/1000000.0;
343  if (verbose) cout << "MFLOPS to build RHS (NRHS = " << NRHS <<") = " << MFLOPS << endl;
344 
345  // Set LHS and RHS and solve
346  BigSolver.SetVectors(LHS, RHS);
347 
348  tstart = Timer.ElapsedTime();
349  ierr = BigSolver.Solve();
350  if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
351  else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
352  assert(ierr>=0);
353  time = Timer.ElapsedTime() - tstart;
354 
355  FLOPS = BigSolver.Flops();
356  MFLOPS = FLOPS/time/1000000.0;
357  if (verbose) cout << "MFLOPS for Solve (NRHS = " << NRHS <<") = " << MFLOPS << endl;
358 
359  double * resid = new double[NRHS];
360  bool OK = Residual(N, NRHS, A, LDA, BigSolver.Transpose(), BigSolver.X(), BigSolver.LDX(),
361  OrigRHS.A(), OrigRHS.LDA(), resid);
362 
363  if (verbose) {
364  if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
365  for (i=0; i<NRHS; i++)
366  cout << "Residual[" << i <<"] = "<< resid[i] << endl;
367  cout << endl;
368  }
369 
370  // Solve again using the Epetra_SerialDenseVector class for LHS and RHS
371 
374  X2.Size(BigMatrix.N());
375  B2.Size(BigMatrix.M());
376  int length = BigMatrix.N();
377  {for (int kk=0; kk<length; kk++) X2[kk] = ((double ) kk)/ ((double) length);} // Define entries of X2
378 
379  RHS_counter.ResetFlops();
380  B2.SetFlopCounter(RHS_counter);
381  tstart = Timer.ElapsedTime();
382  B2.Multiply('N', 'N', 1.0, OrigBigMatrix, X2, 0.0); // Define B2 = A*X2
383  time = Timer.ElapsedTime() - tstart;
384 
385  Epetra_SerialDenseVector OrigB2 = B2;
386 
387  FLOPS = RHS_counter.Flops();
388  MFLOPS = FLOPS/time/1000000.0;
389  if (verbose) cout << "MFLOPS to build single RHS = " << MFLOPS << endl;
390 
391  // Set LHS and RHS and solve
392  BigSolver.SetVectors(X2, B2);
393 
394  tstart = Timer.ElapsedTime();
395  ierr = BigSolver.Solve();
396  time = Timer.ElapsedTime() - tstart;
397  if (ierr==1 && verbose) cout << "LAPACK guidelines suggest this matrix might benefit from equilibration." << endl;
398  else if (ierr!=0 && verbose) cout << "Error in solve = "<<ierr<< endl;
399  assert(ierr>=0);
400 
401  FLOPS = counter.Flops();
402  MFLOPS = FLOPS/time/1000000.0;
403  if (verbose) cout << "MFLOPS to solve single RHS = " << MFLOPS << endl;
404 
405  OK = Residual(N, 1, A, LDA, BigSolver.Transpose(), BigSolver.X(), BigSolver.LDX(), OrigB2.A(),
406  OrigB2.LDA(), resid);
407 
408  if (verbose) {
409  if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
410  cout << "Residual = "<< resid[0] << endl;
411  }
412  delete [] resid;
413  delete [] A;
414  delete [] X;
415 
417  // Now test default constructor and index operators
419 
420  N = 5;
421  Epetra_SerialDenseMatrix C; // Implicit call to default constructor, should not need to call destructor
422  C.Shape(5,5); // Make it 5 by 5
423  double * C1 = new double[N*N];
424  GenerateHilbert(C1, N, N); // Generate Hilber matrix
425 
426  C1[1+2*N] = 1000.0; // Make matrix nonsymmetric
427 
428  // Fill values of C with Hilbert values
429  for (i=0; i<N; i++)
430  for (j=0; j<N; j++)
431  C(i,j) = C1[i+j*N];
432 
433  // Test if values are correctly written and read
434  for (i=0; i<N; i++)
435  for (j=0; j<N; j++) {
436  assert(C(i,j) == C1[i+j*N]);
437  assert(C(i,j) == C[j][i]);
438  }
439 
440  if (verbose)
441  cout << "Default constructor and index operator check OK. Values of Hilbert matrix = "
442  << endl << C << endl
443  << "Values should be 1/(i+j+1), except value (1,2) should be 1000" << endl;
444 
445  delete [] C1;
446 
447  // now test sized/shaped constructor
448  Epetra_SerialDenseMatrix shapedMatrix(10, 12);
449  assert(shapedMatrix.M() == 10);
450  assert(shapedMatrix.N() == 12);
451  for(i = 0; i < 10; i++)
452  for(j = 0; j < 12; j++)
453  assert(shapedMatrix(i, j) == 0.0);
454  Epetra_SerialDenseVector sizedVector(20);
455  assert(sizedVector.Length() == 20);
456  for(i = 0; i < 20; i++)
457  assert(sizedVector(i) == 0.0);
458  if (verbose)
459  cout << "Shaped/sized constructors check OK." << endl;
460 
461  // test Copy/View mode in op= and cpy ctr
462  int temperr = 0;
463  temperr = matrixAssignment(verbose, debug);
464  if(verbose && temperr == 0)
465  cout << "Operator = checked OK." << endl;
466  EPETRA_TEST_ERR(temperr, ierr);
467  temperr = matrixCpyCtr(verbose, debug);
468  if(verbose && temperr == 0)
469  cout << "Copy ctr checked OK." << endl;
470  EPETRA_TEST_ERR(temperr, ierr);
471 
472  // Test some vector methods
473 
475  v1[0] = 1.0;
476  v1[1] = 3.0;
477  v1[2] = 2.0;
478 
480  v2[0] = 2.0;
481  v2[1] = 1.0;
482  v2[2] = -2.0;
483 
484  temperr = 0;
485  if (v1.Norm1()!=6.0) temperr++;
486  if (fabs(sqrt(14.0)-v1.Norm2())>1.0e-6) temperr++;
487  if (v1.NormInf()!=3.0) temperr++;
488  if(verbose && temperr == 0)
489  cout << "Vector Norms checked OK." << endl;
490  temperr = 0;
491  if (v1.Dot(v2)!=1.0) temperr++;
492  if(verbose && temperr == 0)
493  cout << "Vector Dot product checked OK." << endl;
494 
495 #ifdef EPETRA_MPI
496  MPI_Finalize() ;
497 #endif
498 
499 /* end main
500 */
501 return ierr ;
502 }
503 
504 int check(Epetra_SerialDenseSolver &solver, double * A1, int LDA1,
505  int N1, int NRHS1, double OneNorm1,
506  double * B1, int LDB1,
507  double * X1, int LDX1,
508  bool Transpose, bool verbose) {
509 
510  int i;
511  bool OK;
512  // Test query functions
513 
514  int M= solver.M();
515  if (verbose) cout << "\n\nNumber of Rows = " << M << endl<< endl;
516  assert(M==N1);
517 
518  int N= solver.N();
519  if (verbose) cout << "\n\nNumber of Equations = " << N << endl<< endl;
520  assert(N==N1);
521 
522  int LDA = solver.LDA();
523  if (verbose) cout << "\n\nLDA = " << LDA << endl<< endl;
524  assert(LDA==LDA1);
525 
526  int LDB = solver.LDB();
527  if (verbose) cout << "\n\nLDB = " << LDB << endl<< endl;
528  assert(LDB==LDB1);
529 
530  int LDX = solver.LDX();
531  if (verbose) cout << "\n\nLDX = " << LDX << endl<< endl;
532  assert(LDX==LDX1);
533 
534  int NRHS = solver.NRHS();
535  if (verbose) cout << "\n\nNRHS = " << NRHS << endl<< endl;
536  assert(NRHS==NRHS1);
537 
538  assert(solver.ANORM()==-1.0);
539  assert(solver.RCOND()==-1.0);
540  if (!solver.A_Equilibrated() && !solver.B_Equilibrated()) {
541  assert(solver.ROWCND()==-1.0);
542  assert(solver.COLCND()==-1.0);
543  assert(solver.AMAX()==-1.0);
544  }
545 
546 
547  // Other binary tests
548 
549  assert(!solver.Factored());
550  assert(solver.Transpose()==Transpose);
551  assert(!solver.SolutionErrorsEstimated());
552  assert(!solver.Inverted());
553  assert(!solver.ReciprocalConditionEstimated());
554  assert(!solver.Solved());
555 
556  assert(!solver.SolutionRefined());
557 
558 
559  int ierr = solver.Factor();
560  assert(ierr>-1);
561  if (ierr!=0) return(ierr); // Factorization failed due to poor conditioning.
562  double rcond;
563  ierr = solver.ReciprocalConditionEstimate(rcond);
564  assert(ierr==0);
565  if (verbose) {
566 
567  double rcond1 = 1.0/exp(3.5*((double)N));
568  if (N==1) rcond1 = 1.0;
569  cout << "\n\nRCOND = "<< rcond << " should be approx = "
570  << rcond1 << endl << endl;
571  }
572 
573  ierr = solver.Solve();
574  assert(ierr>-1);
575  if (ierr!=0 && verbose) cout << "LAPACK rules suggest system should be equilibrated." << endl;
576 
577  assert(solver.Factored());
578  assert(solver.Transpose()==Transpose);
579  assert(solver.ReciprocalConditionEstimated());
580  assert(solver.Solved());
581 
582  if (solver.SolutionErrorsEstimated()) {
583  if (verbose) {
584  cout << "\n\nFERR[0] = "<< solver.FERR()[0] << endl;
585  cout << "\n\nBERR[0] = "<< solver.BERR()[0] << endl<< endl;
586  }
587  }
588 
589  double * resid = new double[NRHS];
590  OK = Residual(N, NRHS, A1, LDA1, solver.Transpose(), solver.X(), solver.LDX(), B1, LDB1, resid);
591  if (verbose) {
592  if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
593  /*
594  if (solver.A_Equilibrated()) {
595  double * R = solver.R();
596  double * C = solver.C();
597  for (i=0; i<solver.M(); i++)
598  cout << "R[" << i <<"] = "<< R[i] << endl;
599  for (i=0; i<solver.N(); i++)
600  cout << "C[" << i <<"] = "<< C[i] << endl;
601  }
602  */
603  cout << "\n\nResiduals using factorization to solve" << endl;
604  for (i=0; i<NRHS; i++)
605  cout << "Residual[" << i <<"] = "<< resid[i] << endl;
606  cout << endl;
607  }
608 
609 
610  ierr = solver.Invert();
611  assert(ierr>-1);
612 
613  assert(solver.Inverted());
614  assert(!solver.Factored());
615  assert(solver.Transpose()==Transpose);
616 
617 
618  Epetra_SerialDenseMatrix RHS1(Copy, B1, LDB1, N, NRHS);
619  Epetra_SerialDenseMatrix LHS1(Copy, X1, LDX1, N, NRHS);
620  assert(solver.SetVectors(LHS1, RHS1)==0);
621  assert(!solver.Solved());
622 
623  assert(solver.Solve()>-1);
624 
625 
626 
627  OK = Residual(N, NRHS, A1, LDA1, solver.Transpose(), solver.X(), solver.LDX(), B1, LDB1, resid);
628 
629  if (verbose) {
630  if (!OK) cout << "************* Residual do not meet tolerance *************" << endl;
631  cout << "Residuals using inverse to solve" << endl;
632  for (i=0; i<NRHS; i++)
633  cout << "Residual[" << i <<"] = "<< resid[i] << endl;
634  cout << endl;
635  }
636  delete [] resid;
637 
638 
639  return(0);
640 }
641 
642  void GenerateHilbert(double *A, int LDA, int N) {
643  for (int j=0; j<N; j++)
644  for (int i=0; i<N; i++)
645  A[i+j*LDA] = 1.0/((double)(i+j+1));
646  return;
647  }
648 
649 bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
650  double * X, int LDX, double * B, int LDB, double * resid) {
651 
652  Epetra_BLAS Blas;
653  char Transa = 'N';
654  if (Transpose) Transa = 'T';
655  Blas.GEMM(Transa, 'N', N, NRHS, N, -1.0, A, LDA,
656  X, LDX, 1.0, B, LDB);
657  bool OK = true;
658  for (int i=0; i<NRHS; i++) {
659  resid[i] = Blas.NRM2(N, B+i*LDB);
660  if (resid[i]>1.0E-7) OK = false;
661  }
662 
663  return(OK);
664 }
665 
666 
667 //=========================================================================
668 // test matrix operator= (copy & view)
669 int matrixAssignment(bool verbose, bool debug) {
670  int ierr = 0;
671  int returnierr = 0;
672  if(verbose) printHeading("Testing matrix operator=");
673 
674  // each section is in its own block so we can reuse variable names
675  // lhs = left hand side, rhs = right hand side
676 
677  {
678  // copy->copy (more space needed)
679  // orig and dup should have same signature
680  // modifying orig or dup should have no effect on the other
681  if(verbose) cout << "Checking copy->copy (new alloc)" << endl;
682  Epetra_SerialDenseMatrix lhs(2,2);
683  double* rand1 = getRandArray(25);
684  Epetra_SerialDenseMatrix rhs(Copy, rand1, 5, 5, 5);
685  if(debug) {
686  cout << "before assignment:" << endl;
687  printMat("rhs",rhs);
688  printMat("lhs",lhs);
689  }
690  lhs = rhs;
691  if(debug) {
692  cout << "after assignment:" << endl;
693  printMat("rhs",rhs);
694  printMat("lhs",lhs);
695  }
696  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
697  EPETRA_TEST_ERR(!seperateData(rhs,lhs), ierr);
698  delete[] rand1;
699  }
700  returnierr += ierr;
701  if(ierr == 0)
702  if(verbose) cout << "Checked OK." << endl;
703  ierr = 0;
704  {
705  // copy->copy (have enough space)
706  // orig and dup should have same signature
707  // modifying orig or dup should have no effect on the other
708  if(verbose) cout << "\nChecking copy->copy (no alloc)" << endl;
709  double* rand1 = getRandArray(25);
710  double* rand2 = getRandArray(20);
711  Epetra_SerialDenseMatrix lhs(Copy, rand1, 5, 5, 5);
712  Epetra_SerialDenseMatrix rhs(Copy, rand2, 4, 4, 5);
713  double* origA = lhs.A();
714  int origLDA = lhs.LDA();
715  if(debug) {
716  cout << "before assignment:" << endl;
717  printMat("rhs",rhs);
718  printMat("lhs",lhs);
719  }
720  lhs = rhs;
721  if(debug) {
722  cout << "after assignment:" << endl;
723  printMat("rhs",rhs);
724  printMat("lhs",lhs);
725  }
726  // in this case, instead of doing a "normal" LDA test in identSig,
727  // we do our own. Since we had enough space already, A and LDA should
728  // not have been changed by the assignment. (The extra parameter to
729  // identicalSignatures tells it not to test LDA).
730  EPETRA_TEST_ERR((lhs.A() != origA) || (lhs.LDA() != origLDA), ierr);
731  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs,false), ierr);
732  EPETRA_TEST_ERR(!seperateData(rhs,lhs), ierr);
733  delete [] rand1;
734  delete [] rand2;
735  }
736  returnierr += ierr;
737  if(ierr == 0)
738  if(verbose) cout << "Checked OK." << endl;
739  ierr = 0;
740  {
741  // view->copy
742  // orig and dup should have same signature
743  // modifying orig or dup should have no effect on the other
744  if(verbose) cout << "\nChecking view->copy" << endl;
745  double* rand1 = getRandArray(25);
746  double* rand2 = getRandArray(64);
747  Epetra_SerialDenseMatrix lhs(View, rand1, 5, 5, 5);
748  Epetra_SerialDenseMatrix rhs(Copy, rand2, 8, 8, 8);
749  if(debug) {
750  cout << "before assignment:" << endl;
751  printMat("rhs",rhs);
752  printMat("lhs",lhs);
753  }
754  lhs = rhs;
755  if(debug) {
756  cout << "after assignment:" << endl;
757  printMat("rhs",rhs);
758  printMat("lhs",lhs);
759  }
760  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
761  EPETRA_TEST_ERR(!seperateData(rhs,lhs), ierr);
762  delete[] rand1;
763  delete[] rand2;
764  }
765  returnierr += ierr;
766  if(ierr == 0)
767  if(verbose) cout << "Checked OK." << endl;
768  ierr = 0;
769  {
770  // copy->view
771  // orig and dup should have same signature
772  // modifying orig or dup should change the other
773  if(verbose) cout << "\nChecking copy->view" << endl;
774  double* rand1 = getRandArray(10);
775  Epetra_SerialDenseMatrix lhs(4,4);
776  Epetra_SerialDenseMatrix rhs(View, rand1, 2, 2, 5);
777  if(debug) {
778  cout << "before assignment:" << endl;
779  printMat("rhs",rhs);
780  printMat("lhs",lhs);
781  }
782  lhs = rhs;
783  if(debug) {
784  cout << "after assignment:" << endl;
785  printMat("rhs",rhs);
786  printMat("lhs",lhs);
787  }
788  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
789  EPETRA_TEST_ERR(seperateData(rhs,lhs), ierr);
790  delete[] rand1;
791  }
792  returnierr += ierr;
793  if(ierr == 0)
794  if(verbose) cout << "Checked OK." << endl;
795  ierr = 0;
796  {
797  // view->view
798  // orig and dup should have same signature
799  // modifying orig or dup should change the other
800  if(verbose) cout << "\nChecking view->view" << endl;
801  double* rand1 = getRandArray(9);
802  double* rand2 = getRandArray(18);
803  Epetra_SerialDenseMatrix lhs(View, rand1, 3, 3, 3);
804  Epetra_SerialDenseMatrix rhs(View, rand2, 3, 3, 6);
805  if(debug) {
806  cout << "before assignment:" << endl;
807  printMat("rhs",rhs);
808  printMat("lhs",lhs);
809  }
810  lhs = rhs;
811  if(debug) {
812  cout << "after assignment:" << endl;
813  printMat("rhs",rhs);
814  printMat("lhs",lhs);
815  }
816  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
817  EPETRA_TEST_ERR(seperateData(rhs,lhs), ierr);
818  delete[] rand1;
819  delete[] rand2;
820  }
821  returnierr += ierr;
822  if(ierr == 0)
823  if(verbose) cout << "Checked OK." << endl;
824  ierr = 0;
825 
826  return(returnierr);
827 }
828 
829 //=========================================================================
830 // test matrix copy constructor (copy & view)
831 int matrixCpyCtr(bool verbose, bool debug) {
832  const int m1rows = 5;
833  const int m1cols = 4;
834  const int m2rows = 2;
835  const int m2cols = 6;
836 
837  int ierr = 0;
838  int returnierr = 0;
839  if(verbose) printHeading("Testing matrix copy constructors");
840 
841  if(verbose) cout << "checking copy constructor (view)" << endl;
842  double* m1rand = getRandArray(m1rows * m1cols);
843  if(debug) printArray(m1rand, m1rows * m1cols);
844  Epetra_SerialDenseMatrix m1(View, m1rand, m1rows, m1rows, m1cols);
845  if(debug) {
846  cout << "original matrix:" << endl;
847  printMat("m1",m1);
848  }
849  Epetra_SerialDenseMatrix m1clone(m1);
850  if(debug) {
851  cout << "clone matrix:" << endl;
852  printMat("m1clone",m1clone);
853  }
854  if(verbose) cout << "making sure signatures match" << endl;
855  EPETRA_TEST_ERR(!identicalSignatures(m1, m1clone), ierr);
856  delete[] m1rand;
857  returnierr += ierr;
858  if(ierr == 0)
859  if(verbose) cout << "Checked OK." << endl;
860  ierr = 0;
861 
862  if(verbose) cout << "\nchecking copy constructor (copy)" << endl;
863  double* m2rand = getRandArray(m2rows * m2cols);
864  if(debug) printArray(m2rand, m2rows * m2cols);
865  Epetra_SerialDenseMatrix m2(Copy, m2rand, m2rows, m2rows, m2cols);
866  if(debug) {
867  cout << "original matrix:" << endl;
868  printMat("m2",m2);
869  }
870  Epetra_SerialDenseMatrix m2clone(m2);
871  if(debug) {
872  cout << "clone matrix:" << endl;
873  printMat("m2clone",m2clone);
874  }
875  if(verbose) cout << "checking that signatures match" << endl;
876  EPETRA_TEST_ERR(!identicalSignatures(m2, m2clone), ierr);
877  returnierr += ierr;
878  if(ierr == 0)
879  if(verbose) cout << "Checked OK." << endl;
880  ierr = 0;
881 
882  if(verbose) cout << "\nmodifying entry in m2, m2clone should be unchanged" << endl;
883  EPETRA_TEST_ERR(!seperateData(m2, m2clone), ierr);
884  if(debug) {
885  printArray(m2rand, m2rows * m2cols);
886  cout << "orig:" << endl;
887  printMat("m2",m2);
888  cout << "clone:" << endl;
889  printMat("m2clone",m2clone);
890  }
891  delete[] m2rand;
892  returnierr += ierr;
893  if(ierr == 0)
894  if(verbose) cout << "Checked OK." << endl;
895  ierr = 0;
896 
897  return(returnierr);
898 }
899 
900 //=========================================================================
901 // prints section heading with spacers/formatting
902 void printHeading(const char* heading) {
903  cout << "\n==================================================================\n";
904  cout << heading << endl;
905  cout << "==================================================================\n";
906 }
907 
908 //=========================================================================
909 // prints SerialDenseMatrix/Vector with formatting
910 void printMat(const char* name, Epetra_SerialDenseMatrix& matrix) {
911  //cout << "--------------------" << endl;
912  cout << "*** " << name << " ***" << endl;
913  cout << matrix;
914  //cout << "--------------------" << endl;
915 }
916 
917 //=========================================================================
918 // prints double* array with formatting
919 void printArray(double* array, int length) {
920  cout << "user array (size " << length << "): ";
921  for(int i = 0; i < length; i++)
922  cout << array[i] << " ";
923  cout << endl;
924 }
925 
926 //=========================================================================
927 // returns a double* array of a given length, with random values on interval (-1,1).
928 // this is the same generator used in SerialDenseMatrix
929 double* getRandArray(int length) {
930  const double a = 16807.0;
931  const double BigInt = 2147483647.0;
932  const double DbleOne = 1.0;
933  const double DbleTwo = 2.0;
934  double seed = rand();
935 
936  double* array = new double[length];
937 
938  for(int i = 0; i < length; i++) {
939  seed = fmod(a * seed, BigInt);
940  array[i] = DbleTwo * (seed / BigInt) - DbleOne;
941  }
942 
943  return(array);
944 }
945 
946 //=========================================================================
947 // checks the signatures of two matrices
949 
950  if((a.M() != b.M() )|| // check properties first
951  (a.N() != b.N() )||
952  (a.CV() != b.CV() ))
953  return(false);
954 
955  if(testLDA == true) // if we are coming from op= c->c #2 (have enough space)
956  if(a.LDA() != b.LDA()) // then we don't check LDA (but we do check it in the test function)
957  return(false);
958 
959  if(a.CV() == View) { // if we're still here, we need to check the data
960  if(a.A() != b.A()) // for a view, this just means checking the pointers
961  return(false); // for a copy, this means checking each element
962  }
963  else { // CV == Copy
964  const int m = a.M();
965  const int n = a.N();
966  for(int i = 0; i < m; i++)
967  for(int j = 0; j < n; j++) {
968  if(a(i,j) != b(i,j))
969  return(false);
970  }
971  }
972 
973  return(true); // if we're still here, signatures are identical
974 }
975 
976 //=========================================================================
977 // checks if two matrices are independent or not
979  bool seperate;
980 
981  int r = EPETRA_MIN(a.M(),b.M()) / 2; // ensures (r,c) is valid
982  int c = EPETRA_MIN(a.N(),b.N()) / 2; // in both matrices
983 
984  double orig_a = a(r,c);
985  double new_value = a(r,c) + 1;
986  if(b(r,c) == new_value) // there's a chance b could be independent, but
987  new_value++; // already have new_value in (r,c).
988 
989  a(r,c) = new_value;
990  if(b(r,c) == new_value)
991  seperate = false;
992  else
993  seperate = true;
994 
995  a(r,c) = orig_a; // undo change we made to a
996 
997  return(seperate);
998 }
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
virtual double NormInf() const
Computes the Infinity-Norm of the this matrix.
bool Transpose()
Returns true if transpose of this matrix has and will be used.
void GEMM(const char TRANSA, const char TRANSB, const int M, const int N, const int K, const float ALPHA, const float *A, const int LDA, const float *B, const int LDB, const float BETA, float *C, const int LDC) const
Epetra_BLAS matrix-matrix multiply function (SGEMM)
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
double ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed)...
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)
int matrixAssignment(bool verbose, bool debug)
double Norm2() const
Compute 2-norm of each vector in multi-vector.
float NRM2(const int N, const float *X, const int INCX=1) const
Epetra_BLAS norm function (SNRM2).
Definition: Epetra_BLAS.cpp:81
#define EPETRA_TEST_ERR(a, b)
bool B_Equilibrated()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
double ElapsedTime(void) const
Epetra_Time elapsed time function.
bool SolutionRefined()
Returns true if the current set of vectors has been refined.
double COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
double * BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
Epetra_DataAccess CV() const
Returns the data access mode of the this matrix.
bool seperateData(Epetra_IntSerialDenseMatrix &a, Epetra_IntSerialDenseMatrix &b)
virtual double NormOne() const
Computes the 1-Norm of the this matrix.
double NormInf() const
Compute Inf-norm of each vector in multi-vector.
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
#define EPETRA_MIN(x, y)
int M() const
Returns row dimension of system.
void SolveWithTranspose(bool Flag)
If Flag is true, causes all subsequent function calls to work with the transpose of this matrix...
bool SolutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
std::string Epetra_Version()
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
double ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
void printHeading(const char *heading)
double RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed)...
double * A() const
Returns pointer to the this matrix.
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)
int N() const
Returns column dimension of system.
int matrixCpyCtr(bool verbose, bool debug)
int Length() const
Returns length of vector.
void printArray(int *array, int length)
int LDB() const
Returns the leading dimension of the RHS.
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
double * X() const
Returns pointer to current solution.
int LDA() const
Returns the leading dimension of the this matrix.
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
double Norm1() const
Compute 1-norm of each vector in multi-vector.
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
virtual int Invert(void)
Inverts the this matrix.
void SolveToRefinedSolution(bool Flag)
Causes all solves to compute solution to best ability using iterative refinement. ...
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
int SetMatrix(Epetra_SerialDenseMatrix &A)
Sets the pointers for coefficient matrix.
int LDA() const
Returns the leading dimension of the this matrix.
int LDX() const
Returns the leading dimension of the solution.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Definition: Epetra_Flops.h:74
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_SerialDenseMatrix &A, const Epetra_SerialDenseMatrix &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
bool A_Equilibrated()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
double Flops() const
Returns the number of floating point operations with this multi-vector.
void FactorWithEquilibration(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor...
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
double AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed)...
double * FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
bool Solved()
Returns true if the current set of vectors has been solved.
int main(int argc, char *argv[])
double Dot(const Epetra_SerialDenseVector &x) const
Compute 1-norm of each vector in multi-vector.
int N() const
Returns column dimension of system.
int NRHS() const
Returns the number of current right hand sides and solution vectors.
int * getRandArray(int length)
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
bool identicalSignatures(Epetra_IntSerialDenseMatrix &a, Epetra_IntSerialDenseMatrix &b, bool testLDA=true)
virtual int ReciprocalConditionEstimate(double &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
int n
virtual int Factor(void)
Computes the in-place LU factorization of the matrix using the LAPACK routine DGETRF.
int M() const
Returns row dimension of system.
void GenerateHilbert(double *A, int LDA, int N)
Epetra_SerialDenseSolver: A class for solving dense linear problems.