Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/CrsRectMatrix_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_LocalMap.h"
44 #include "Epetra_Map.h"
45 #include "Epetra_Time.h"
46 #include "Epetra_Vector.h"
47 #include "Epetra_CrsMatrix.h"
48 #include "Epetra_Flops.h"
49 #include "Epetra_Export.h"
50 #include "Epetra_Import.h"
53 #ifdef EPETRA_MPI
54 #include "Epetra_MpiComm.h"
55 #include "mpi.h"
56 #else
57 #include "Epetra_SerialComm.h"
58 #endif
59 #include "../epetra_test_err.h"
60 #include "Epetra_Version.h"
61 
62 int main(int argc, char *argv[])
63 {
64  int ierr = 0, i, j, forierr = 0;
65  int ntrials = 1;
66 
67 #ifdef EPETRA_MPI
68 
69  // Initialize MPI
70 
71  MPI_Init(&argc,&argv);
72  Epetra_MpiComm Comm( MPI_COMM_WORLD );
73 
74 #else
75 
76  Epetra_SerialComm Comm;
77 
78 #endif
79 
80  bool verbose = false;
81 
82  // Check if we should print results to standard out
83  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
84 
85 
86 
87  // char tmp; if (Comm.MyPID()==0) { cout << "Press any key to continue..."<< endl; cin >> tmp;} Comm.Barrier();
88 
89  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
90  int MyPID = Comm.MyPID();
91  int NumProc = Comm.NumProc();
92 
93  if (verbose && Comm.MyPID()==0)
94  cout << Epetra_Version() << endl << endl;
95 
96  if (verbose)
97  cout << "Processor " << MyPID << " of " << NumProc << " is alive." << endl;
98 
99  bool verbose1 = verbose;
100 
101  // Redefine verbose to only print on PE 0
102  if (verbose && Comm.MyPID()!=0) verbose = false;
103 
104  int NumMyEquations = 10000;
105 
106  long long NumGlobalEquations = (long long)(NumMyEquations) * NumProc;
107  long long NumGlobalVariables = 2 * NumGlobalEquations+1;
108 
109  // Construct a Map that puts approximately the same Number of equations on each processor
110 
111  Epetra_Map RowMap(NumGlobalEquations, 0LL, Comm);
112  Epetra_Map XMap(NumGlobalVariables, 0LL, Comm);
113  Epetra_Map& YMap = RowMap;
114 
115  // Get update list and number of local equations from newly created Map
116  long long * MyGlobalElements = new long long[RowMap.NumMyElements()];
117  RowMap.MyGlobalElements(MyGlobalElements);
118 
119  // Get update list and number of local equations from newly created XMap
120  long long * XGlobalElements = new long long[XMap.NumMyElements()];
121  XMap.MyGlobalElements(XGlobalElements);
122 
123  // Get update list and number of local variables from newly created YMap
124  long long * YGlobalElements = new long long[YMap.NumMyElements()];
125  YMap.MyGlobalElements(YGlobalElements);
126 
127  // We need vectors to compute:
128  // X = A^T*Y
129  // AATY = A*A^T*Y = A*X
130  // and
131  // BY = B*Y
132 
133  Epetra_Vector Y(YMap);
134  Epetra_Vector X(XMap);
135  Epetra_Vector AATY(YMap);
136  Epetra_Vector BY(YMap);
137 
138 
139  // Fill Y Vector
140  Y.Random();
141  //Y.PutScalar(1.0);
142 
143  // To create A^T explicitly we need an assembly map that is two elements longer than
144  // the XMap, because each processor will be making contributions to two rows beyond what
145  // it will own.
146  long long ATAssemblyNumMyElements = 2*MyGlobalElements[NumMyEquations-1] + 2 - 2*MyGlobalElements[0] + 1;
147  long long * ATAssemblyGlobalElements = new long long[ATAssemblyNumMyElements];
148 
149  for (i=0; i<ATAssemblyNumMyElements; i++) ATAssemblyGlobalElements[i] = 2*MyGlobalElements[0] + i;
150  Epetra_Map ATAssemblyMap((long long)-1, ATAssemblyNumMyElements, ATAssemblyGlobalElements, 0LL, Comm);
151 
152  // Create a Epetra_Matrix with the values of A
153  // A is a simple 1D weighted average operator that mimics a restriction operator
154  // that might be found in a multigrid code.
155  // Also create A^T explicitly
156 
157  Epetra_CrsMatrix A(Copy, RowMap, 3);
158  Epetra_CrsMatrix ATAssembly(Copy, ATAssemblyMap, 2);
159  Epetra_CrsMatrix AT(Copy, XMap, 2);
160 
161  //cout << "ATAssemblyMap = "<< endl<< ATAssemblyMap << endl
162  // << "XMap = " << endl << XMap << endl
163  // << "RowMap = " << endl << RowMap << endl;
164  // Add rows one-at-a-time
165  // Need some vectors to help
166  // Off diagonal Values will always be -1
167 
168 
169  double *Values = new double[3];
170  long long *Indices = new long long[3];
171  int NumEntries;
172  /*
173  Values[0] = 0.25;
174  Values[1] = 0.5;
175  Values[2] = 0.25;
176  */
177  Values[0] = 0.5;
178  Values[1] = 0.25;
179  Values[2] = 0.25;
180  forierr = 0;
181  for (i=0; i<NumMyEquations; i++)
182  {
183  /*
184  Indices[0] = 2*MyGlobalElements[i];
185  Indices[1] = 2*MyGlobalElements[i]+1;
186  Indices[2] = 2*MyGlobalElements[i]+2;
187  */
188  Indices[0] = 2*MyGlobalElements[i]+1;
189  Indices[1] = 2*MyGlobalElements[i]+2;
190  Indices[2] = 2*MyGlobalElements[i];
191  NumEntries = 3;
192  forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
193  for (j=0; j<3; j++)
194  forierr += !(ATAssembly.InsertGlobalValues(Indices[j],1, &(Values[j]), &(MyGlobalElements[i]))>=0);
195  }
196  EPETRA_TEST_ERR(forierr,ierr);
197 
198 
199  EPETRA_TEST_ERR(!(ATAssembly.FillComplete()==0),ierr);
200  // Gather AT values from ATAssembly matrix
201  Epetra_Export Exporter(ATAssemblyMap, XMap);
202  EPETRA_TEST_ERR(!(AT.Export(ATAssembly, Exporter, Add)==0),ierr);
203 
204  // Finish up
205  EPETRA_TEST_ERR(!(A.FillComplete(XMap, YMap)==0),ierr);
206  EPETRA_TEST_ERR(!(AT.FillComplete(YMap, XMap)==0),ierr);
207 
208 
209  if (verbose1 && NumGlobalEquations<20) {
210  if (verbose) cout << "\n\n Matrix A\n" << endl;
211  cout << A << endl;
212  if (verbose) cout << " \n\n Matrix A Transpose\n" << endl;
213  cout << AT << endl;
214  }
215 
216 
217  // Create a Epetra_Matrix containing B = A*A^T.
218  // This matrix will be a square tridiagonal matrix. We will use it to compare the results
219  // of A*(A^T*X) using two methods: (1) with two calls to Multiply using A^T and then A and
220  // (2) using B directly.
221 
222  Epetra_CrsMatrix B(Copy, RowMap, 3);
223 
224  Values[0] = 1.0/16.0;
225  Values[1] = 3.0/8.0;
226  Values[2] = 1.0/16.0;
227  int Valstart;
228  forierr = 0;
229  for (i=0; i<NumMyEquations; i++)
230  {
231  if (MyGlobalElements[i] == 0) {
232  Indices[0] = MyGlobalElements[i];
233  Indices[1] = MyGlobalElements[i]+1;
234  NumEntries = 2;
235  Valstart = 1;
236  }
237  else {
238  Indices[0] = MyGlobalElements[i]-1;
239  Indices[1] = MyGlobalElements[i];
240  Indices[2] = MyGlobalElements[i]+1;
241  NumEntries = 3;
242  Valstart = 0;
243  }
244  if (MyGlobalElements[i] == NumGlobalEquations-1) NumEntries--;
245  forierr += !(B.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values+Valstart, Indices)==0);
246  }
247  EPETRA_TEST_ERR(forierr,ierr);
248 
249  // Finish up
250  EPETRA_TEST_ERR(!(B.FillComplete()==0),ierr);
251  if (verbose && NumGlobalEquations<20) cout << "\n\nMatrix B \n" << endl;
252  if (verbose1 && NumGlobalEquations<20) cout << B << endl;
253 
254 
255  Epetra_Flops counter;
256  A.SetFlopCounter(counter);
257  B.SetFlopCounter(A);
258  Epetra_Time timer(Comm);
259  for (i=0; i<ntrials; i++) {
260  EPETRA_TEST_ERR(!(B.Multiply(false, Y, BY)==0),ierr); // Compute BY = B*Y
261  }
262  double elapsed_time = timer.ElapsedTime();
263  double total_flops = B.Flops();
264  counter.ResetFlops();
265 
266  double MFLOPs = total_flops/elapsed_time/1000000.0;
267 
268  if (verbose) cout << "\n\nTotal MFLOPs for B*Y = " << MFLOPs << endl<< endl;
269  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = B*Y \n";
270  if (verbose1 && NumGlobalEquations<20) cout << BY << endl;
271 
272 
274 
275  timer.ResetStartTime();
276  for (i=0; i<ntrials; i++) {
277  EPETRA_TEST_ERR(!(A.Multiply(true, Y, X)==0),ierr); // Compute X = A^T*Y
278  }
279  elapsed_time = timer.ElapsedTime();
280  total_flops = A.Flops();
281  counter.ResetFlops();
282  MFLOPs = total_flops/elapsed_time/1000000.0;
283 
284  if (verbose) cout << "\n\nTotal MFLOPs for A^T*Y using A and trans=true = " << MFLOPs << endl<< endl;
285  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = AT*Y \n";
286  if (verbose1 && NumGlobalEquations<20) cout << X << endl;
287 
289 
290  timer.ResetStartTime();
291  EPETRA_TEST_ERR(!(A.Multiply(false, X, AATY)==0),ierr); // Compute AATY = A*X
292  elapsed_time = timer.ElapsedTime();
293  total_flops = A.Flops();
294  MFLOPs = total_flops/elapsed_time/1000000.0;
295  counter.ResetFlops();
296  Epetra_Vector resid(YMap);
297  resid.Update(1.0, BY, -1.0, AATY, 0.0);
298  double residual;
299  resid.Norm2(&residual);
300 
301  if (verbose) cout << "\n\nTotal MFLOPs for A*X using A and trans=false = " << MFLOPs << endl<< endl;
302  if (verbose) cout << "Residual = " << residual << endl<< endl;
303  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = A*ATY \n";
304  if (verbose1 && NumGlobalEquations<20) cout << AATY << endl;
305 
307 
308  AT.SetFlopCounter(counter);
309  timer.ResetStartTime();
310  for (i=0; i<ntrials; i++) {
311  EPETRA_TEST_ERR(!(AT.Multiply(false, Y, X)==0),ierr); // Compute X = A^T*Y
312  }
313  elapsed_time = timer.ElapsedTime();
314  total_flops = AT.Flops();
315  counter.ResetFlops();
316  MFLOPs = total_flops/elapsed_time/1000000.0;
317 
318  if (verbose) cout << "\n\nTotal MFLOPs for A^T*Y using AT and trans=false = " << MFLOPs << endl<< endl;
319  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = AT*Y \n";
320  if (verbose1 && NumGlobalEquations<20) cout << X << endl;
321 
323 
324  timer.ResetStartTime();
325  for (i=0; i<ntrials; i++) {
326  EPETRA_TEST_ERR(!(AT.Multiply(true, X, AATY)==0),ierr); // Compute AATY = A*X
327  }
328  elapsed_time = timer.ElapsedTime();
329  total_flops = AT.Flops();
330  MFLOPs = total_flops/elapsed_time/1000000.0;
331  counter.ResetFlops();
332  resid.Update(1.0, BY, -1.0, AATY, 0.0);
333  resid.Norm2(&residual);
334 
335  if (verbose) cout << "\n\nTotal MFLOPs for A*X using AT and trans=true = " << MFLOPs << endl<< endl;
336  if (verbose) cout << "Residual = " << residual << endl<< endl;
337  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Z = A*ATY \n";
338  if (verbose1 && NumGlobalEquations<20) cout <<AATY << endl;
339 
340 
342  // Now test use of Epetra_LocalMap vectors: First test case of local replicated range vector
343 
344  {
345  Epetra_CrsMatrix AL(Copy, RowMap, 3);
346  for (i=0; i<NumMyEquations; i++)
347  {
348  forierr += !(A.ExtractGlobalRowCopy(MyGlobalElements[i], 3, NumEntries, Values, Indices)==0);
349  forierr += !(AL.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
350  }
351  EPETRA_TEST_ERR(forierr,ierr);
352 
353  Epetra_LocalMap YLMap(NumGlobalEquations, 0LL, Comm);
354  EPETRA_TEST_ERR(!(AL.FillComplete(XMap, YLMap)==0),ierr);
355  AL.SetFlopCounter(A);
356  Epetra_Vector YL(YLMap);
357  Epetra_Vector ALX(YLMap);
358 
359  timer.ResetStartTime();
360  for (i=0; i<ntrials; i++) {
361  EPETRA_TEST_ERR(!(A.Multiply(false, X, Y)==0),ierr); // Compute Y= A*X
362  }
363  elapsed_time = timer.ElapsedTime();
364  total_flops = A.Flops();
365  counter.ResetFlops();
366  MFLOPs = total_flops/elapsed_time/1000000.0;
367 
368 
369 
370  if (verbose) cout << "\n\nTotal MFLOPs for Y=A*X using global distributed Y = " << MFLOPs << endl<< endl;
371  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Y = A*X using distributed Y \n";
372  if (verbose1 && NumGlobalEquations<20) cout << Y << endl;
373  if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist Y range map\n";
374  if (verbose1 && NumGlobalEquations<20) cout << A << endl;
375 
376  timer.ResetStartTime();
377  for (i=0; i<ntrials; i++) {
378  EPETRA_TEST_ERR(!(AL.Multiply(false, X, ALX)==0),ierr); // Compute YL= A*X
379  }
380  elapsed_time = timer.ElapsedTime();
381  total_flops = AL.Flops();
382  counter.ResetFlops();
383  MFLOPs = total_flops/elapsed_time/1000000.0;
384 
385  if (verbose) cout << "\n\nTotal MFLOPs for Y=A*X using Local replicated Y = " << MFLOPs << endl<< endl;
386  if (verbose && NumGlobalEquations<20) cout << "\n\nVector YL = AL*X using local replicated Y \n";
387  if (verbose1 && NumGlobalEquations<20) cout << ALX << endl;
388  if (verbose && NumGlobalEquations<20) cout << "\n\nA using local Y range map\n";
389  if (verbose1 && NumGlobalEquations<20) cout << AL << endl;
390 
391  // Now gather Y values from the distributed Y and compare them to the local replicated Y values
392  Epetra_Import g2limporter(YLMap, YMap); // Import from distributed Y map to local Y map
393  EPETRA_TEST_ERR(!(YL.Import(Y, g2limporter, Insert)==0),ierr);
394  if (verbose && NumGlobalEquations<20) cout << "\n\nVector YL = imported from distributed Y \n";
395  if (verbose1 && NumGlobalEquations<20) cout << YL << endl;
396  EPETRA_TEST_ERR(!(YL.Update(-1.0, ALX, 1.0)==0),ierr);
397  EPETRA_TEST_ERR(!(YL.Norm2(&residual)==0),ierr);
398  if (verbose) cout << "Residual = " << residual << endl<< endl;
399 
400 
401  //
402  // Multiply by transpose
403  //
404 
405  timer.ResetStartTime();
406  for (i=0; i<ntrials; i++) {
407  EPETRA_TEST_ERR(!(A.Multiply(true, Y, X)==0),ierr); // Compute X = A^TY
408  }
409  elapsed_time = timer.ElapsedTime();
410  total_flops = A.Flops();
411  counter.ResetFlops();
412  MFLOPs = total_flops/elapsed_time/1000000.0;
413 
414 
415 
416  if (verbose) cout << "\n\nTotal MFLOPs for X=A^TY using global distributed Y = " << MFLOPs << endl<< endl;
417  if (verbose && NumGlobalEquations<20) cout << "\n\nVector X using distributed Y \n";
418  if (verbose1 && NumGlobalEquations<20) cout << Y << endl;
419  if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist Y range map\n";
420  if (verbose1 && NumGlobalEquations<20) cout << A << endl;
421 
422  Epetra_Vector X1(XMap);
423 
424  timer.ResetStartTime();
425  for (i=0; i<ntrials; i++) {
426  EPETRA_TEST_ERR(!(AL.Multiply(true, ALX, X1)==0),ierr); // Compute X1 = AL^T*Y
427  }
428  elapsed_time = timer.ElapsedTime();
429  total_flops = AL.Flops();
430  counter.ResetFlops();
431  MFLOPs = total_flops/elapsed_time/1000000.0;
432 
433  if (verbose) cout << "\n\nTotal MFLOPs for X1=A^T*Y using Local replicated Y = " << MFLOPs << endl<< endl;
434  if (verbose && NumGlobalEquations<20) cout << "\n\nVector X1 using local replicated Y \n";
435  if (verbose1 && NumGlobalEquations<20) cout << X1 << endl;
436 
437  EPETRA_TEST_ERR(!(X1.Update(-1.0, X, 1.0)==0),ierr);
438  EPETRA_TEST_ERR(!(X1.Norm2(&residual)==0),ierr);
439  if (verbose) cout << "Residual = " << residual << endl<< endl;
440  }
441 
443  // Finally test use of Epetra_LocalMap vectors using local replicated domain vector
444 
445  {
446  Epetra_CrsMatrix AL(Copy, RowMap, 3);
447  for (i=0; i<NumMyEquations; i++)
448  {
449  forierr += !(A.ExtractGlobalRowCopy(MyGlobalElements[i], 3, NumEntries, Values, Indices)==0);
450  forierr += !(AL.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
451  }
452  EPETRA_TEST_ERR(forierr,ierr);
453 
454  Epetra_LocalMap XLMap(NumGlobalVariables, 0LL, Comm);
455  EPETRA_TEST_ERR(!(AL.FillComplete(XLMap, YMap)==0),ierr);
456  AL.SetFlopCounter(A);
457  Epetra_Vector XL(XLMap);
458  Epetra_Vector ALX(XLMap);
459 
460  timer.ResetStartTime();
461  for (i=0; i<ntrials; i++) {
462  EPETRA_TEST_ERR(!(A.Multiply(false, X, Y)==0),ierr); // Compute Y= A*X
463  }
464  elapsed_time = timer.ElapsedTime();
465  total_flops = A.Flops();
466  counter.ResetFlops();
467  MFLOPs = total_flops/elapsed_time/1000000.0;
468 
469 
470 
471  if (verbose) cout << "\n\nTotal MFLOPs for Y=A*X using global distributed X = " << MFLOPs << endl<< endl;
472  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Y = A*X using distributed X \n";
473  if (verbose1 && NumGlobalEquations<20) cout << Y << endl;
474  //if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist X range map\n";
475  //if (verbose1 && NumGlobalEquations<20) cout << A << endl;
476 
477  // Now gather X values from the distributed X
478  Epetra_Import g2limporter(XLMap, XMap); // Import from distributed X map to local X map
479  EPETRA_TEST_ERR(!(XL.Import(X, g2limporter, Insert)==0),ierr);
480  if (verbose && NumGlobalEquations<20) cout << "\n\nVector XL = imported from distributed X \n";
481  if (verbose1 && NumGlobalEquations<20) cout << XL << endl;
482  Epetra_Vector Y1(Y);
483  timer.ResetStartTime();
484  for (i=0; i<ntrials; i++) {
485  EPETRA_TEST_ERR(!(AL.Multiply(false, XL, Y1)==0),ierr); // Compute Y1= AL*XL
486  }
487  elapsed_time = timer.ElapsedTime();
488  total_flops = AL.Flops();
489  counter.ResetFlops();
490  MFLOPs = total_flops/elapsed_time/1000000.0;
491 
492  if (verbose) cout << "\n\nTotal MFLOPs for Y1=A*XL using Local replicated X = " << MFLOPs << endl<< endl;
493  if (verbose && NumGlobalEquations<20) cout << "\n\nVector Y1 = AL*XL using local replicated X \n";
494  if (verbose1 && NumGlobalEquations<20) cout << Y1 << endl;
495  //if (verbose && NumGlobalEquations<20) cout << "\n\nA using local X domain map\n";
496  //if (verbose1 && NumGlobalEquations<20) cout << AL << endl;
497 
498  EPETRA_TEST_ERR(!(Y1.Update(-1.0, Y, 1.0)==0),ierr);
499  EPETRA_TEST_ERR(!(Y1.Norm2(&residual)==0),ierr);
500  if (verbose) cout << "Residual = " << residual << endl<< endl;
501 
502 
503  //
504  // Multiply by transpose
505  //
506  timer.ResetStartTime();
507  for (i=0; i<ntrials; i++) {
508  EPETRA_TEST_ERR(!(A.Multiply(true, Y, X)==0),ierr); // Compute X = A^TY
509  }
510  elapsed_time = timer.ElapsedTime();
511  total_flops = A.Flops();
512  counter.ResetFlops();
513  MFLOPs = total_flops/elapsed_time/1000000.0;
514 
515 
516  if (verbose) cout << "\n\nTotal MFLOPs for X=A^TY using global distributed X = " << MFLOPs << endl<< endl;
517  if (verbose && NumGlobalEquations<20) cout << "\n\nVector X using distributed X \n";
518  if (verbose1 && NumGlobalEquations<20) cout << X << endl;
519  //if (verbose && NumGlobalEquations<20) cout << "\n\nA using dist X domain map\n";
520  //if (verbose1 && NumGlobalEquations<20) cout << A << endl;
521 
522  timer.ResetStartTime();
523  for (i=0; i<ntrials; i++) {
524  EPETRA_TEST_ERR(!(AL.Multiply(true, Y, XL)==0),ierr); // Compute XL = AL^T*Y
525  }
526  elapsed_time = timer.ElapsedTime();
527  total_flops = AL.Flops();
528  counter.ResetFlops();
529  MFLOPs = total_flops/elapsed_time/1000000.0;
530 
531  if (verbose) cout << "\n\nTotal MFLOPs for XL=A^T*Y1 using Local replicated X = " << MFLOPs << endl<< endl;
532  if (verbose && NumGlobalEquations<20) cout << "\n\nVector XL using local replicated X \n";
533  if (verbose1 && NumGlobalEquations<20) cout << XL << endl;
534 
535  Epetra_Vector XL1(XLMap);
536  EPETRA_TEST_ERR(!(XL1.Import(X, g2limporter, Insert)==0),ierr);
537  EPETRA_TEST_ERR(!(XL1.Update(-1.0, XL, 1.0)==0),ierr);
538  EPETRA_TEST_ERR(!(XL1.Norm2(&residual)==0),ierr);
539  if (verbose) cout << "Residual = " << residual << endl<< endl;
540  }
541  // Release all objects
542  delete [] Values;
543  delete [] Indices;
544  delete [] MyGlobalElements;
545  delete [] XGlobalElements;
546  delete [] YGlobalElements;
547  delete [] ATAssemblyGlobalElements;
548 
549 #ifdef EPETRA_MPI
550  MPI_Finalize() ;
551 #endif
552 
553  /* end main
554  */
555  return ierr ;
556 }
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:127
int Random()
Set multi-vector values to random numbers.
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.
#define EPETRA_TEST_ERR(a, b)
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
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.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements...
Definition: Epetra_Export.h:70
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
std::string Epetra_Version()
Epetra_Import: This class builds an import object for efficient importing of off-processor elements...
Definition: Epetra_Import.h:71
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:83
int NumMyElements() const
Number of elements on the calling processor.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified global row in user-provided arrays.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:85
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:66
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
double Flops() const
Returns the number of floating point operations with this multi-vector.
int main(int argc, char *argv[])
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.