Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
58int 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
64void GenerateHilbert(double *A, int LDA, int N);
65
66bool Residual( int N, int NRHS, double * A, int LDA, bool Transpose,
67 double * X, int LDX, double * B, int LDB, double * resid);
68
69int matrixCpyCtr(bool verbose, bool debug);
70int matrixAssignment(bool verbose, bool debug);
71void printHeading(const char* heading);
72double* getRandArray(int length);
73void printMat(const char* name, Epetra_SerialDenseMatrix& matrix);
74void printArray(double* array, int length);
77
78
79int 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
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
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*/
501return ierr ;
502}
503
504int 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
649bool 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)
669int 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;
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);
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)
831int 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
902void printHeading(const char* heading) {
903 cout << "\n==================================================================\n";
904 cout << heading << endl;
905 cout << "==================================================================\n";
906}
907
908//=========================================================================
909// prints SerialDenseMatrix/Vector with formatting
910void 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
919void 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
929double* 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}
#define EPETRA_MIN(x, y)
std::string Epetra_Version()
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition Epetra_BLAS.h:70
float NRM2(const int N, const float *X, const int INCX=1) const
Epetra_BLAS norm function (SNRM2).
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)
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
double Flops() const
Returns the number of floating point operations with this multi-vector.
Epetra_Flops: The Epetra Floating Point Operations Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Epetra_MpiComm: The Epetra MPI Communication Class.
int MyPID() const
Return my process ID.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_SerialDenseMatrix: A class for constructing and using real double precision general dense matr...
double * A() const
Returns pointer to the this matrix.
int LDA() const
Returns the leading dimension of the this matrix.
virtual double NormInf() const
Computes the Infinity-Norm of the this matrix.
int Scale(double ScalarA)
Inplace scalar-matrix product A = a A.
int Shape(int NumRows, int NumCols)
Set dimensions of a Epetra_SerialDenseMatrix object; init values to zero.
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.
int M() const
Returns row dimension of system.
int N() const
Returns column dimension of system.
Epetra_DataAccess CV() const
Returns the data access mode of the this matrix.
virtual double NormOne() const
Computes the 1-Norm of the this matrix.
Epetra_SerialDenseSolver: A class for solving dense linear problems.
int M() const
Returns row dimension of system.
bool ReciprocalConditionEstimated()
Returns true if the condition number of the this matrix has been computed (value available via Recipr...
int LDB() const
Returns the leading dimension of the RHS.
virtual int Solve(void)
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
int LDX() const
Returns the leading dimension of the solution.
double AMAX() const
Returns the absolute value of the largest entry of the this matrix (returns -1 if not yet computed).
double COLCND() const
Ratio of smallest to largest column scale factors for the this matrix (returns -1 if not yet computed...
bool SolutionRefined()
Returns true if the current set of vectors has been refined.
void SolveToRefinedSolution(bool Flag)
Causes all solves to compute solution to best ability using iterative refinement.
int SetMatrix(Epetra_SerialDenseMatrix &A)
Sets the pointers for coefficient matrix.
bool Transpose()
Returns true if transpose of this matrix has and will be used.
virtual int ReciprocalConditionEstimate(double &Value)
Returns the reciprocal of the 1-norm condition number of the this matrix.
double * X() const
Returns pointer to current solution.
bool SolutionErrorsEstimated()
Returns true if forward and backward error estimated have been computed (available via FERR() and BER...
int N() const
Returns column dimension of system.
double ROWCND() const
Ratio of smallest to largest row scale factors for the this matrix (returns -1 if not yet computed).
void SolveWithTranspose(bool Flag)
If Flag is true, causes all subsequent function calls to work with the transpose of this matrix,...
bool A_Equilibrated()
Returns true if factor is equilibrated (factor available via AF() and LDAF()).
double * FERR() const
Returns a pointer to the forward error estimates computed by LAPACK.
int NRHS() const
Returns the number of current right hand sides and solution vectors.
double ANORM() const
Returns the 1-Norm of the this matrix (returns -1 if not yet computed).
double RCOND() const
Returns the reciprocal of the condition number of the this matrix (returns -1 if not yet computed).
virtual int Factor(void)
Computes the in-place LU factorization of the matrix using the LAPACK routine DGETRF.
int LDA() const
Returns the leading dimension of the this matrix.
bool Factored()
Returns true if matrix is factored (factor available via AF() and LDAF()).
bool Inverted()
Returns true if matrix inverse has been computed (inverse available via AF() and LDAF()).
virtual int Invert(void)
Inverts the this matrix.
bool B_Equilibrated()
Returns true if RHS is equilibrated (RHS available via B() and LDB()).
double * BERR() const
Returns a pointer to the backward error estimates computed by LAPACK.
bool Solved()
Returns true if the current set of vectors has been solved.
void FactorWithEquilibration(bool Flag)
Causes equilibration to be called just before the matrix factorization as part of the call to Factor.
int SetVectors(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &B)
Sets the pointers for left and right hand side vector(s).
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
int Size(int Length_in)
Set length of a Epetra_SerialDenseVector object; init values to zero.
double Norm1() const
Compute 1-norm of each vector in multi-vector.
double NormInf() const
Compute Inf-norm of each vector in multi-vector.
double Norm2() const
Compute 2-norm of each vector in multi-vector.
double Dot(const Epetra_SerialDenseVector &x) const
Compute 1-norm of each vector in multi-vector.
int Length() const
Returns length of vector.
Epetra_Time: The Epetra Timing Class.
Definition Epetra_Time.h:75
double ElapsedTime(void) const
Epetra_Time elapsed time function.
#define EPETRA_TEST_ERR(a, b)
void printMat(const char *name, Epetra_SerialDenseMatrix &matrix)
int main(int argc, char *argv[])
int matrixAssignment(bool verbose, bool debug)
int matrixCpyCtr(bool verbose, bool debug)
double * getRandArray(int length)
void GenerateHilbert(double *A, int LDA, int N)
void printArray(double *array, int length)
bool seperateData(Epetra_SerialDenseMatrix &a, Epetra_SerialDenseMatrix &b)
bool identicalSignatures(Epetra_SerialDenseMatrix &a, Epetra_SerialDenseMatrix &b, bool testLDA=true)
int check(Epetra_SerialDenseSolver &solver, double *A1, int LDA, int N1, int NRHS1, double OneNorm1, double *B1, int LDB1, double *X1, int LDX1, bool Transpose, bool verbose)
void printHeading(const char *heading)
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)