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  }
734  returnierr += ierr;
735  if(ierr == 0)
736  if(verbose) cout << "Checked OK." << endl;
737  ierr = 0;
738  {
739  // view->copy
740  // orig and dup should have same signature
741  // modifying orig or dup should have no effect on the other
742  if(verbose) cout << "\nChecking view->copy" << endl;
743  double* rand1 = getRandArray(25);
744  double* rand2 = getRandArray(64);
745  Epetra_SerialDenseMatrix lhs(View, rand1, 5, 5, 5);
746  Epetra_SerialDenseMatrix rhs(Copy, rand2, 8, 8, 8);
747  if(debug) {
748  cout << "before assignment:" << endl;
749  printMat("rhs",rhs);
750  printMat("lhs",lhs);
751  }
752  lhs = rhs;
753  if(debug) {
754  cout << "after assignment:" << endl;
755  printMat("rhs",rhs);
756  printMat("lhs",lhs);
757  }
758  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
759  EPETRA_TEST_ERR(!seperateData(rhs,lhs), ierr);
760  delete[] rand1;
761  delete[] rand2;
762  }
763  returnierr += ierr;
764  if(ierr == 0)
765  if(verbose) cout << "Checked OK." << endl;
766  ierr = 0;
767  {
768  // copy->view
769  // orig and dup should have same signature
770  // modifying orig or dup should change the other
771  if(verbose) cout << "\nChecking copy->view" << endl;
772  double* rand1 = getRandArray(10);
773  Epetra_SerialDenseMatrix lhs(4,4);
774  Epetra_SerialDenseMatrix rhs(View, rand1, 2, 2, 5);
775  if(debug) {
776  cout << "before assignment:" << endl;
777  printMat("rhs",rhs);
778  printMat("lhs",lhs);
779  }
780  lhs = rhs;
781  if(debug) {
782  cout << "after assignment:" << endl;
783  printMat("rhs",rhs);
784  printMat("lhs",lhs);
785  }
786  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
787  EPETRA_TEST_ERR(seperateData(rhs,lhs), ierr);
788  delete[] rand1;
789  }
790  returnierr += ierr;
791  if(ierr == 0)
792  if(verbose) cout << "Checked OK." << endl;
793  ierr = 0;
794  {
795  // view->view
796  // orig and dup should have same signature
797  // modifying orig or dup should change the other
798  if(verbose) cout << "\nChecking view->view" << endl;
799  double* rand1 = getRandArray(9);
800  double* rand2 = getRandArray(18);
801  Epetra_SerialDenseMatrix lhs(View, rand1, 3, 3, 3);
802  Epetra_SerialDenseMatrix rhs(View, rand2, 3, 3, 6);
803  if(debug) {
804  cout << "before assignment:" << endl;
805  printMat("rhs",rhs);
806  printMat("lhs",lhs);
807  }
808  lhs = rhs;
809  if(debug) {
810  cout << "after assignment:" << endl;
811  printMat("rhs",rhs);
812  printMat("lhs",lhs);
813  }
814  EPETRA_TEST_ERR(!identicalSignatures(rhs,lhs), ierr);
815  EPETRA_TEST_ERR(seperateData(rhs,lhs), ierr);
816  delete[] rand1;
817  delete[] rand2;
818  }
819  returnierr += ierr;
820  if(ierr == 0)
821  if(verbose) cout << "Checked OK." << endl;
822  ierr = 0;
823 
824  return(returnierr);
825 }
826 
827 //=========================================================================
828 // test matrix copy constructor (copy & view)
829 int matrixCpyCtr(bool verbose, bool debug) {
830  const int m1rows = 5;
831  const int m1cols = 4;
832  const int m2rows = 2;
833  const int m2cols = 6;
834 
835  int ierr = 0;
836  int returnierr = 0;
837  if(verbose) printHeading("Testing matrix copy constructors");
838 
839  if(verbose) cout << "checking copy constructor (view)" << endl;
840  double* m1rand = getRandArray(m1rows * m1cols);
841  if(debug) printArray(m1rand, m1rows * m1cols);
842  Epetra_SerialDenseMatrix m1(View, m1rand, m1rows, m1rows, m1cols);
843  if(debug) {
844  cout << "original matrix:" << endl;
845  printMat("m1",m1);
846  }
847  Epetra_SerialDenseMatrix m1clone(m1);
848  if(debug) {
849  cout << "clone matrix:" << endl;
850  printMat("m1clone",m1clone);
851  }
852  if(verbose) cout << "making sure signatures match" << endl;
853  EPETRA_TEST_ERR(!identicalSignatures(m1, m1clone), ierr);
854  delete[] m1rand;
855  returnierr += ierr;
856  if(ierr == 0)
857  if(verbose) cout << "Checked OK." << endl;
858  ierr = 0;
859 
860  if(verbose) cout << "\nchecking copy constructor (copy)" << endl;
861  double* m2rand = getRandArray(m2rows * m2cols);
862  if(debug) printArray(m2rand, m2rows * m2cols);
863  Epetra_SerialDenseMatrix m2(Copy, m2rand, m2rows, m2rows, m2cols);
864  if(debug) {
865  cout << "original matrix:" << endl;
866  printMat("m2",m2);
867  }
868  Epetra_SerialDenseMatrix m2clone(m2);
869  if(debug) {
870  cout << "clone matrix:" << endl;
871  printMat("m2clone",m2clone);
872  }
873  if(verbose) cout << "checking that signatures match" << endl;
874  EPETRA_TEST_ERR(!identicalSignatures(m2, m2clone), ierr);
875  returnierr += ierr;
876  if(ierr == 0)
877  if(verbose) cout << "Checked OK." << endl;
878  ierr = 0;
879 
880  if(verbose) cout << "\nmodifying entry in m2, m2clone should be unchanged" << endl;
881  EPETRA_TEST_ERR(!seperateData(m2, m2clone), ierr);
882  if(debug) {
883  printArray(m2rand, m2rows * m2cols);
884  cout << "orig:" << endl;
885  printMat("m2",m2);
886  cout << "clone:" << endl;
887  printMat("m2clone",m2clone);
888  }
889  delete[] m2rand;
890  returnierr += ierr;
891  if(ierr == 0)
892  if(verbose) cout << "Checked OK." << endl;
893  ierr = 0;
894 
895  return(returnierr);
896 }
897 
898 //=========================================================================
899 // prints section heading with spacers/formatting
900 void printHeading(const char* heading) {
901  cout << "\n==================================================================\n";
902  cout << heading << endl;
903  cout << "==================================================================\n";
904 }
905 
906 //=========================================================================
907 // prints SerialDenseMatrix/Vector with formatting
908 void printMat(const char* name, Epetra_SerialDenseMatrix& matrix) {
909  //cout << "--------------------" << endl;
910  cout << "*** " << name << " ***" << endl;
911  cout << matrix;
912  //cout << "--------------------" << endl;
913 }
914 
915 //=========================================================================
916 // prints double* array with formatting
917 void printArray(double* array, int length) {
918  cout << "user array (size " << length << "): ";
919  for(int i = 0; i < length; i++)
920  cout << array[i] << " ";
921  cout << endl;
922 }
923 
924 //=========================================================================
925 // returns a double* array of a given length, with random values on interval (-1,1).
926 // this is the same generator used in SerialDenseMatrix
927 double* getRandArray(int length) {
928  const double a = 16807.0;
929  const double BigInt = 2147483647.0;
930  const double DbleOne = 1.0;
931  const double DbleTwo = 2.0;
932  double seed = rand();
933 
934  double* array = new double[length];
935 
936  for(int i = 0; i < length; i++) {
937  seed = fmod(a * seed, BigInt);
938  array[i] = DbleTwo * (seed / BigInt) - DbleOne;
939  }
940 
941  return(array);
942 }
943 
944 //=========================================================================
945 // checks the signatures of two matrices
947 
948  if((a.M() != b.M() )|| // check properties first
949  (a.N() != b.N() )||
950  (a.CV() != b.CV() ))
951  return(false);
952 
953  if(testLDA == true) // if we are coming from op= c->c #2 (have enough space)
954  if(a.LDA() != b.LDA()) // then we don't check LDA (but we do check it in the test function)
955  return(false);
956 
957  if(a.CV() == View) { // if we're still here, we need to check the data
958  if(a.A() != b.A()) // for a view, this just means checking the pointers
959  return(false); // for a copy, this means checking each element
960  }
961  else { // CV == Copy
962  const int m = a.M();
963  const int n = a.N();
964  for(int i = 0; i < m; i++)
965  for(int j = 0; j < n; j++) {
966  if(a(i,j) != b(i,j))
967  return(false);
968  }
969  }
970 
971  return(true); // if we're still here, signatures are identical
972 }
973 
974 //=========================================================================
975 // checks if two matrices are independent or not
977  bool seperate;
978 
979  int r = EPETRA_MIN(a.M(),b.M()) / 2; // ensures (r,c) is valid
980  int c = EPETRA_MIN(a.N(),b.N()) / 2; // in both matrices
981 
982  double orig_a = a(r,c);
983  double new_value = a(r,c) + 1;
984  if(b(r,c) == new_value) // there's a chance b could be independent, but
985  new_value++; // already have new_value in (r,c).
986 
987  a(r,c) = new_value;
988  if(b(r,c) == new_value)
989  seperate = false;
990  else
991  seperate = true;
992 
993  a(r,c) = orig_a; // undo change we made to a
994 
995  return(seperate);
996 }
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.