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