Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/CrsMatrix/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"
45 #include "Epetra_CrsMatrix.h"
46 #include "Epetra_Vector.h"
47 #include "Epetra_Flops.h"
48 #ifdef EPETRA_MPI
49 #include "Epetra_MpiComm.h"
50 #include "mpi.h"
51 #else
52 #include "Epetra_SerialComm.h"
53 #endif
54 #include "../epetra_test_err.h"
55 #include "Epetra_Version.h"
56 
57 // prototypes
58 
59 int check(Epetra_CrsMatrix& A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1,
60  int NumGlobalNonzeros1, int * MyGlobalElements, bool verbose);
61 
62 int power_method(bool TransA, Epetra_CrsMatrix& A,
63  Epetra_Vector& q,
64  Epetra_Vector& z,
65  Epetra_Vector& resid,
66  double * lambda, int niters, double tolerance,
67  bool verbose);
68 
70 
71 int main(int argc, char *argv[])
72 {
73  int ierr = 0, forierr = 0;
74  bool debug = false;
75 
76 #ifdef EPETRA_MPI
77 
78  // Initialize MPI
79 
80  MPI_Init(&argc,&argv);
81  int rank; // My process ID
82 
83  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
84  Epetra_MpiComm Comm( MPI_COMM_WORLD );
85 
86 #else
87 
88  int rank = 0;
89  Epetra_SerialComm Comm;
90 
91 #endif
92 
93  bool verbose = false;
94 
95  // Check if we should print results to standard out
96  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
97 
98  int verbose_int = verbose ? 1 : 0;
99  Comm.Broadcast(&verbose_int, 1, 0);
100  verbose = verbose_int==1 ? true : false;
101 
102 
103  // char tmp;
104  // if (rank==0) cout << "Press any key to continue..."<< std::endl;
105  // if (rank==0) cin >> tmp;
106  // Comm.Barrier();
107 
108  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
109  int MyPID = Comm.MyPID();
110  int NumProc = Comm.NumProc();
111 
112  if(verbose && MyPID==0)
113  cout << Epetra_Version() << std::endl << std::endl;
114 
115  if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
116  << " is alive."<<endl;
117 
118  bool verbose1 = verbose;
119 
120  // Redefine verbose to only print on PE 0
121  if(verbose && rank!=0)
122  verbose = false;
123 
124  int NumMyEquations = 10000;
125  int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);
126  if(MyPID < 3)
127  NumMyEquations++;
128 
129  // Construct a Map that puts approximately the same Number of equations on each processor
130 
131  Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
132 
133  // Get update list and number of local equations from newly created Map
134  int* MyGlobalElements = new int[Map.NumMyElements()];
135  Map.MyGlobalElements(MyGlobalElements);
136 
137  // Create an integer vector NumNz that is used to build the Petra Matrix.
138  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
139 
140  int* NumNz = new int[NumMyEquations];
141 
142  // We are building a tridiagonal matrix where each row has (-1 2 -1)
143  // So we need 2 off-diagonal terms (except for the first and last equation)
144 
145  for (int i = 0; i < NumMyEquations; i++)
146  if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
147  NumNz[i] = 1;
148  else
149  NumNz[i] = 2;
150 
151  // Create a Epetra_Matrix
152 
153  Epetra_CrsMatrix A(Copy, Map, NumNz);
156 
157  // Add rows one-at-a-time
158  // Need some vectors to help
159  // Off diagonal Values will always be -1
160 
161 
162  double* Values = new double[2];
163  Values[0] = -1.0;
164  Values[1] = -1.0;
165  int* Indices = new int[2];
166  double two = 2.0;
167  int NumEntries;
168 
169  forierr = 0;
170  for (int i = 0; i < NumMyEquations; i++) {
171  if(MyGlobalElements[i] == 0) {
172  Indices[0] = 1;
173  NumEntries = 1;
174  }
175  else if (MyGlobalElements[i] == NumGlobalEquations-1) {
176  Indices[0] = NumGlobalEquations-2;
177  NumEntries = 1;
178  }
179  else {
180  Indices[0] = MyGlobalElements[i]-1;
181  Indices[1] = MyGlobalElements[i]+1;
182  NumEntries = 2;
183  }
184  forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices)==0);
185  forierr += !(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i)>0); // Put in the diagonal entry
186  }
187  EPETRA_TEST_ERR(forierr,ierr);
188 
189  int * indexOffsetTmp;
190  int * indicesTmp;
191  double * valuesTmp;
192  // Finish up
193  EPETRA_TEST_ERR(!(A.IndicesAreGlobal()),ierr);
194  EPETRA_TEST_ERR(!(A.ExtractCrsDataPointers(indexOffsetTmp, indicesTmp, valuesTmp)==-1),ierr); // Should fail
195  EPETRA_TEST_ERR(!(A.FillComplete(false)==0),ierr);
196  EPETRA_TEST_ERR(!(A.ExtractCrsDataPointers(indexOffsetTmp, indicesTmp, valuesTmp)==-1),ierr); // Should fail
197  EPETRA_TEST_ERR(!(A.IndicesAreLocal()),ierr);
199  A.OptimizeStorage();
200  EPETRA_TEST_ERR(!(A.StorageOptimized()),ierr);
201  EPETRA_TEST_ERR(!(A.ExtractCrsDataPointers(indexOffsetTmp, indicesTmp, valuesTmp)==0),ierr); // Should succeed
202  const Epetra_CrsGraph & GofA = A.Graph();
203  EPETRA_TEST_ERR((indicesTmp!=GofA[0] || valuesTmp!=A[0]),ierr); // Extra check to see if operator[] is consistent
206 
207  int NumMyNonzeros = 3 * NumMyEquations;
208  if(A.LRID(0) >= 0)
209  NumMyNonzeros--; // If I own first global row, then there is one less nonzero
210  if(A.LRID(NumGlobalEquations-1) >= 0)
211  NumMyNonzeros--; // If I own last global row, then there is one less nonzero
212  EPETRA_TEST_ERR(check(A, NumMyEquations, NumGlobalEquations, NumMyNonzeros, 3*NumGlobalEquations-2,
213  MyGlobalElements, verbose),ierr);
214  forierr = 0;
215  for (int i = 0; i < NumMyEquations; i++)
216  forierr += !(A.NumGlobalEntries(MyGlobalElements[i])==NumNz[i]+1);
217  EPETRA_TEST_ERR(forierr,ierr);
218  forierr = 0;
219  for (int i = 0; i < NumMyEquations; i++)
220  forierr += !(A.NumMyEntries(i)==NumNz[i]+1);
221  EPETRA_TEST_ERR(forierr,ierr);
222 
223  if (verbose) cout << "\n\nNumEntries function check OK" << std::endl<< std::endl;
224 
226 
227  // Create vectors for Power method
228 
229  Epetra_Vector q(Map);
230  Epetra_Vector z(Map);
231  Epetra_Vector resid(Map);
232 
233  // variable needed for iteration
234  double lambda = 0.0;
235  // int niters = 10000;
236  int niters = 200;
237  double tolerance = 1.0e-1;
238 
240 
241  // Iterate
242 
243  Epetra_Flops flopcounter;
244  A.SetFlopCounter(flopcounter);
245  q.SetFlopCounter(A);
246  z.SetFlopCounter(A);
247  resid.SetFlopCounter(A);
248 
249 
250  Epetra_Time timer(Comm);
251  EPETRA_TEST_ERR(power_method(false, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
252  double elapsed_time = timer.ElapsedTime();
253  double total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
254  double MFLOPs = total_flops/elapsed_time/1000000.0;
255 
256  if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
257 
259 
260  // Solve transpose problem
261 
262  if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
263  << std::endl;
264  // Iterate
265  lambda = 0.0;
266  flopcounter.ResetFlops();
267  timer.ResetStartTime();
268  EPETRA_TEST_ERR(power_method(true, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
269  elapsed_time = timer.ElapsedTime();
270  total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
271  MFLOPs = total_flops/elapsed_time/1000000.0;
272 
273  if (verbose) cout << "\n\nTotal MFLOPs for transpose solve = " << MFLOPs << std::endl<< endl;
274 
276 
277  // Increase diagonal dominance
278 
279  if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
280  << endl;
281 
282 
283  if (A.MyGlobalRow(0)) {
284  int numvals = A.NumGlobalEntries(0);
285  double * Rowvals = new double [numvals];
286  int * Rowinds = new int [numvals];
287  A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds); // Get A[0,0]
288 
289  for (int i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
290 
291  A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
292  delete [] Rowvals;
293  delete [] Rowinds;
294  }
295  // Iterate (again)
296  lambda = 0.0;
297  flopcounter.ResetFlops();
298  timer.ResetStartTime();
299  EPETRA_TEST_ERR(power_method(false, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
300  elapsed_time = timer.ElapsedTime();
301  total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
302  MFLOPs = total_flops/elapsed_time/1000000.0;
303 
304  if (verbose) cout << "\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
305 
307 
308  // Solve transpose problem
309 
310  if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
311  << endl;
312 
313  // Iterate (again)
314  lambda = 0.0;
315  flopcounter.ResetFlops();
316  timer.ResetStartTime();
317  EPETRA_TEST_ERR(power_method(true, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
318  elapsed_time = timer.ElapsedTime();
319  total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
320  MFLOPs = total_flops/elapsed_time/1000000.0;
321 
322 
323  if (verbose) cout << "\n\nTotal MFLOPs for tranpose of second solve = " << MFLOPs << endl<< endl;
324 
325  if (verbose) cout << "\n\n*****Testing constant entry constructor" << endl<< endl;
326 
327  Epetra_CrsMatrix AA(Copy, Map, 5);
328 
329  if (debug) Comm.Barrier();
330 
331  double dble_one = 1.0;
332  for (int i=0; i< NumMyEquations; i++) AA.InsertGlobalValues(MyGlobalElements[i], 1, &dble_one, MyGlobalElements+i);
333 
334  // Note: All processors will call the following Insert routines, but only the processor
335  // that owns it will actually do anything
336 
337  int One = 1;
338  if (AA.MyGlobalRow(0)) {
339  EPETRA_TEST_ERR(!(AA.InsertGlobalValues(0, 0, &dble_one, &One)==0),ierr);
340  }
341  else EPETRA_TEST_ERR(!(AA.InsertGlobalValues(0, 1, &dble_one, &One)==-1),ierr);
342  EPETRA_TEST_ERR(!(AA.FillComplete(false)==0),ierr);
344  EPETRA_TEST_ERR(!(AA.UpperTriangular()),ierr);
345  EPETRA_TEST_ERR(!(AA.LowerTriangular()),ierr);
346 
347  if (debug) Comm.Barrier();
348  EPETRA_TEST_ERR(check(AA, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
349  MyGlobalElements, verbose),ierr);
350 
351  if (debug) Comm.Barrier();
352 
353  forierr = 0;
354  for (int i=0; i<NumMyEquations; i++) forierr += !(AA.NumGlobalEntries(MyGlobalElements[i])==1);
355  EPETRA_TEST_ERR(forierr,ierr);
356 
357  if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
358 
359  if (debug) Comm.Barrier();
360 
361  if (verbose) cout << "\n\n*****Testing copy constructor" << endl<< endl;
362 
363  Epetra_CrsMatrix B(AA);
364  EPETRA_TEST_ERR(check(B, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
365  MyGlobalElements, verbose),ierr);
366 
367  forierr = 0;
368  for (int i=0; i<NumMyEquations; i++) forierr += !(B.NumGlobalEntries(MyGlobalElements[i])==1);
369  EPETRA_TEST_ERR(forierr,ierr);
370 
371  if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
372 
373  if (debug) Comm.Barrier();
374 
375  if (verbose) cout << "\n\n*****Testing local view constructor" << endl<< endl;
376 
377  Epetra_CrsMatrix BV(View, AA.RowMap(), AA.ColMap(), 0);
378 
379  forierr = 0;
380  int* Inds;
381  double* Vals;
382  for (int i = 0; i < NumMyEquations; i++) {
383  forierr += !(AA.ExtractMyRowView(i, NumEntries, Vals, Inds)==0);
384  forierr += !(BV.InsertMyValues(i, NumEntries, Vals, Inds)==0);
385  }
386  BV.FillComplete(false);
387  EPETRA_TEST_ERR(check(BV, NumMyEquations, NumGlobalEquations, NumMyEquations, NumGlobalEquations,
388  MyGlobalElements, verbose),ierr);
389 
390  forierr = 0;
391  for (int i=0; i<NumMyEquations; i++) forierr += !(BV.NumGlobalEntries(MyGlobalElements[i])==1);
392  EPETRA_TEST_ERR(forierr,ierr);
393 
394  if (verbose) cout << "\n\nNumEntries function check OK" << endl<< endl;
395 
396  if (debug) Comm.Barrier();
397  if (verbose) cout << "\n\n*****Testing post construction modifications" << endl<< endl;
398 
399  EPETRA_TEST_ERR(!(B.InsertGlobalValues(0, 1, &dble_one, &One)==-2),ierr);
400 
401 
402  // Release all objects
403  delete [] NumNz;
404  delete [] Values;
405  delete [] Indices;
406  delete [] MyGlobalElements;
407 
408 
409  if (verbose1) {
410  // Test ostream << operator (if verbose1)
411  // Construct a Map that puts 2 equations on each PE
412 
413  int NumMyElements1 = 2;
414  int NumMyEquations1 = NumMyElements1;
415  int NumGlobalEquations1 = NumMyEquations1*NumProc;
416 
417  Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
418 
419  // Get update list and number of local equations from newly created Map
420  int * MyGlobalElements1 = new int[Map1.NumMyElements()];
421  Map1.MyGlobalElements(MyGlobalElements1);
422 
423  // Create an integer vector NumNz that is used to build the Petra Matrix.
424  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
425 
426  int * NumNz1 = new int[NumMyEquations1];
427 
428  // We are building a tridiagonal matrix where each row has (-1 2 -1)
429  // So we need 2 off-diagonal terms (except for the first and last equation)
430 
431  for (int i=0; i<NumMyEquations1; i++)
432  if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalEquations1-1)
433  NumNz1[i] = 1;
434  else
435  NumNz1[i] = 2;
436 
437  // Create a Epetra_Matrix
438 
439  Epetra_CrsMatrix A1(Copy, Map1, NumNz1);
440 
441  // Add rows one-at-a-time
442  // Need some vectors to help
443  // Off diagonal Values will always be -1
444 
445 
446  double *Values1 = new double[2];
447  Values1[0] = -1.0; Values1[1] = -1.0;
448  int *Indices1 = new int[2];
449  double two1 = 2.0;
450  int NumEntries1;
451 
452  forierr = 0;
453  for (int i=0; i<NumMyEquations1; i++)
454  {
455  if (MyGlobalElements1[i]==0)
456  {
457  Indices1[0] = 1;
458  NumEntries1 = 1;
459  }
460  else if (MyGlobalElements1[i] == NumGlobalEquations1-1)
461  {
462  Indices1[0] = NumGlobalEquations1-2;
463  NumEntries1 = 1;
464  }
465  else
466  {
467  Indices1[0] = MyGlobalElements1[i]-1;
468  Indices1[1] = MyGlobalElements1[i]+1;
469  NumEntries1 = 2;
470  }
471  forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1)==0);
472  forierr += !(A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i)>0); // Put in the diagonal entry
473  }
474  EPETRA_TEST_ERR(forierr,ierr);
475  delete [] Indices1;
476  delete [] Values1;
477 
478  // Finish up
479  EPETRA_TEST_ERR(!(A1.FillComplete(false)==0),ierr);
480 
481  // Test diagonal extraction function
482 
483  Epetra_Vector checkDiag(Map1);
484  EPETRA_TEST_ERR(!(A1.ExtractDiagonalCopy(checkDiag)==0),ierr);
485 
486  forierr = 0;
487  for (int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==two1);
488  EPETRA_TEST_ERR(forierr,ierr);
489 
490  // Test diagonal replacement method
491 
492  forierr = 0;
493  for (int i=0; i<NumMyEquations1; i++) checkDiag[i]=two1*two1;
494  EPETRA_TEST_ERR(forierr,ierr);
495 
496  EPETRA_TEST_ERR(!(A1.ReplaceDiagonalValues(checkDiag)==0),ierr);
497 
498  Epetra_Vector checkDiag1(Map1);
499  EPETRA_TEST_ERR(!(A1.ExtractDiagonalCopy(checkDiag1)==0),ierr);
500 
501  forierr = 0;
502  for (int i=0; i<NumMyEquations1; i++) forierr += !(checkDiag[i]==checkDiag1[i]);
503  EPETRA_TEST_ERR(forierr,ierr);
504 
505  if (verbose) cout << "\n\nDiagonal extraction and replacement OK.\n\n" << endl;
506 
507  double orignorm = A1.NormOne();
508  EPETRA_TEST_ERR(!(A1.Scale(4.0)==0),ierr);
509  EPETRA_TEST_ERR(!(A1.NormOne()!=orignorm),ierr);
510 
511  if (verbose) cout << "\n\nMatrix scale OK.\n\n" << endl;
512 
513  if (verbose) cout << "\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
514  cout << A1 << endl;
515 
516 
517  // Release all objects
518  delete [] NumNz1;
519  delete [] MyGlobalElements1;
520 
521  }
522 
523  if (verbose) cout << "\n\n*****Testing LeftScale and RightScale" << endl << endl;
524 
525  int NumMyElements2 = 7;
526  int NumMyRows2 = 1;//This value should not be changed without editing the
527  // code below.
528  Epetra_Map RowMap(-1,NumMyRows2,0,Comm);
529  Epetra_Map ColMap(NumMyElements2,NumMyElements2,0,Comm);
530  // The DomainMap needs to be different from the ColMap for the test to
531  // be meaningful.
532  Epetra_Map DomainMap(NumMyElements2,0,Comm);
533  int NumMyRangeElements2 = 0;
534  // We need to distribute the elements differently for the range map also.
535  if (MyPID % 2 == 0)
536  NumMyRangeElements2 = NumMyRows2*2; //put elements on even number procs
537  if (NumProc % 2 == 1 && MyPID == NumProc-1)
538  NumMyRangeElements2 = NumMyRows2; //If number of procs is odd, put
539  // the last NumMyElements2 elements on the last proc
540  Epetra_Map RangeMap(-1,NumMyRangeElements2,0,Comm);
541  Epetra_CrsMatrix A2(Copy,RowMap,ColMap,NumMyElements2);
542  double * Values2 = new double[NumMyElements2];
543  int * Indices2 = new int[NumMyElements2];
544 
545  for (int i=0; i<NumMyElements2; i++) {
546  Values2[i] = i+MyPID;
547  Indices2[i]=i;
548  }
549 
550  A2.InsertMyValues(0,NumMyElements2,Values2,Indices2);
551  A2.FillComplete(DomainMap,RangeMap,false);
552  Epetra_CrsMatrix A2copy(A2);
553 
554  double * RowLeftScaleValues = new double[NumMyRows2];
555  double * ColRightScaleValues = new double[NumMyElements2];
556  int RowLoopLength = RowMap.MaxMyGID()-RowMap.MinMyGID()+1;
557  for (int i=0; i<RowLoopLength; i++)
558  RowLeftScaleValues[i] = (i + RowMap.MinMyGID() ) % 2 + 1;
559  // For the column map, all procs own all elements
560  for (int i=0; i<NumMyElements2;i++)
561  ColRightScaleValues[i] = i % 2 + 1;
562 
563  int RangeLoopLength = RangeMap.MaxMyGID()-RangeMap.MinMyGID()+1;
564  double * RangeLeftScaleValues = new double[RangeLoopLength];
565  int DomainLoopLength = DomainMap.MaxMyGID()-DomainMap.MinMyGID()+1;
566  double * DomainRightScaleValues = new double[DomainLoopLength];
567  for (int i=0; i<RangeLoopLength; i++)
568  RangeLeftScaleValues[i] = 1.0/((i + RangeMap.MinMyGID() ) % 2 + 1);
569  for (int i=0; i<DomainLoopLength;i++)
570  DomainRightScaleValues[i] = 1.0/((i + DomainMap.MinMyGID() ) % 2 + 1);
571 
572  Epetra_Vector xRow(View,RowMap,RowLeftScaleValues);
573  Epetra_Vector xCol(View,ColMap,ColRightScaleValues);
574  Epetra_Vector xRange(View,RangeMap,RangeLeftScaleValues);
575  Epetra_Vector xDomain(View,DomainMap,DomainRightScaleValues);
576 
577  double A2infNorm = A2.NormInf();
578  double A2oneNorm = A2.NormOne();
579 
580  if (verbose1) cout << A2;
581  EPETRA_TEST_ERR(A2.LeftScale(xRow),ierr);
582  double A2infNorm1 = A2.NormInf();
583  double A2oneNorm1 = A2.NormOne();
584  bool ScalingBroke = false;
585  if (A2infNorm1>2*A2infNorm||A2infNorm1<A2infNorm) {
586  EPETRA_TEST_ERR(-31,ierr);
587  ScalingBroke = true;
588  }
589  if (A2oneNorm1>2*A2oneNorm||A2oneNorm1<A2oneNorm) {
590 
591  EPETRA_TEST_ERR(-32,ierr);
592  ScalingBroke = true;
593  }
594  if (verbose1) cout << A2;
595  EPETRA_TEST_ERR(A2.RightScale(xCol),ierr);
596  double A2infNorm2 = A2.NormInf();
597  double A2oneNorm2 = A2.NormOne();
598  if (A2infNorm2>=2*A2infNorm1||A2infNorm2<=A2infNorm1) {
599  EPETRA_TEST_ERR(-33,ierr);
600  ScalingBroke = true;
601  }
602  if (A2oneNorm2>2*A2oneNorm1||A2oneNorm2<=A2oneNorm1) {
603  EPETRA_TEST_ERR(-34,ierr);
604  ScalingBroke = true;
605  }
606  if (verbose1) cout << A2;
607  EPETRA_TEST_ERR(A2.RightScale(xDomain),ierr);
608  double A2infNorm3 = A2.NormInf();
609  double A2oneNorm3 = A2.NormOne();
610  // The last two scaling ops cancel each other out
611  if (A2infNorm3!=A2infNorm1) {
612  EPETRA_TEST_ERR(-35,ierr)
613  ScalingBroke = true;
614  }
615  if (A2oneNorm3!=A2oneNorm1) {
616  EPETRA_TEST_ERR(-36,ierr)
617  ScalingBroke = true;
618  }
619  if (verbose1) cout << A2;
620  EPETRA_TEST_ERR(A2.LeftScale(xRange),ierr);
621  double A2infNorm4 = A2.NormInf();
622  double A2oneNorm4 = A2.NormOne();
623  // The 4 scaling ops all cancel out
624  if (A2infNorm4!=A2infNorm) {
625  EPETRA_TEST_ERR(-37,ierr)
626  ScalingBroke = true;
627  }
628  if (A2oneNorm4!=A2oneNorm) {
629  EPETRA_TEST_ERR(-38,ierr)
630  ScalingBroke = true;
631  }
632 
633  //
634  // Now try changing the values underneath and make sure that
635  // telling one process about the change causes NormInf() and
636  // NormOne() to recompute the norm on all processes.
637  //
638 
639  double *values;
640  int num_my_rows = A2.NumMyRows() ;
641  int num_entries;
642 
643  for ( int i=0 ; i< num_my_rows; i++ ) {
644  EPETRA_TEST_ERR( A2.ExtractMyRowView( i, num_entries, values ), ierr );
645  for ( int j = 0 ; j <num_entries; j++ ) {
646  values[j] *= 2.0;
647  }
648  }
649 
650 
651  if ( MyPID == 0 )
652  A2.SumIntoGlobalValues( 0, 0, 0, 0 ) ;
653 
654  double A2infNorm5 = A2.NormInf();
655  double A2oneNorm5 = A2.NormOne();
656 
657  if (A2infNorm5!=2.0 * A2infNorm4) {
658  EPETRA_TEST_ERR(-39,ierr)
659  ScalingBroke = true;
660  }
661  if (A2oneNorm5!= 2.0 * A2oneNorm4) {
662  EPETRA_TEST_ERR(-40,ierr)
663  ScalingBroke = true;
664  }
665 
666  //
667  // Restore the values underneath
668  //
669  for ( int i=0 ; i< num_my_rows; i++ ) {
670  EPETRA_TEST_ERR( A2.ExtractMyRowView( i, num_entries, values ), ierr );
671  for ( int j = 0 ; j <num_entries; j++ ) {
672  values[j] /= 2.0;
673  }
674  }
675 
676  if (verbose1) cout << A2;
677 
678  if (ScalingBroke) {
679  if (verbose) cout << endl << "LeftScale and RightScale tests FAILED" << endl << endl;
680  }
681  else {
682  if (verbose) cout << endl << "LeftScale and RightScale tests PASSED" << endl << endl;
683  }
684 
685  Comm.Barrier();
686 
687  if (verbose) cout << "\n\n*****Testing InvRowMaxs and InvColMaxs" << endl << endl;
688 
689  if (verbose1) cout << A2 << endl;
690  EPETRA_TEST_ERR(A2.InvRowMaxs(xRow),ierr);
691  EPETRA_TEST_ERR(A2.InvRowMaxs(xRange),ierr);
692  if (verbose1) cout << xRow << endl << xRange << endl;
693 
694  if (verbose) cout << "\n\n*****Testing InvRowSums and InvColSums" << endl << endl;
695  bool InvSumsBroke = false;
696 // Works!
697  EPETRA_TEST_ERR(A2.InvRowSums(xRow),ierr);
698  if (verbose1) cout << xRow;
699  EPETRA_TEST_ERR(A2.LeftScale(xRow),ierr);
700  float A2infNormFloat = A2.NormInf();
701  if (verbose1) cout << A2 << endl;
702  if (fabs(1.0-A2infNormFloat) > 1.e-5) {
703  EPETRA_TEST_ERR(-41,ierr);
704  InvSumsBroke = true;
705  }
706 
707  // Works
708  int expectedcode = 1;
709  if (Comm.NumProc()>1) expectedcode = 0;
710  EPETRA_TEST_ERR(!(A2.InvColSums(xDomain)==expectedcode),ierr); // This matrix has a single row, the first column has a zero, so a warning is issued.
711  if (verbose1) cout << xDomain << endl;
712  EPETRA_TEST_ERR(A2.RightScale(xDomain),ierr);
713  float A2oneNormFloat2 = A2.NormOne();
714  if (verbose1) cout << A2;
715  if (fabs(1.0-A2oneNormFloat2)>1.e-5) {
716  EPETRA_TEST_ERR(-42,ierr)
717  InvSumsBroke = true;
718  }
719 
720 // Works!
721  EPETRA_TEST_ERR(A2.InvRowSums(xRange),ierr);
722 
723  if (verbose1) cout << xRange;
724  EPETRA_TEST_ERR(A2.LeftScale(xRange),ierr);
725  float A2infNormFloat2 = A2.NormInf(); // We use a float so that rounding error
726  // will not prevent the sum from being 1.0.
727  if (verbose1) cout << A2;
728  if (fabs(1.0-A2infNormFloat2)>1.e-5) {
729  cout << "InfNorm should be = 1, but InfNorm = " << A2infNormFloat2 << endl;
730  EPETRA_TEST_ERR(-43,ierr);
731  InvSumsBroke = true;
732  }
733 
734  // Doesn't work - may not need this test because column ownership is not unique
735  /* EPETRA_TEST_ERR(A2.InvColSums(xCol),ierr);
736 cout << xCol;
737  EPETRA_TEST_ERR(A2.RightScale(xCol),ierr);
738  float A2oneNormFloat = A2.NormOne();
739 cout << A2;
740  if (fabs(1.0-A2oneNormFloat)>1.e-5) {
741  EPETRA_TEST_ERR(-44,ierr);
742  InvSumsBroke = true;
743  }
744  */
745  delete [] ColRightScaleValues;
746  delete [] DomainRightScaleValues;
747  if (verbose) cout << "Begin partial sum testing." << endl;
748  // Test with a matrix that has partial sums for a subset of the rows
749  // on multiple processors. (Except for the serial case, of course.)
750  int NumMyRows3 = 2; // Changing this requires further changes below
751  int * myGlobalElements = new int[NumMyRows3];
752  for (int i=0; i<NumMyRows3; i++) myGlobalElements[i] = MyPID+i;
753  Epetra_Map RowMap3(NumProc*2, NumMyRows3, myGlobalElements, 0, Comm);
754  int NumMyElements3 = 5;
755  Epetra_CrsMatrix A3(Copy, RowMap3, NumMyElements3);
756  double * Values3 = new double[NumMyElements3];
757  int * Indices3 = new int[NumMyElements3];
758  for (int i=0; i < NumMyElements3; i++) {
759  Values3[i] = (int) (MyPID + (i+1));
760  Indices3[i]=i;
761  }
762  for (int i=0; i<NumMyRows3; i++) {
763  A3.InsertGlobalValues(myGlobalElements[i],NumMyElements3,Values3,Indices3);
764  }
765  Epetra_Map RangeMap3(NumProc+1, 0, Comm);
766  Epetra_Map DomainMap3(NumMyElements3, 0, Comm);
767  EPETRA_TEST_ERR(A3.FillComplete(DomainMap3, RangeMap3,false),ierr);
768  if (verbose1) cout << A3;
769  Epetra_Vector xRange3(RangeMap3,false);
770  Epetra_Vector xDomain3(DomainMap3,false);
771 
772  EPETRA_TEST_ERR(A3.InvRowSums(xRange3),ierr);
773 
774  if (verbose1) cout << xRange3;
775  EPETRA_TEST_ERR(A3.LeftScale(xRange3),ierr);
776  float A3infNormFloat = A3.NormInf();
777  if (verbose1) cout << A3;
778  if (1.0!=A3infNormFloat) {
779  cout << "InfNorm should be = 1, but InfNorm = " << A3infNormFloat <<endl;
780  EPETRA_TEST_ERR(-61,ierr);
781  InvSumsBroke = true;
782  }
783  // we want to take the transpose of our matrix and fill in different values.
784  int NumMyColumns3 = NumMyRows3;
785  Epetra_Map ColMap3cm(RowMap3);
786  Epetra_Map RowMap3cm(A3.ColMap());
787 
788  Epetra_CrsMatrix A3cm(Copy,RowMap3cm,ColMap3cm,NumProc+1);
789  double *Values3cm = new double[NumMyColumns3];
790  int * Indices3cm = new int[NumMyColumns3];
791  for (int i=0; i<NumMyColumns3; i++) {
792  Values3cm[i] = MyPID + i + 1;
793  Indices3cm[i]= i + MyPID;
794  }
795  for (int ii=0; ii<NumMyElements3; ii++) {
796  A3cm.InsertGlobalValues(ii, NumMyColumns3, Values3cm, Indices3cm);
797  }
798 
799  // The DomainMap and the RangeMap from the last test will work fine for
800  // the RangeMap and DomainMap, respectively, but I will make copies to
801  // avaoid confusion when passing what looks like a DomainMap where we
802  // need a RangeMap and vice vera.
803  Epetra_Map RangeMap3cm(DomainMap3);
804  Epetra_Map DomainMap3cm(RangeMap3);
805  EPETRA_TEST_ERR(A3cm.FillComplete(DomainMap3cm,RangeMap3cm),ierr);
806  if (verbose1) cout << A3cm << endl;
807 
808  // Again, we can copy objects from the last example.
809  //Epetra_Vector xRange3cm(xDomain3); //Don't use at this time
810  Epetra_Vector xDomain3cm(DomainMap3cm,false);
811 
812  EPETRA_TEST_ERR(A3cm.InvColSums(xDomain3cm),ierr);
813 
814  if (verbose1) cout << xDomain3cm << endl;
815 
816  EPETRA_TEST_ERR(A3cm.RightScale(xDomain3cm),ierr);
817  float A3cmOneNormFloat = A3cm.NormOne();
818  if (verbose1) cout << A3cm << endl;
819  if (1.0!=A3cmOneNormFloat) {
820  cout << "OneNorm should be = 1, but OneNorm = " << A3cmOneNormFloat << endl;
821  EPETRA_TEST_ERR(-62,ierr);
822  InvSumsBroke = true;
823  }
824 
825  if (verbose) cout << "End partial sum testing" << endl;
826  if (verbose) cout << "Begin replicated testing" << endl;
827 
828  // We will now view the shared row as a repliated row, rather than one
829  // that has partial sums of its entries on mulitple processors.
830  // We will reuse much of the data used for the partial sum tesitng.
831  Epetra_Vector xRow3(RowMap3,false);
832  Epetra_CrsMatrix A4(Copy, RowMap3, NumMyElements3);
833  for (int ii=0; ii < NumMyElements3; ii++) {
834  Values3[ii] = (int)((ii*.6)+1.0);
835  }
836  for (int ii=0; ii<NumMyRows3; ii++) {
837  A4.InsertGlobalValues(myGlobalElements[ii],NumMyElements3,Values3,Indices3);
838  }
839  EPETRA_TEST_ERR(A4.FillComplete(DomainMap3, RangeMap3,false),ierr);
840  if (verbose1) cout << A4 << endl;
841  // The next two lines should be expanded into a verifiable test.
842  EPETRA_TEST_ERR(A4.InvRowMaxs(xRow3),ierr);
843  EPETRA_TEST_ERR(A4.InvRowMaxs(xRange3),ierr);
844  if (verbose1) cout << xRow3 << xRange3;
845 
846  EPETRA_TEST_ERR(A4.InvRowSums(xRow3),ierr);
847  if (verbose1) cout << xRow3;
848  EPETRA_TEST_ERR(A4.LeftScale(xRow3),ierr);
849  float A4infNormFloat = A4.NormInf();
850  if (verbose1) cout << A4;
851  if (2.0!=A4infNormFloat && NumProc != 1) {
852  if (verbose1) cout << "InfNorm should be = 2 (because one column is replicated on two processors and NormOne() does not handle replication), but InfNorm = " << A4infNormFloat <<endl;
853  EPETRA_TEST_ERR(-63,ierr);
854  InvSumsBroke = true;
855  }
856  else if (1.0!=A4infNormFloat && NumProc == 1) {
857  if (verbose1) cout << "InfNorm should be = 1, but InfNorm = " << A4infNormFloat <<endl;
858  EPETRA_TEST_ERR(-63,ierr);
859  InvSumsBroke = true;
860  }
861 
862  Epetra_Vector xCol3cm(ColMap3cm,false);
863  Epetra_CrsMatrix A4cm(Copy, RowMap3cm, ColMap3cm, NumProc+1);
864  //Use values from A3cm
865  for (int ii=0; ii<NumMyElements3; ii++) {
866  A4cm.InsertGlobalValues(ii,NumMyColumns3,Values3cm,Indices3cm);
867  }
868  EPETRA_TEST_ERR(A4cm.FillComplete(DomainMap3cm, RangeMap3cm,false),ierr);
869  if (verbose1) cout << A4cm << endl;
870  // The next two lines should be expanded into a verifiable test.
871  EPETRA_TEST_ERR(A4cm.InvColMaxs(xCol3cm),ierr);
872  EPETRA_TEST_ERR(A4cm.InvColMaxs(xDomain3cm),ierr);
873  if (verbose1) cout << xCol3cm << xDomain3cm;
874 
875  EPETRA_TEST_ERR(A4cm.InvColSums(xCol3cm),ierr);
876 
877  if (verbose1) cout << xCol3cm << endl;
878  EPETRA_TEST_ERR(A4cm.RightScale(xCol3cm),ierr);
879  float A4cmOneNormFloat = A4cm.NormOne();
880  if (verbose1) cout << A4cm << endl;
881  if (2.0!=A4cmOneNormFloat && NumProc != 1) {
882  if (verbose1) cout << "OneNorm should be = 2 (because one column is replicated on two processors and NormOne() does not handle replication), but OneNorm = " << A4cmOneNormFloat << endl;
883  EPETRA_TEST_ERR(-64,ierr);
884  InvSumsBroke = true;
885  }
886  else if (1.0!=A4cmOneNormFloat && NumProc == 1) {
887  if (verbose1) cout << "OneNorm should be = 1, but OneNorm = " << A4infNormFloat <<endl;
888  EPETRA_TEST_ERR(-64,ierr);
889  InvSumsBroke = true;
890  }
891 
892  if (verbose) cout << "End replicated testing" << endl;
893 
894  if (InvSumsBroke) {
895  if (verbose) cout << endl << "InvRowSums tests FAILED" << endl << endl;
896  }
897  else
898  if (verbose) cout << endl << "InvRowSums tests PASSED" << endl << endl;
899 
900  A3cm.PutScalar(2.0);
901  int nnz_A3cm = A3cm.Graph().NumGlobalNonzeros();
902  double check_frobnorm = sqrt(nnz_A3cm*4.0);
903  double frobnorm = A3cm.NormFrobenius();
904 
905  bool frobnorm_test_failed = false;
906  if (fabs(check_frobnorm-frobnorm) > 5.e-5) {
907  frobnorm_test_failed = true;
908  }
909 
910  if (frobnorm_test_failed) {
911  if (verbose) std::cout << "Frobenius-norm test FAILED."<<std::endl;
912  EPETRA_TEST_ERR(-65, ierr);
913  }
914 
915  // Subcommunicator test - only processor 1 has unknowns
916  {
917  int rv=0;
918  int NumMyRows = (MyPID==0) ? NumGlobalEquations : 0;
919  Epetra_Map Map1(-1,NumMyRows, 0, Comm);
920 
921  Epetra_CrsMatrix *A1 = new Epetra_CrsMatrix(Copy, Map1,0);
922  double value = 1.0;
923  for(int i=0; i<NumMyRows; i++) {
924  int GID = Map1.GID(i);
925  EPETRA_TEST_ERR(A1->InsertGlobalValues(GID, 1,&value,&GID),ierr);
926  }
927  EPETRA_TEST_ERR(A1->FillComplete(),ierr);
928 
929  Epetra_BlockMap *Map2 = Map1.RemoveEmptyProcesses();
930  rv=A1->RemoveEmptyProcessesInPlace(Map2);
931  if(rv!=0) {
932  if (verbose) std::cout << "Subcommunicator test FAILED."<<std::endl;
933  EPETRA_TEST_ERR(-66, ierr);
934  }
935 
936  delete Map2;
937  }
938 
939  delete [] Values2;
940  delete [] Indices2;
941  delete [] myGlobalElements;
942  delete [] Values3;
943  delete [] Indices3;
944  delete [] Values3cm;
945  delete [] Indices3cm;
946  delete [] RangeLeftScaleValues;
947  delete [] RowLeftScaleValues;
948 #ifdef EPETRA_MPI
949  MPI_Finalize() ;
950 #endif
951 
952 /* end main
953 */
954 return ierr ;
955 }
956 
958  Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose)
959 {
960 
961  // Fill z with random Numbers
962  z.Random();
963 
964  // variable needed for iteration
965  double normz, residual;
966 
967  int ierr = 1;
968 
969  for(int iter = 0; iter < niters; iter++) {
970  z.Norm2(&normz); // Compute 2-norm of z
971  q.Scale(1.0/normz, z);
972  A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
973  q.Dot(z, lambda); // Approximate maximum eigenvaluE
974  if(iter%100==0 || iter+1==niters) {
975  resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
976  resid.Norm2(&residual);
977  if(verbose) cout << "Iter = " << iter << " Lambda = " << *lambda
978  << " Residual of A*q - lambda*q = " << residual << endl;
979  }
980  if(residual < tolerance) {
981  ierr = 0;
982  break;
983  }
984  }
985  return(ierr);
986 }
987 
988 int check(Epetra_CrsMatrix& A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1,
989  int NumGlobalNonzeros1, int* MyGlobalElements, bool verbose)
990 {
991  (void)MyGlobalElements;
992  int ierr = 0, forierr = 0;
993  int NumGlobalIndices;
994  int NumMyIndices;
995  int* MyViewIndices = 0;
996  int* GlobalViewIndices = 0;
997  double* MyViewValues = 0;
998  double* GlobalViewValues = 0;
999  int MaxNumIndices = A.Graph().MaxNumIndices();
1000  int* MyCopyIndices = new int[MaxNumIndices];
1001  int* GlobalCopyIndices = new int[MaxNumIndices];
1002  double* MyCopyValues = new double[MaxNumIndices];
1003  double* GlobalCopyValues = new double[MaxNumIndices];
1004 
1005  // Test query functions
1006 
1007  int NumMyRows = A.NumMyRows();
1008  if (verbose) cout << "\n\nNumber of local Rows = " << NumMyRows << endl<< endl;
1009 
1010  EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);
1011 
1012  int NumMyNonzeros = A.NumMyNonzeros();
1013  if (verbose) cout << "\n\nNumber of local Nonzero entries = " << NumMyNonzeros << endl<< endl;
1014 
1015  EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);
1016 
1017  int NumGlobalRows = A.NumGlobalRows();
1018  if (verbose) cout << "\n\nNumber of global Rows = " << NumGlobalRows << endl<< endl;
1019 
1020  EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);
1021 
1022  int NumGlobalNonzeros = A.NumGlobalNonzeros();
1023  if (verbose) cout << "\n\nNumber of global Nonzero entries = " << NumGlobalNonzeros << endl<< endl;
1024 
1025  EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);
1026 
1027  // GlobalRowView should be illegal (since we have local indices)
1028 
1029  EPETRA_TEST_ERR(!(A.ExtractGlobalRowView(A.RowMap().MaxMyGID(), NumGlobalIndices, GlobalViewValues, GlobalViewIndices)==-2),ierr);
1030 
1031  // Other binary tests
1032 
1033  EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
1034  EPETRA_TEST_ERR(!(A.Filled()),ierr);
1035  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID())),ierr);
1036  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID())),ierr);
1037  EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID()),ierr);
1038  EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID()),ierr);
1039  EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
1040  EPETRA_TEST_ERR(!(A.MyLRID(NumMyRows-1)),ierr);
1041  EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
1042  EPETRA_TEST_ERR(A.MyLRID(NumMyRows),ierr);
1043 
1044  forierr = 0;
1045  for (int i = 0; i < NumMyRows; i++) {
1046  int Row = A.GRID(i);
1047  A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyValues, GlobalCopyIndices);
1048  A.ExtractMyRowView(i, NumMyIndices, MyViewValues, MyViewIndices); // this is where the problem comes from
1049  forierr += !(NumGlobalIndices == NumMyIndices);
1050  for(int j = 1; j < NumMyIndices; j++) {
1051  forierr += !(MyViewIndices[j-1] < MyViewIndices[j]); // this is where the test fails
1052  }
1053  for(int j = 0; j < NumGlobalIndices; j++) {
1054  forierr += !(GlobalCopyIndices[j] == A.GCID(MyViewIndices[j]));
1055  forierr += !(A.LCID(GlobalCopyIndices[j]) == MyViewIndices[j]);
1056  forierr += !(GlobalCopyValues[j] == MyViewValues[j]);
1057  }
1058  }
1059  EPETRA_TEST_ERR(forierr,ierr);
1060 
1061  forierr = 0;
1062  for (int i = 0; i < NumMyRows; i++) {
1063  int Row = A.GRID(i);
1064  A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyValues, GlobalCopyIndices);
1065  A.ExtractMyRowCopy(i, MaxNumIndices, NumMyIndices, MyCopyValues, MyCopyIndices);
1066  forierr += !(NumGlobalIndices == NumMyIndices);
1067  for (int j = 1; j < NumMyIndices; j++)
1068  forierr += !(MyCopyIndices[j-1] < MyCopyIndices[j]);
1069  for (int j = 0; j < NumGlobalIndices; j++) {
1070  forierr += !(GlobalCopyIndices[j] == A.GCID(MyCopyIndices[j]));
1071  forierr += !(A.LCID(GlobalCopyIndices[j]) == MyCopyIndices[j]);
1072  forierr += !(GlobalCopyValues[j] == MyCopyValues[j]);
1073  }
1074 
1075  }
1076  EPETRA_TEST_ERR(forierr,ierr);
1077 
1078  delete [] MyCopyIndices;
1079  delete [] GlobalCopyIndices;
1080  delete [] MyCopyValues;
1081  delete [] GlobalCopyValues;
1082 
1083  if (verbose) cout << "\n\nRows sorted check OK" << endl<< endl;
1084 
1085  return (ierr);
1086 }
1087 
1089 {
1090  int numLocalElems = 5;
1091  int localProc = Comm.MyPID();
1092  int firstElem = localProc*numLocalElems;
1093  int err;
1094  Epetra_Map map(-1, numLocalElems, 0, Comm);
1095 
1096  Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, map, 1);
1097 
1098  for (int i=0; i<numLocalElems; ++i) {
1099  int row = firstElem+i;
1100  int col = row;
1101  double val = 1.0;
1102 
1103  err = A->InsertGlobalValues(row, 1, &val, &col);
1104  if (err != 0) {
1105  cerr << "A->InsertGlobalValues("<<row<<") returned err="<<err<<endl;
1106  return(err);
1107  }
1108  }
1109 
1110  A->FillComplete(false);
1111 
1112  Epetra_CrsMatrix B(Copy, A->Graph());
1113 
1114  delete A;
1115 
1116  for (int i=0; i<numLocalElems; ++i) {
1117  int row = firstElem+i;
1118  int col = row;
1119  double val = 1.0;
1120 
1121  err = B.ReplaceGlobalValues(row, 1, &val, &col);
1122  if (err != 0) {
1123  cerr << "B.InsertGlobalValues("<<row<<") returned err="<<err<<endl;
1124  return(err);
1125  }
1126  }
1127 
1128  return(0);
1129 }
1130 
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
int LRID(int GRID_in) const
Returns the local row index for given global row index, returns -1 if no local row for this global ro...
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
bool MyGRID(int GRID_in) const
Returns true if the GRID passed in belongs to the calling processor in this map, otherwise returns fa...
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
int Random()
Set multi-vector values to random numbers.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false...
int NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
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.
int NumGlobalRows() const
Returns the number of global matrix rows.
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.
double NormInf() const
Returns the infinity norm of the global matrix.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
int power_method(Epetra_CrsMatrix &A, double &lambda, int niters, double tolerance, bool verbose)
double NormOne() const
Returns the one norm of the global matrix.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
#define EPETRA_MIN(x, y)
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor&#39;...
std::string Epetra_Version()
bool IndicesAreLocal() const
If matrix indices has been transformed to local, this query returns true, otherwise it returns false...
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
bool UpperTriangular() const
If matrix is upper triangular in local index space, this query returns true, otherwise it returns fal...
int InvColSums(Epetra_Vector &x) const
Computes the inverse of the sum of absolute values of the columns of the Epetra_CrsMatrix, results returned in x.
virtual int MyPID() const =0
Return my process ID.
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:75
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
bool MyLRID(int LRID_in) const
Returns true if the LRID passed in belongs to the calling processor in this map, otherwise returns fa...
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Returns a view of the specified local row values via pointers to internal data.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int NumMyElements() const
Number of elements on the calling processor.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
Epetra_Map * RemoveEmptyProcesses() const
Return a new BlockMap with processes with zero elements removed.
Definition: Epetra_Map.cpp:179
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
int ExtractGlobalRowView(int GlobalRow, int &NumEntries, double *&Values, int *&Indices) const
Returns a view of the specified global row values via pointers to internal data.
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 OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor. ...
bool NoDiagonal() const
If matrix has no diagonal entries in global index space, this query returns true, otherwise it return...
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int RemoveEmptyProcessesInPlace(const Epetra_BlockMap *NewMap)
Remove processes owning zero rows from the Maps and their communicator.
int InvRowSums(Epetra_Vector &x) const
Computes the inverse of the sum of absolute values of the rows of the Epetra_CrsMatrix, results returned in x.
bool LowerTriangular() const
If matrix is lower triangular in local index space, this query returns true, otherwise it returns fal...
int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the left with a Epetra_Vector x.
int MaxMyGID() const
Returns the maximum global ID owned by this processor.
int MinMyGID() const
Returns the minimum global ID owned by this processor.
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
bool IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
int InvRowMaxs(Epetra_Vector &x) const
Computes the inverse of the max of absolute values of the rows of the Epetra_CrsMatrix, results returned in x.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int check_graph_sharing(Epetra_Comm &Comm)
void Barrier() const
Epetra_SerialComm Barrier function.
int GRID(int LRID_in) const
Returns the global row index for give local row index, returns IndexBase-1 if we don&#39;t have this loca...
int InvColMaxs(Epetra_Vector &x) const
Computes the max of absolute values of the columns of the Epetra_CrsMatrix, results returned in x...
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
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 MaxNumIndices() const
Returns the maximum number of nonzero entries across all rows on this processor.
int main(int argc, char *argv[])
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
int Broadcast(double *MyVals, int Count, int Root) const
Epetra_SerialComm Broadcast function.
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A &lt;- ScalarConstant * A)...
int InsertMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given local row of the matrix.
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
int RightScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the right with a Epetra_Vector x.
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Replace specified existing values with this list of entries for a given global row of the matrix...
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
Returns internal data pointers associated with Crs matrix format.
int GCID(int LCID_in) const
Returns the global column index for give local column index, returns IndexBase-1 if we don&#39;t have thi...
int NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor&#39;s portion of the matrix.