Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/BasicPerfTest/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 #define EPETRA_HAVE_JADMATRIX
44 #define EPETRA_VERY_SHORT_PERFTEST
45 #define EPETRA_HAVE_STATICPROFILE
46 #include "Epetra_Map.h"
47 #include "Epetra_LocalMap.h"
48 #include "Epetra_BlockMap.h"
49 #include "Epetra_Time.h"
50 #include "Epetra_CrsMatrix.h"
51 #include "Epetra_VbrMatrix.h"
52 #include "Epetra_Vector.h"
53 #include "Epetra_IntVector.h"
54 #include "Epetra_MultiVector.h"
55 #include "Epetra_IntSerialDenseVector.h"
56 #include "Epetra_SerialDenseVector.h"
57 #include "Epetra_Flops.h"
58 #ifdef EPETRA_MPI
59 #include "Epetra_MpiComm.h"
60 #include "mpi.h"
61 #else
62 #include "Epetra_SerialComm.h"
63 #endif
64 #include "../epetra_test_err.h"
65 #include "Epetra_Version.h"
66 #ifdef EPETRA_HAVE_JADMATRIX
67 #include "Epetra_JadMatrix.h"
68 #endif
69 
70 // prototypes
71 
72 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
73  int * xoff, int * yoff,
74  const Epetra_Comm &comm, bool verbose, bool summary,
75  Epetra_Map *& map,
76  Epetra_CrsMatrix *& A,
77  Epetra_Vector *& b,
78  Epetra_Vector *& bt,
79  Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly);
80 
81 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
82  int * xoff, int * yoff, int nrhs,
83  const Epetra_Comm &comm, bool verbose, bool summary,
84  Epetra_Map *& map,
85  Epetra_CrsMatrix *& A,
86  Epetra_MultiVector *& b,
87  Epetra_MultiVector *& bt,
88  Epetra_MultiVector *&xexact, bool StaticProfile, bool MakeLocalOnly);
89 
90 void GenerateVbrProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
91  int * xoff, int * yoff,
92  int nsizes, int * sizes,
93  const Epetra_Comm &comm, bool verbose, bool summary,
94  Epetra_BlockMap *& map,
95  Epetra_VbrMatrix *& A,
96  Epetra_Vector *& b,
97  Epetra_Vector *& bt,
98  Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly);
99 
100 void GenerateVbrProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
101  int * xoff, int * yoff,
102  int nsizes, int * sizes, int nrhs,
103  const Epetra_Comm &comm, bool verbose, bool summary,
104  Epetra_BlockMap *& map,
105  Epetra_VbrMatrix *& A,
106  Epetra_MultiVector *& b,
107  Epetra_MultiVector *& bt,
108  Epetra_MultiVector *&xexact, bool StaticProfile, bool MakeLocalOnly);
109 
110 void GenerateMyGlobalElements(int numNodesX, int numNodesY, int numProcsX, int numProcs,
111  int myPID, int * & myGlobalElements);
112 
113 void runMatrixTests(Epetra_CrsMatrix * A, Epetra_MultiVector * b, Epetra_MultiVector * bt,
114  Epetra_MultiVector * xexact, bool StaticProfile, bool verbose, bool summary);
115 #ifdef EPETRA_HAVE_JADMATRIX
116 void runJadMatrixTests(Epetra_JadMatrix * A, Epetra_MultiVector * b, Epetra_MultiVector * bt,
117  Epetra_MultiVector * xexact, bool StaticProfile, bool verbose, bool summary);
118 #endif
119 void runLUMatrixTests(Epetra_CrsMatrix * L, Epetra_MultiVector * bL, Epetra_MultiVector * btL, Epetra_MultiVector * xexactL,
120  Epetra_CrsMatrix * U, Epetra_MultiVector * bU, Epetra_MultiVector * btU, Epetra_MultiVector * xexactU,
121  bool StaticProfile, bool verbose, bool summary);
122 int main(int argc, char *argv[])
123 {
124  int ierr = 0;
125  double elapsed_time;
126  double total_flops;
127  double MFLOPs;
128 
129 
130 #ifdef EPETRA_MPI
131 
132  // Initialize MPI
133  MPI_Init(&argc,&argv);
134  Epetra_MpiComm comm( MPI_COMM_WORLD );
135 #else
136  Epetra_SerialComm comm;
137 #endif
138 
139  bool verbose = false;
140  bool summary = false;
141 
142  // Check if we should print verbose results to standard out
143  if (argc>6) if (argv[6][0]=='-' && argv[6][1]=='v') verbose = true;
144 
145  // Check if we should print verbose results to standard out
146  if (argc>6) if (argv[6][0]=='-' && argv[6][1]=='s') summary = true;
147 
148  if(argc < 6) {
149  cerr << "Usage: " << argv[0]
150  << " NumNodesX NumNodesY NumProcX NumProcY NumPoints [-v|-s]" << endl
151  << "where:" << endl
152  << "NumNodesX - Number of mesh nodes in X direction per processor" << endl
153  << "NumNodesY - Number of mesh nodes in Y direction per processor" << endl
154  << "NumProcX - Number of processors to use in X direction" << endl
155  << "NumProcY - Number of processors to use in Y direction" << endl
156  << "NumPoints - Number of points to use in stencil (5, 9 or 25 only)" << endl
157  << "-v|-s - (Optional) Run in verbose mode if -v present or summary mode if -s present" << endl
158  << " NOTES: NumProcX*NumProcY must equal the number of processors used to run the problem." << endl << endl
159  << " Serial example:" << endl
160  << argv[0] << " 16 12 1 1 25 -v" << endl
161  << " Run this program in verbose mode on 1 processor using a 16 X 12 grid with a 25 point stencil."<< endl <<endl
162  << " MPI example:" << endl
163  << "mpirun -np 32 " << argv[0] << " 10 12 4 8 9 -v" << endl
164  << " Run this program in verbose mode on 32 processors putting a 10 X 12 subgrid on each processor using 4 processors "<< endl
165  << " in the X direction and 8 in the Y direction. Total grid size is 40 points in X and 96 in Y with a 9 point stencil."<< endl
166  << endl;
167  return(1);
168 
169  }
170  //char tmp;
171  //if (comm.MyPID()==0) cout << "Press any key to continue..."<< endl;
172  //if (comm.MyPID()==0) cin >> tmp;
173  //comm.Barrier();
174 
175  comm.SetTracebackMode(0); // This should shut down any error traceback reporting
176  if (verbose && comm.MyPID()==0)
177  cout << Epetra_Version() << endl << endl;
178  if (summary && comm.MyPID()==0) {
179  if (comm.NumProc()==1)
180  cout << Epetra_Version() << endl << endl;
181  else
182  cout << endl << endl; // Print two blank line to keep output columns lined up
183  }
184 
185  if (verbose) cout << comm <<endl;
186 
187 
188  // Redefine verbose to only print on PE 0
189 
190  if (verbose && comm.MyPID()!=0) verbose = false;
191  if (summary && comm.MyPID()!=0) summary = false;
192 
193  int numNodesX = atoi(argv[1]);
194  int numNodesY = atoi(argv[2]);
195  int numProcsX = atoi(argv[3]);
196  int numProcsY = atoi(argv[4]);
197  int numPoints = atoi(argv[5]);
198 
199  if (verbose || (summary && comm.NumProc()==1)) {
200  cout << " Number of local nodes in X direction = " << numNodesX << endl
201  << " Number of local nodes in Y direction = " << numNodesY << endl
202  << " Number of global nodes in X direction = " << numNodesX*numProcsX << endl
203  << " Number of global nodes in Y direction = " << numNodesY*numProcsY << endl
204  << " Number of local nonzero entries = " << numNodesX*numNodesY*numPoints << endl
205  << " Number of global nonzero entries = " << numNodesX*numNodesY*numPoints*numProcsX*numProcsY << endl
206  << " Number of Processors in X direction = " << numProcsX << endl
207  << " Number of Processors in Y direction = " << numProcsY << endl
208  << " Number of Points in stencil = " << numPoints << endl << endl;
209  }
210  // Print blank line to keep output columns lined up
211  if (summary && comm.NumProc()>1)
212  cout << endl << endl << endl << endl << endl << endl << endl << endl<< endl << endl;
213 
214  if (numProcsX*numProcsY!=comm.NumProc()) {
215  cerr << "Number of processors = " << comm.NumProc() << endl
216  << " is not the product of " << numProcsX << " and " << numProcsY << endl << endl;
217  return(1);
218  }
219 
220  if (numPoints!=5 && numPoints!=9 && numPoints!=25) {
221  cerr << "Number of points specified = " << numPoints << endl
222  << " is not 5, 9, 25" << endl << endl;
223  return(1);
224  }
225 
226  if (numNodesX*numNodesY<=0) {
227  cerr << "Product of number of nodes is <= zero" << endl << endl;
228  return(1);
229  }
230 
231  Epetra_IntSerialDenseVector Xoff, XLoff, XUoff;
232  Epetra_IntSerialDenseVector Yoff, YLoff, YUoff;
233  if (numPoints==5) {
234 
235  // Generate a 5-point 2D Finite Difference matrix
236  Xoff.Size(5);
237  Yoff.Size(5);
238  Xoff[0] = -1; Xoff[1] = 1; Xoff[2] = 0; Xoff[3] = 0; Xoff[4] = 0;
239  Yoff[0] = 0; Yoff[1] = 0; Yoff[2] = 0; Yoff[3] = -1; Yoff[4] = 1;
240 
241  // Generate a 2-point 2D Lower triangular Finite Difference matrix
242  XLoff.Size(2);
243  YLoff.Size(2);
244  XLoff[0] = -1; XLoff[1] = 0;
245  YLoff[0] = 0; YLoff[1] = -1;
246 
247  // Generate a 3-point 2D upper triangular Finite Difference matrix
248  XUoff.Size(3);
249  YUoff.Size(3);
250  XUoff[0] = 0; XUoff[1] = 1; XUoff[2] = 0;
251  YUoff[0] = 0; YUoff[1] = 0; YUoff[2] = 1;
252  }
253  else if (numPoints==9) {
254  // Generate a 9-point 2D Finite Difference matrix
255  Xoff.Size(9);
256  Yoff.Size(9);
257  Xoff[0] = -1; Xoff[1] = 0; Xoff[2] = 1;
258  Yoff[0] = -1; Yoff[1] = -1; Yoff[2] = -1;
259  Xoff[3] = -1; Xoff[4] = 0; Xoff[5] = 1;
260  Yoff[3] = 0; Yoff[4] = 0; Yoff[5] = 0;
261  Xoff[6] = -1; Xoff[7] = 0; Xoff[8] = 1;
262  Yoff[6] = 1; Yoff[7] = 1; Yoff[8] = 1;
263 
264  // Generate a 5-point lower triangular 2D Finite Difference matrix
265  XLoff.Size(5);
266  YLoff.Size(5);
267  XLoff[0] = -1; XLoff[1] = 0; Xoff[2] = 1;
268  YLoff[0] = -1; YLoff[1] = -1; Yoff[2] = -1;
269  XLoff[3] = -1; XLoff[4] = 0;
270  YLoff[3] = 0; YLoff[4] = 0;
271 
272  // Generate a 4-point upper triangular 2D Finite Difference matrix
273  XUoff.Size(4);
274  YUoff.Size(4);
275  XUoff[0] = 1;
276  YUoff[0] = 0;
277  XUoff[1] = -1; XUoff[2] = 0; XUoff[3] = 1;
278  YUoff[1] = 1; YUoff[2] = 1; YUoff[3] = 1;
279 
280  }
281  else {
282  // Generate a 25-point 2D Finite Difference matrix
283  Xoff.Size(25);
284  Yoff.Size(25);
285  int xi = 0, yi = 0;
286  int xo = -2, yo = -2;
287  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
288  Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ;
289  xo = -2, yo++;
290  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
291  Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ;
292  xo = -2, yo++;
293  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
294  Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ;
295  xo = -2, yo++;
296  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
297  Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ;
298  xo = -2, yo++;
299  Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++; Xoff[xi++] = xo++;
300  Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ; Yoff[yi++] = yo ;
301 
302  // Generate a 13-point lower triangular 2D Finite Difference matrix
303  XLoff.Size(13);
304  YLoff.Size(13);
305  xi = 0, yi = 0;
306  xo = -2, yo = -2;
307  XLoff[xi++] = xo++; XLoff[xi++] = xo++; XLoff[xi++] = xo++; XLoff[xi++] = xo++; XLoff[xi++] = xo++;
308  YLoff[yi++] = yo ; YLoff[yi++] = yo ; YLoff[yi++] = yo ; YLoff[yi++] = yo ; YLoff[yi++] = yo ;
309  xo = -2, yo++;
310  XLoff[xi++] = xo++; XLoff[xi++] = xo++; XLoff[xi++] = xo++; XLoff[xi++] = xo++; XLoff[xi++] = xo++;
311  YLoff[yi++] = yo ; YLoff[yi++] = yo ; YLoff[yi++] = yo ; YLoff[yi++] = yo ; YLoff[yi++] = yo ;
312  xo = -2, yo++;
313  XLoff[xi++] = xo++; XLoff[xi++] = xo++; XLoff[xi++] = xo++;
314  YLoff[yi++] = yo ; YLoff[yi++] = yo ; YLoff[yi++] = yo ;
315 
316  // Generate a 13-point upper triangular 2D Finite Difference matrix
317  XUoff.Size(13);
318  YUoff.Size(13);
319  xi = 0, yi = 0;
320  xo = 0, yo = 0;
321  XUoff[xi++] = xo++; XUoff[xi++] = xo++; XUoff[xi++] = xo++;
322  YUoff[yi++] = yo ; YUoff[yi++] = yo ; YUoff[yi++] = yo ;
323  xo = -2, yo++;
324  XUoff[xi++] = xo++; XUoff[xi++] = xo++; XUoff[xi++] = xo++; XUoff[xi++] = xo++; XUoff[xi++] = xo++;
325  YUoff[yi++] = yo ; YUoff[yi++] = yo ; YUoff[yi++] = yo ; YUoff[yi++] = yo ; YUoff[yi++] = yo ;
326  xo = -2, yo++;
327  XUoff[xi++] = xo++; XUoff[xi++] = xo++; XUoff[xi++] = xo++; XUoff[xi++] = xo++; XUoff[xi++] = xo++;
328  YUoff[yi++] = yo ; YUoff[yi++] = yo ; YUoff[yi++] = yo ; YUoff[yi++] = yo ; YUoff[yi++] = yo ;
329 
330  }
331 
332  Epetra_Map * map;
333  Epetra_Map * mapL;
334  Epetra_Map * mapU;
335  Epetra_CrsMatrix * A;
336  Epetra_CrsMatrix * L;
337  Epetra_CrsMatrix * U;
338  Epetra_MultiVector * b;
339  Epetra_MultiVector * bt;
340  Epetra_MultiVector * xexact;
341  Epetra_MultiVector * bL;
342  Epetra_MultiVector * btL;
343  Epetra_MultiVector * xexactL;
344  Epetra_MultiVector * bU;
345  Epetra_MultiVector * btU;
346  Epetra_MultiVector * xexactU;
347  Epetra_SerialDenseVector resvec(0);
348 
349  //Timings
350  Epetra_Flops flopcounter;
351  Epetra_Time timer(comm);
352 
353 #ifdef EPETRA_VERY_SHORT_PERFTEST
354  int jstop = 1;
355 #elif EPETRA_SHORT_PERFTEST
356  int jstop = 1;
357 #else
358  int jstop = 2;
359 #endif
360  for (int j=0; j<jstop; j++) {
361  for (int k=1; k<17; k++) {
362 #ifdef EPETRA_VERY_SHORT_PERFTEST
363  if (k<3 || (k%4==0 && k<9)) {
364 #elif EPETRA_SHORT_PERFTEST
365  if (k<6 || k%4==0) {
366 #else
367  if (k<7 || k%2==0) {
368 #endif
369  int nrhs=k;
370  if (verbose) cout << "\n*************** Results for " << nrhs << " RHS with ";
371 
372  bool StaticProfile = (j!=0);
373  if (verbose) {
374  if (StaticProfile) cout << " static profile\n";
375  else cout << " dynamic profile\n";
376  }
377 
378  GenerateCrsProblem(numNodesX, numNodesY, numProcsX, numProcsY, numPoints,
379  Xoff.Values(), Yoff.Values(), nrhs, comm, verbose, summary,
380  map, A, b, bt, xexact, StaticProfile, false);
381 
382 
383 #ifdef EPETRA_HAVE_JADMATRIX
384 
385  timer.ResetStartTime();
386  Epetra_JadMatrix JA(*A);
387  elapsed_time = timer.ElapsedTime();
388  if (verbose) cout << "Time to create Jagged diagonal matrix = " << elapsed_time << endl;
389 
390  //cout << "A = " << *A << endl;
391  //cout << "JA = " << JA << endl;
392 
393  runJadMatrixTests(&JA, b, bt, xexact, StaticProfile, verbose, summary);
394 
395 #endif
396  runMatrixTests(A, b, bt, xexact, StaticProfile, verbose, summary);
397 
398  delete A;
399  delete b;
400  delete bt;
401  delete xexact;
402 
403  GenerateCrsProblem(numNodesX, numNodesY, numProcsX, numProcsY, XLoff.Length(),
404  XLoff.Values(), YLoff.Values(), nrhs, comm, verbose, summary,
405  mapL, L, bL, btL, xexactL, StaticProfile, true);
406 
407 
408  GenerateCrsProblem(numNodesX, numNodesY, numProcsX, numProcsY, XUoff.Length(),
409  XUoff.Values(), YUoff.Values(), nrhs, comm, verbose, summary,
410  mapU, U, bU, btU, xexactU, StaticProfile, true);
411 
412 
413  runLUMatrixTests(L, bL, btL, xexactL, U, bU, btU, xexactU, StaticProfile, verbose, summary);
414 
415  delete L;
416  delete bL;
417  delete btL;
418  delete xexactL;
419  delete mapL;
420 
421  delete U;
422  delete bU;
423  delete btU;
424  delete xexactU;
425  delete mapU;
426 
427  Epetra_MultiVector q(*map, nrhs);
428  Epetra_MultiVector z(q);
429  Epetra_MultiVector r(q);
430 
431  delete map;
432  q.SetFlopCounter(flopcounter);
433  z.SetFlopCounter(q);
434  r.SetFlopCounter(q);
435 
436  resvec.Resize(nrhs);
437 
438 
439  flopcounter.ResetFlops();
440  timer.ResetStartTime();
441 
442  //10 norms
443  for( int i = 0; i < 10; ++i )
444  q.Norm2( resvec.Values() );
445 
446  elapsed_time = timer.ElapsedTime();
447  total_flops = q.Flops();
448  MFLOPs = total_flops/elapsed_time/1000000.0;
449  if (verbose) cout << "\nTotal MFLOPs for 10 Norm2's= " << MFLOPs << endl;
450 
451  if (summary) {
452  if (comm.NumProc()==1) cout << "Norm2" << '\t';
453  cout << MFLOPs << endl;
454  }
455 
456  flopcounter.ResetFlops();
457  timer.ResetStartTime();
458 
459  //10 dot's
460  for( int i = 0; i < 10; ++i )
461  q.Dot(z, resvec.Values());
462 
463  elapsed_time = timer.ElapsedTime();
464  total_flops = q.Flops();
465  MFLOPs = total_flops/elapsed_time/1000000.0;
466  if (verbose) cout << "Total MFLOPs for 10 Dot's = " << MFLOPs << endl;
467 
468  if (summary) {
469  if (comm.NumProc()==1) cout << "DotProd" << '\t';
470  cout << MFLOPs << endl;
471  }
472 
473  flopcounter.ResetFlops();
474  timer.ResetStartTime();
475 
476  //10 dot's
477  for( int i = 0; i < 10; ++i )
478  q.Update(1.0, z, 1.0, r, 0.0);
479 
480  elapsed_time = timer.ElapsedTime();
481  total_flops = q.Flops();
482  MFLOPs = total_flops/elapsed_time/1000000.0;
483  if (verbose) cout << "Total MFLOPs for 10 Updates= " << MFLOPs << endl;
484 
485  if (summary) {
486  if (comm.NumProc()==1) cout << "Update" << '\t';
487  cout << MFLOPs << endl;
488  }
489  }
490  }
491  }
492 #ifdef EPETRA_MPI
493  MPI_Finalize() ;
494 #endif
495 
496 return ierr ;
497 }
498 
499 // Constructs a 2D PDE finite difference matrix using the list of x and y offsets.
500 //
501 // nx (In) - number of grid points in x direction
502 // ny (In) - number of grid points in y direction
503 // The total number of equations will be nx*ny ordered such that the x direction changes
504 // most rapidly:
505 // First equation is at point (0,0)
506 // Second at (1,0)
507 // ...
508 // nx equation at (nx-1,0)
509 // nx+1st equation at (0,1)
510 
511 // numPoints (In) - number of points in finite difference stencil
512 // xoff (In) - stencil offsets in x direction (of length numPoints)
513 // yoff (In) - stencil offsets in y direction (of length numPoints)
514 // A standard 5-point finite difference stencil would be described as:
515 // numPoints = 5
516 // xoff = [-1, 1, 0, 0, 0]
517 // yoff = [ 0, 0, 0, -1, 1]
518 
519 // nrhs - Number of rhs to generate. (First interface produces vectors, so nrhs is not needed
520 
521 // comm (In) - an Epetra_Comm object describing the parallel machine (numProcs and my proc ID)
522 // map (Out) - Epetra_Map describing distribution of matrix and vectors/multivectors
523 // A (Out) - Epetra_CrsMatrix constructed for nx by ny grid using prescribed stencil
524 // Off-diagonal values are random between 0 and 1. If diagonal is part of stencil,
525 // diagonal will be slightly diag dominant.
526 // b (Out) - Generated RHS. Values satisfy b = A*xexact
527 // bt (Out) - Generated RHS. Values satisfy b = A'*xexact
528 // xexact (Out) - Generated exact solution to Ax = b and b' = A'xexact
529 
530 // Note: Caller of this function is responsible for deleting all output objects.
531 
532 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
533  int * xoff, int * yoff,
534  const Epetra_Comm &comm, bool verbose, bool summary,
535  Epetra_Map *& map,
536  Epetra_CrsMatrix *& A,
537  Epetra_Vector *& b,
538  Epetra_Vector *& bt,
539  Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly) {
540 
541  Epetra_MultiVector * b1, * bt1, * xexact1;
542 
543  GenerateCrsProblem(numNodesX, numNodesY, numProcsX, numProcsY, numPoints,
544  xoff, yoff, 1, comm, verbose, summary,
545  map, A, b1, bt1, xexact1, StaticProfile, MakeLocalOnly);
546 
547  b = dynamic_cast<Epetra_Vector *>(b1);
548  bt = dynamic_cast<Epetra_Vector *>(bt1);
549  xexact = dynamic_cast<Epetra_Vector *>(xexact1);
550 
551  return;
552 }
553 
554 void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
555  int * xoff, int * yoff, int nrhs,
556  const Epetra_Comm &comm, bool verbose, bool summary,
557  Epetra_Map *& map,
558  Epetra_CrsMatrix *& A,
559  Epetra_MultiVector *& b,
560  Epetra_MultiVector *& bt,
561  Epetra_MultiVector *&xexact, bool StaticProfile, bool MakeLocalOnly) {
562 
563  Epetra_Time timer(comm);
564  // Determine my global IDs
565  int * myGlobalElements;
566  GenerateMyGlobalElements(numNodesX, numNodesY, numProcsX, numProcsY, comm.MyPID(), myGlobalElements);
567 
568  int numMyEquations = numNodesX*numNodesY;
569 
570  map = new Epetra_Map(-1, numMyEquations, myGlobalElements, 0, comm); // Create map with 2D block partitioning.
571  delete [] myGlobalElements;
572 
573  int numGlobalEquations = map->NumGlobalElements();
574 
575  int profile = 0; if (StaticProfile) profile = numPoints;
576 
577 #ifdef EPETRA_HAVE_STATICPROFILE
578 
579  if (MakeLocalOnly)
580  A = new Epetra_CrsMatrix(Copy, *map, *map, profile, StaticProfile); // Construct matrix with rowmap=colmap
581  else
582  A = new Epetra_CrsMatrix(Copy, *map, profile, StaticProfile); // Construct matrix
583 
584 #else
585 
586  if (MakeLocalOnly)
587  A = new Epetra_CrsMatrix(Copy, *map, *map, profile); // Construct matrix with rowmap=colmap
588  else
589  A = new Epetra_CrsMatrix(Copy, *map, profile); // Construct matrix
590 
591 #endif
592 
593  int * indices = new int[numPoints];
594  double * values = new double[numPoints];
595 
596  double dnumPoints = (double) numPoints;
597  int nx = numNodesX*numProcsX;
598 
599  for (int i=0; i<numMyEquations; i++) {
600 
601  int rowID = map->GID(i);
602  int numIndices = 0;
603 
604  for (int j=0; j<numPoints; j++) {
605  int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
606  if (colID>-1 && colID<numGlobalEquations) {
607  indices[numIndices] = colID;
608  double value = - ((double) rand())/ ((double) RAND_MAX);
609  if (colID==rowID)
610  values[numIndices++] = dnumPoints - value; // Make diagonal dominant
611  else
612  values[numIndices++] = value;
613  }
614  }
615  //cout << "Building row " << rowID << endl;
616  A->InsertGlobalValues(rowID, numIndices, values, indices);
617  }
618 
619  delete [] indices;
620  delete [] values;
621  double insertTime = timer.ElapsedTime();
622  timer.ResetStartTime();
623  A->FillComplete(false);
624  double fillCompleteTime = timer.ElapsedTime();
625 
626  if (verbose)
627  cout << "Time to insert matrix values = " << insertTime << endl
628  << "Time to complete fill = " << fillCompleteTime << endl;
629  if (summary) {
630  if (comm.NumProc()==1) cout << "InsertTime" << '\t';
631  cout << insertTime << endl;
632  if (comm.NumProc()==1) cout << "FillCompleteTime" << '\t';
633  cout << fillCompleteTime << endl;
634  }
635 
636  if (nrhs<=1) {
637  b = new Epetra_Vector(*map);
638  bt = new Epetra_Vector(*map);
639  xexact = new Epetra_Vector(*map);
640  }
641  else {
642  b = new Epetra_MultiVector(*map, nrhs);
643  bt = new Epetra_MultiVector(*map, nrhs);
644  xexact = new Epetra_MultiVector(*map, nrhs);
645  }
646 
647  xexact->Random(); // Fill xexact with random values
648 
649  A->Multiply(false, *xexact, *b);
650  A->Multiply(true, *xexact, *bt);
651 
652  return;
653 }
654 
655 
656 // Constructs a 2D PDE finite difference matrix using the list of x and y offsets.
657 //
658 // nx (In) - number of grid points in x direction
659 // ny (In) - number of grid points in y direction
660 // The total number of equations will be nx*ny ordered such that the x direction changes
661 // most rapidly:
662 // First equation is at point (0,0)
663 // Second at (1,0)
664 // ...
665 // nx equation at (nx-1,0)
666 // nx+1st equation at (0,1)
667 
668 // numPoints (In) - number of points in finite difference stencil
669 // xoff (In) - stencil offsets in x direction (of length numPoints)
670 // yoff (In) - stencil offsets in y direction (of length numPoints)
671 // A standard 5-point finite difference stencil would be described as:
672 // numPoints = 5
673 // xoff = [-1, 1, 0, 0, 0]
674 // yoff = [ 0, 0, 0, -1, 1]
675 
676 // nsizes (In) - Length of element size list used to create variable block map and matrix
677 // sizes (In) - integer list of element sizes of length nsizes
678 // The map associated with this VbrMatrix will be created by cycling through the sizes list.
679 // For example, if nsize = 3 and sizes = [ 2, 4, 3], the block map will have elementsizes
680 // of 2, 4, 3, 2, 4, 3, ...
681 
682 // nrhs - Number of rhs to generate. (First interface produces vectors, so nrhs is not needed
683 
684 // comm (In) - an Epetra_Comm object describing the parallel machine (numProcs and my proc ID)
685 // map (Out) - Epetra_Map describing distribution of matrix and vectors/multivectors
686 // A (Out) - Epetra_VbrMatrix constructed for nx by ny grid using prescribed stencil
687 // Off-diagonal values are random between 0 and 1. If diagonal is part of stencil,
688 // diagonal will be slightly diag dominant.
689 // b (Out) - Generated RHS. Values satisfy b = A*xexact
690 // bt (Out) - Generated RHS. Values satisfy b = A'*xexact
691 // xexact (Out) - Generated exact solution to Ax = b and b' = A'xexact
692 
693 // Note: Caller of this function is responsible for deleting all output objects.
694 
695 void GenerateVbrProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
696  int * xoff, int * yoff,
697  int nsizes, int * sizes,
698  const Epetra_Comm &comm, bool verbose, bool summary,
699  Epetra_BlockMap *& map,
700  Epetra_VbrMatrix *& A,
701  Epetra_Vector *& b,
702  Epetra_Vector *& bt,
703  Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly) {
704 
705  Epetra_MultiVector * b1, * bt1, * xexact1;
706 
707  GenerateVbrProblem(numNodesX, numNodesY, numProcsX, numProcsY, numPoints,
708  xoff, yoff, nsizes, sizes,
709  1, comm, verbose, summary, map, A, b1, bt1, xexact1, StaticProfile, MakeLocalOnly);
710 
711  b = dynamic_cast<Epetra_Vector *>(b1);
712  bt = dynamic_cast<Epetra_Vector *>(bt1);
713  xexact = dynamic_cast<Epetra_Vector *>(xexact1);
714 
715  return;
716 }
717 
718 void GenerateVbrProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints,
719  int * xoff, int * yoff,
720  int nsizes, int * sizes, int nrhs,
721  const Epetra_Comm &comm, bool verbose, bool summary,
722  Epetra_BlockMap *& map,
723  Epetra_VbrMatrix *& A,
724  Epetra_MultiVector *& b,
725  Epetra_MultiVector *& bt,
726  Epetra_MultiVector *&xexact, bool StaticProfile, bool MakeLocalOnly) {
727 
728  int i, j;
729 
730  // Determine my global IDs
731  int * myGlobalElements;
732  GenerateMyGlobalElements(numNodesX, numNodesY, numProcsX, numProcsY, comm.MyPID(), myGlobalElements);
733 
734  int numMyElements = numNodesX*numNodesY;
735 
736  Epetra_Map ptMap(-1, numMyElements, myGlobalElements, 0, comm); // Create map with 2D block partitioning.
737  delete [] myGlobalElements;
738 
739  int numGlobalEquations = ptMap.NumGlobalElements();
740 
741  Epetra_IntVector elementSizes(ptMap); // This vector will have the list of element sizes
742  for (i=0; i<numMyElements; i++)
743  elementSizes[i] = sizes[ptMap.GID(i)%nsizes]; // cycle through sizes array
744 
745  map = new Epetra_BlockMap(-1, numMyElements, ptMap.MyGlobalElements(), elementSizes.Values(),
746  ptMap.IndexBase(), ptMap.Comm());
747 
748  int profile = 0; if (StaticProfile) profile = numPoints;
749 
750  if (MakeLocalOnly)
751  A = new Epetra_VbrMatrix(Copy, *map, *map, profile); // Construct matrix rowmap=colmap
752  else
753  A = new Epetra_VbrMatrix(Copy, *map, profile); // Construct matrix
754 
755  int * indices = new int[numPoints];
756 
757  // This section of code creates a vector of random values that will be used to create
758  // light-weight dense matrices to pass into the VbrMatrix construction process.
759 
760  int maxElementSize = 0;
761  for (i=0; i< nsizes; i++) maxElementSize = EPETRA_MAX(maxElementSize, sizes[i]);
762 
763  Epetra_LocalMap lmap(maxElementSize*maxElementSize, ptMap.IndexBase(), ptMap.Comm());
764  Epetra_Vector randvec(lmap);
765  randvec.Random();
766  randvec.Scale(-1.0); // Make value negative
767  int nx = numNodesX*numProcsX;
768 
769 
770  for (i=0; i<numMyElements; i++) {
771  int rowID = map->GID(i);
772  int numIndices = 0;
773  int rowDim = sizes[rowID%nsizes];
774  for (j=0; j<numPoints; j++) {
775  int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
776  if (colID>-1 && colID<numGlobalEquations)
777  indices[numIndices++] = colID;
778  }
779 
780  A->BeginInsertGlobalValues(rowID, numIndices, indices);
781 
782  for (j=0; j < numIndices; j++) {
783  int colDim = sizes[indices[j]%nsizes];
784  A->SubmitBlockEntry(&(randvec[0]), rowDim, rowDim, colDim);
785  }
786  A->EndSubmitEntries();
787  }
788 
789  delete [] indices;
790 
791  A->FillComplete();
792 
793  // Compute the InvRowSums of the matrix rows
794  Epetra_Vector invRowSums(A->RowMap());
795  Epetra_Vector rowSums(A->RowMap());
796  A->InvRowSums(invRowSums);
797  rowSums.Reciprocal(invRowSums);
798 
799  // Jam the row sum values into the diagonal of the Vbr matrix (to make it diag dominant)
800  int numBlockDiagonalEntries;
801  int * rowColDims;
802  int * diagoffsets = map->FirstPointInElementList();
803  A->BeginExtractBlockDiagonalView(numBlockDiagonalEntries, rowColDims);
804  for (i=0; i< numBlockDiagonalEntries; i++) {
805  double * diagVals;
806  int diagLDA;
807  A->ExtractBlockDiagonalEntryView(diagVals, diagLDA);
808  int rowDim = map->ElementSize(i);
809  for (j=0; j<rowDim; j++) diagVals[j+j*diagLDA] = rowSums[diagoffsets[i]+j];
810  }
811 
812  if (nrhs<=1) {
813  b = new Epetra_Vector(*map);
814  bt = new Epetra_Vector(*map);
815  xexact = new Epetra_Vector(*map);
816  }
817  else {
818  b = new Epetra_MultiVector(*map, nrhs);
819  bt = new Epetra_MultiVector(*map, nrhs);
820  xexact = new Epetra_MultiVector(*map, nrhs);
821  }
822 
823  xexact->Random(); // Fill xexact with random values
824 
825  A->Multiply(false, *xexact, *b);
826  A->Multiply(true, *xexact, *bt);
827 
828  return;
829 }
830 
831 void GenerateMyGlobalElements(int numNodesX, int numNodesY, int numProcsX, int numProcs,
832  int myPID, int * & myGlobalElements) {
833 
834  myGlobalElements = new int[numNodesX*numNodesY];
835  int myProcX = myPID%numProcsX;
836  int myProcY = myPID/numProcsX;
837  int curGID = myProcY*(numProcsX*numNodesX)*numNodesY+myProcX*numNodesX;
838  for (int j=0; j<numNodesY; j++) {
839  for (int i=0; i<numNodesX; i++) {
840  myGlobalElements[j*numNodesX+i] = curGID+i;
841  }
842  curGID+=numNodesX*numProcsX;
843  }
844  //for (int i=0; i<numNodesX*numNodesY; i++) cout << "MYPID " << myPID <<" GID "<< myGlobalElements[i] << endl;
845 
846  return;
847 }
848 
849 void runMatrixTests(Epetra_CrsMatrix * A, Epetra_MultiVector * b, Epetra_MultiVector * bt,
850  Epetra_MultiVector * xexact, bool StaticProfile, bool verbose, bool summary) {
851 
852  Epetra_MultiVector z(*b);
853  Epetra_MultiVector r(*b);
854  Epetra_SerialDenseVector resvec(b->NumVectors());
855 
856  //Timings
857  Epetra_Flops flopcounter;
858  A->SetFlopCounter(flopcounter);
859  Epetra_Time timer(A->Comm());
860  std::string statdyn = "dynamic";
861  if (StaticProfile) statdyn = "static ";
862 
863  for (int j=0; j<4; j++) { // j = 0/2 is notrans, j = 1/3 is trans
864 
865  bool TransA = (j==1 || j==3);
866  std::string contig = "without";
867  if (j>1) contig = "with ";
868 
869 #ifdef EPETRA_SHORT_PERFTEST
870  int kstart = 1;
871 #else
872  int kstart = 0;
873 #endif
874  for (int k=kstart; k<2; k++) { // Loop over old multiply vs. new multiply
875 
876  std::string oldnew = "old";
877  if (k>0) oldnew = "new";
878 
879  if (j==2) A->OptimizeStorage();
880 
881  flopcounter.ResetFlops();
882  timer.ResetStartTime();
883 
884  if (k==0) {
885  //10 matvecs
886 #ifndef EPETRA_SHORT_PERFTEST
887  for( int i = 0; i < 10; ++i )
888  A->Multiply1(TransA, *xexact, z); // Compute z = A*xexact or z = A'*xexact using old Multiply method
889 #endif
890  }
891  else {
892  //10 matvecs
893  for( int i = 0; i < 10; ++i )
894  A->Multiply(TransA, *xexact, z); // Compute z = A*xexact or z = A'*xexact
895  }
896 
897  double elapsed_time = timer.ElapsedTime();
898  double total_flops = A->Flops();
899 
900  // Compute residual
901  if (TransA)
902  r.Update(-1.0, z, 1.0, *bt, 0.0); // r = bt - z
903  else
904  r.Update(-1.0, z, 1.0, *b, 0.0); // r = b - z
905 
906  r.Norm2(resvec.Values());
907 
908  if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
909  double MFLOPs = total_flops/elapsed_time/1000000.0;
910  if (verbose) cout << "Total MFLOPs for 10 " << oldnew << " MatVec's with " << statdyn << " Profile (Trans = " << TransA
911  << ") and " << contig << " optimized storage = " << MFLOPs << " (" << elapsed_time << " s)" <<endl;
912  if (summary) {
913  if (A->Comm().NumProc()==1) {
914  if (TransA) cout << "TransMv" << statdyn<< "Prof" << contig << "OptStor" << '\t';
915  else cout << "NoTransMv" << statdyn << "Prof" << contig << "OptStor" << '\t';
916  }
917  cout << MFLOPs << endl;
918  }
919  }
920  }
921  return;
922 }
923 #ifdef EPETRA_HAVE_JADMATRIX
924 void runJadMatrixTests(Epetra_JadMatrix * A, Epetra_MultiVector * b, Epetra_MultiVector * bt,
925  Epetra_MultiVector * xexact, bool StaticProfile, bool verbose, bool summary) {
926 
927  Epetra_MultiVector z(*b);
928  Epetra_MultiVector r(*b);
929  Epetra_SerialDenseVector resvec(b->NumVectors());
930 
931  //Timings
932  Epetra_Flops flopcounter;
933  A->SetFlopCounter(flopcounter);
934  Epetra_Time timer(A->Comm());
935 
936  for (int j=0; j<2; j++) { // j = 0 is notrans, j = 1 is trans
937 
938  bool TransA = (j==1);
939  A->SetUseTranspose(TransA);
940  flopcounter.ResetFlops();
941  timer.ResetStartTime();
942 
943  //10 matvecs
944  for( int i = 0; i < 10; ++i )
945  A->Apply(*xexact, z); // Compute z = A*xexact or z = A'*xexact
946 
947  double elapsed_time = timer.ElapsedTime();
948  double total_flops = A->Flops();
949 
950  // Compute residual
951  if (TransA)
952  r.Update(-1.0, z, 1.0, *bt, 0.0); // r = bt - z
953  else
954  r.Update(-1.0, z, 1.0, *b, 0.0); // r = b - z
955 
956  r.Norm2(resvec.Values());
957 
958  if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
959  double MFLOPs = total_flops/elapsed_time/1000000.0;
960  if (verbose) cout << "Total MFLOPs for 10 " << " Jagged Diagonal MatVec's with (Trans = " << TransA
961  << ") " << MFLOPs << " (" << elapsed_time << " s)" <<endl;
962  if (summary) {
963  if (A->Comm().NumProc()==1) {
964  if (TransA) cout << "TransMv" << '\t';
965  else cout << "NoTransMv" << '\t';
966  }
967  cout << MFLOPs << endl;
968  }
969  }
970  return;
971 }
972 #endif
973 //=========================================================================================
974 void runLUMatrixTests(Epetra_CrsMatrix * L, Epetra_MultiVector * bL, Epetra_MultiVector * btL, Epetra_MultiVector * xexactL,
975  Epetra_CrsMatrix * U, Epetra_MultiVector * bU, Epetra_MultiVector * btU, Epetra_MultiVector * xexactU,
976  bool StaticProfile, bool verbose, bool summary) {
977 
978  if (L->NoDiagonal()) {
979  bL->Update(1.0, *xexactL, 1.0); // Add contribution of a unit diagonal to bL
980  btL->Update(1.0, *xexactL, 1.0); // Add contribution of a unit diagonal to btL
981  }
982  if (U->NoDiagonal()) {
983  bU->Update(1.0, *xexactU, 1.0); // Add contribution of a unit diagonal to bU
984  btU->Update(1.0, *xexactU, 1.0); // Add contribution of a unit diagonal to btU
985  }
986 
987  Epetra_MultiVector z(*bL);
988  Epetra_MultiVector r(*bL);
989  Epetra_SerialDenseVector resvec(bL->NumVectors());
990 
991  //Timings
992  Epetra_Flops flopcounter;
993  L->SetFlopCounter(flopcounter);
994  U->SetFlopCounter(flopcounter);
995  Epetra_Time timer(L->Comm());
996  std::string statdyn = "dynamic";
997  if (StaticProfile) statdyn = "static ";
998 
999  for (int j=0; j<4; j++) { // j = 0/2 is notrans, j = 1/3 is trans
1000 
1001  bool TransA = (j==1 || j==3);
1002  std::string contig = "without";
1003  if (j>1) contig = "with ";
1004 
1005  if (j==2) {
1006  L->OptimizeStorage();
1007  U->OptimizeStorage();
1008  }
1009 
1010  flopcounter.ResetFlops();
1011  timer.ResetStartTime();
1012 
1013  //10 lower solves
1014  bool Upper = false;
1015  bool UnitDiagonal = L->NoDiagonal(); // If no diagonal, then unit must be used
1016  Epetra_MultiVector * b = TransA ? btL : bL; // solve with the appropriate b vector
1017  for( int i = 0; i < 10; ++i )
1018  L->Solve(Upper, TransA, UnitDiagonal, *b, z); // Solve Lz = bL or L'z = bLt
1019 
1020  double elapsed_time = timer.ElapsedTime();
1021  double total_flops = L->Flops();
1022 
1023  // Compute residual
1024  r.Update(-1.0, z, 1.0, *xexactL, 0.0); // r = bt - z
1025  r.Norm2(resvec.Values());
1026 
1027  if (resvec.NormInf()>0.000001) {
1028  cout << "resvec = " << resvec << endl;
1029  cout << "z = " << z << endl;
1030  cout << "xexactL = " << *xexactL << endl;
1031  cout << "r = " << r << endl;
1032  }
1033 
1034  if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
1035  double MFLOPs = total_flops/elapsed_time/1000000.0;
1036  if (verbose) cout << "Total MFLOPs for 10 " << " Lower solves " << statdyn << " Profile (Trans = " << TransA
1037  << ") and " << contig << " opt storage = " << MFLOPs << " (" << elapsed_time << " s)" <<endl;
1038  if (summary) {
1039  if (L->Comm().NumProc()==1) {
1040  if (TransA) cout << "TransLSv" << statdyn<< "Prof" << contig << "OptStor" << '\t';
1041  else cout << "NoTransLSv" << statdyn << "Prof" << contig << "OptStor" << '\t';
1042  }
1043  cout << MFLOPs << endl;
1044  }
1045  flopcounter.ResetFlops();
1046  timer.ResetStartTime();
1047 
1048  //10 upper solves
1049  Upper = true;
1050  UnitDiagonal = U->NoDiagonal(); // If no diagonal, then unit must be used
1051  b = TransA ? btU : bU; // solve with the appropriate b vector
1052  for( int i = 0; i < 10; ++i )
1053  U->Solve(Upper, TransA, UnitDiagonal, *b, z); // Solve Lz = bL or L'z = bLt
1054 
1055  elapsed_time = timer.ElapsedTime();
1056  total_flops = U->Flops();
1057 
1058  // Compute residual
1059  r.Update(-1.0, z, 1.0, *xexactU, 0.0); // r = bt - z
1060  r.Norm2(resvec.Values());
1061 
1062  if (resvec.NormInf()>0.001) {
1063  cout << "U = " << *U << endl;
1064  //cout << "resvec = " << resvec << endl;
1065  cout << "z = " << z << endl;
1066  cout << "xexactU = " << *xexactU << endl;
1067  //cout << "r = " << r << endl;
1068  cout << "b = " << *b << endl;
1069  }
1070 
1071 
1072  if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
1073  MFLOPs = total_flops/elapsed_time/1000000.0;
1074  if (verbose) cout << "Total MFLOPs for 10 " << " Upper solves " << statdyn << " Profile (Trans = " << TransA
1075  << ") and " << contig << " opt storage = " << MFLOPs <<endl;
1076  if (summary) {
1077  if (L->Comm().NumProc()==1) {
1078  if (TransA) cout << "TransUSv" << statdyn<< "Prof" << contig << "OptStor" << '\t';
1079  else cout << "NoTransUSv" << statdyn << "Prof" << contig << "OptStor" << '\t';
1080  }
1081  cout << MFLOPs << endl;
1082  }
1083  }
1084  return;
1085 }
Epetra_VbrMatrix::ExtractBlockDiagonalEntryView
int ExtractBlockDiagonalEntryView(double *&Values, int &LDA) const
Extract a view of a block diagonal entry from the matrix.
Definition: Epetra_VbrMatrix.cpp:1486
Epetra_MultiVector
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
Definition: Epetra_MultiVector.h:184
Epetra_CompObject::SetFlopCounter
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
Definition: Epetra_CompObject.h:78
Epetra_BlockMap::NumGlobalElements
int NumGlobalElements() const
Number of elements across all processors.
Definition: Epetra_BlockMap.h:546
Epetra_BlockMap::FirstPointInElementList
int * FirstPointInElementList() const
Pointer to internal array containing a mapping between the local elements and the first local point n...
Definition: Epetra_BlockMap.cpp:971
Epetra_Map
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_BasicRowMatrix::SetUseTranspose
virtual int SetUseTranspose(bool use_transpose)
If set true, transpose of this operator will be applied.
Definition: Epetra_BasicRowMatrix.h:447
Epetra_MultiVector::Random
int Random()
Set multi-vector values to random numbers.
Definition: Epetra_MultiVector.cpp:490
Epetra_VbrMatrix::Multiply
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_VbrMatrix multiplied by a Epetra_MultiVector X in Y.
Definition: Epetra_VbrMatrix.cpp:3333
Epetra_IntSerialDenseVector
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
Definition: Epetra_IntSerialDenseVector.h:87
Epetra_BlockMap::ElementSize
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
Definition: Epetra_BlockMap.h:573
Epetra_IntSerialDenseVector::Length
int Length() const
Returns length of vector.
Definition: Epetra_IntSerialDenseVector.h:208
Epetra_SerialDenseVector.h
Copy
Definition: Epetra_DataAccess.h:55
Epetra_CrsMatrix::Multiply
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
Definition: Epetra_CrsMatrix.cpp:3028
Epetra_IntVector
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer...
Definition: Epetra_IntVector.h:124
Epetra_VbrMatrix::BeginExtractBlockDiagonalView
int BeginExtractBlockDiagonalView(int &NumBlockDiagonalEntries, int *&RowColDims) const
Initiates a view of the block diagonal entries.
Definition: Epetra_VbrMatrix.cpp:1476
Epetra_BlockMap::MyGlobalElements
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Definition: Epetra_BlockMap.cpp:848
Epetra_BasicRowMatrix::Comm
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
Definition: Epetra_BasicRowMatrix.h:424
Epetra_BlockMap.h
Epetra_Time::ElapsedTime
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Definition: Epetra_Time.cpp:135
Epetra_IntSerialDenseVector::Size
int Size(int Length_in)
Set length of a Epetra_IntSerialDenseVector object; init values to zero.
Definition: Epetra_IntSerialDenseVector.h:139
Epetra_CrsMatrix::InsertGlobalValues
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
Definition: Epetra_CrsMatrix.cpp:540
Epetra_Object::SetTracebackMode
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Definition: Epetra_Object.cpp:77
GenerateVbrProblem
void GenerateVbrProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints, int *xoff, int *yoff, int nsizes, int *sizes, const Epetra_Comm &comm, bool verbose, bool summary, Epetra_BlockMap *&map, Epetra_VbrMatrix *&A, Epetra_Vector *&b, Epetra_Vector *&bt, Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly)
Definition: test/BasicPerfTest/cxx_main.cpp:695
Epetra_Flops.h
Epetra_VbrMatrix
Epetra_VbrMatrix: A class for the construction and use of real-valued double-precision variable block...
Definition: Epetra_VbrMatrix.h:172
Epetra_SerialDenseVector::NormInf
double NormInf() const
Compute Inf-norm of each vector in multi-vector.
Definition: Epetra_SerialDenseVector.cpp:128
Epetra_Vector
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
Definition: Epetra_Vector.h:142
runJadMatrixTests
void runJadMatrixTests(Epetra_JadMatrix *A, Epetra_MultiVector *b, Epetra_MultiVector *bt, Epetra_MultiVector *xexact, bool StaticProfile, bool verbose, bool summary)
Definition: test/BasicPerfTest/cxx_main.cpp:924
Epetra_SerialComm::MyPID
int MyPID() const
Return my process ID.
Definition: Epetra_SerialComm.h:432
Epetra_VbrMatrix.h
Epetra_MpiComm
Epetra_MpiComm: The Epetra MPI Communication Class.
Definition: Epetra_MpiComm.h:64
Epetra_Version
std::string Epetra_Version()
Definition: Epetra_Version.h:47
Epetra_IntVector::Values
int * Values() const
Returns a pointer to an array containing the values of this vector.
Definition: Epetra_IntVector.h:247
Epetra_SerialComm.h
Epetra_MultiVector::NumVectors
int NumVectors() const
Returns the number of vectors in the multi-vector.
Definition: Epetra_MultiVector.h:925
Epetra_JadMatrix
Epetra_JadMatrix: A class for constructing matrix objects optimized for common kernels.
Definition: Epetra_JadMatrix.h:68
Epetra_Comm::MyPID
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsMatrix::FillComplete
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
Definition: Epetra_CrsMatrix.cpp:1140
Epetra_SerialDenseVector::Resize
int Resize(int Length_in)
Resize a Epetra_SerialDenseVector object.
Definition: Epetra_SerialDenseVector.h:171
Epetra_Time
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
Epetra_MultiVector::Dot
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
Definition: Epetra_MultiVector.cpp:1123
Epetra_BlockMap::IndexBase
int IndexBase() const
Index base for this map.
Definition: Epetra_BlockMap.h:586
Epetra_SerialDenseVector
Epetra_SerialDenseVector: A class for constructing and using dense vectors.
Definition: Epetra_SerialDenseVector.h:95
Epetra_JadMatrix.h
Epetra_Comm
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
GenerateCrsProblem
void GenerateCrsProblem(int numNodesX, int numNodesY, int numProcsX, int numProcsY, int numPoints, int *xoff, int *yoff, const Epetra_Comm &comm, bool verbose, bool summary, Epetra_Map *&map, Epetra_CrsMatrix *&A, Epetra_Vector *&b, Epetra_Vector *&bt, Epetra_Vector *&xexact, bool StaticProfile, bool MakeLocalOnly)
Definition: test/BasicPerfTest/cxx_main.cpp:532
Epetra_BasicRowMatrix::Apply
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_RowMatrix applied to a Epetra_MultiVector X in Y.
Definition: Epetra_BasicRowMatrix.h:459
Epetra_CrsMatrix::Solve
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a local solve using the Epetra_CrsMatrix on a Epetra_Vector x in y...
Definition: Epetra_CrsMatrix.cpp:1629
Epetra_CrsMatrix::OptimizeStorage
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
Definition: Epetra_CrsMatrix.cpp:1279
Epetra_BlockMap::GID
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
Definition: Epetra_BlockMap.cpp:1282
Epetra_CrsMatrix::NoDiagonal
bool NoDiagonal() const
If matrix has no diagonal entries in global index space, this query returns true, otherwise it return...
Definition: Epetra_CrsMatrix.h:1028
Epetra_BlockMap
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
Definition: Epetra_BlockMap.h:194
Epetra_VbrMatrix::RowMap
const Epetra_BlockMap & RowMap() const
Returns the RowMap object as an Epetra_BlockMap (the Epetra_Map base class) needed for implementing E...
Definition: Epetra_VbrMatrix.h:1051
Epetra_VbrMatrix::EndSubmitEntries
int EndSubmitEntries()
Completes processing of all data passed in for the current block row.
Definition: Epetra_VbrMatrix.cpp:705
Epetra_SerialComm::NumProc
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
Definition: Epetra_SerialComm.h:435
Epetra_BlockMap::Comm
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
Definition: Epetra_BlockMap.h:770
Epetra_SerialComm
Epetra_SerialComm: The Epetra Serial Communication Class.
Definition: Epetra_SerialComm.h:61
Epetra_Flops::ResetFlops
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
Epetra_Flops
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
Epetra_Time.h
Epetra_MpiComm.h
Epetra_IntVector.h
Epetra_IntSerialDenseVector.h
Epetra_MultiVector::Update
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
Definition: Epetra_MultiVector.cpp:1276
Epetra_Comm::NumProc
virtual int NumProc() const =0
Returns total number of processes.
Epetra_VbrMatrix::InvRowSums
int InvRowSums(Epetra_Vector &x) const
Computes the sum of absolute values of the rows of the Epetra_VbrMatrix, results returned in x...
Definition: Epetra_VbrMatrix.cpp:2479
Epetra_LocalMap.h
Epetra_CrsMatrix
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
Definition: Epetra_CrsMatrix.h:173
Epetra_MultiVector::Norm2
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
Definition: Epetra_MultiVector.cpp:1547
Epetra_CompObject::Flops
double Flops() const
Returns the number of floating point operations with this multi-vector.
Definition: Epetra_CompObject.h:93
Epetra_Map.h
Epetra_MultiVector.h
Epetra_VbrMatrix::FillComplete
int FillComplete()
Signal that data entry is complete, perform transformations to local index space. ...
Definition: Epetra_VbrMatrix.cpp:875
Epetra_Vector.h
Epetra_VbrMatrix::BeginInsertGlobalValues
int BeginInsertGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices)
Initiate insertion of a list of elements in a given global row of the matrix, values are inserted via...
Definition: Epetra_VbrMatrix.cpp:522
Epetra_VbrMatrix::SubmitBlockEntry
int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols)
Submit a block entry to the indicated block row and column specified in the Begin routine...
Definition: Epetra_VbrMatrix.cpp:682
main
int main(int argc, char *argv[])
Definition: example/C_wrappers/c_main.c:49
runMatrixTests
void runMatrixTests(Epetra_CrsMatrix *A, Epetra_MultiVector *b, Epetra_MultiVector *bt, Epetra_MultiVector *xexact, bool StaticProfile, bool verbose, bool summary)
Definition: test/BasicPerfTest/cxx_main.cpp:849
Epetra_Version.h
GenerateMyGlobalElements
void GenerateMyGlobalElements(int numNodesX, int numNodesY, int numProcsX, int numProcs, int myPID, int *&myGlobalElements)
Definition: test/BasicPerfTest/cxx_main.cpp:831
Epetra_SerialDenseVector::Values
double * Values() const
Returns pointer to the values in vector.
Definition: Epetra_SerialDenseVector.h:274
Epetra_IntSerialDenseVector::Values
int * Values()
Returns pointer to the values in vector.
Definition: Epetra_IntSerialDenseVector.h:211
Epetra_CrsMatrix.h
runLUMatrixTests
void runLUMatrixTests(Epetra_CrsMatrix *L, Epetra_MultiVector *bL, Epetra_MultiVector *btL, Epetra_MultiVector *xexactL, Epetra_CrsMatrix *U, Epetra_MultiVector *bU, Epetra_MultiVector *btU, Epetra_MultiVector *xexactU, bool StaticProfile, bool verbose, bool summary)
Definition: test/BasicPerfTest/cxx_main.cpp:974
Epetra_Time::ResetStartTime
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
Definition: Epetra_Time.cpp:129
Epetra_CrsMatrix::Multiply1
int Multiply1(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Definition: Epetra_CrsMatrix.cpp:4050
Epetra_CrsMatrix::Comm
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
Definition: Epetra_CrsMatrix.h:1251
EPETRA_MAX
#define EPETRA_MAX(x, y)
Definition: Epetra_ConfigDefs.h:62
Epetra_LocalMap
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
Definition: Epetra_LocalMap.h:89
Generated by   doxygen 1.8.5