Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/CrsMatrix_LL/cxx_main.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Epetra: Linear Algebra Services Package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
42 
43 #include "Epetra_Map.h"
44 #include "Epetra_Time.h"
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, long long NumGlobalRows1, int NumMyNonzeros1,
60  long long NumGlobalNonzeros1, long long * 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  long long 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, 0LL, Comm);
132 
133  // Get update list and number of local equations from newly created Map
134  long long* MyGlobalElements = new long long[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  long long* Indices = new long long[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  long long * Rowinds = new long long[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  long long 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((long long)-1, NumMyElements1, 0LL, Comm);
418 
419  // Get update list and number of local equations from newly created Map
420  long long* MyGlobalElements1 = new long long[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  long long *Indices1 = new long long[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((long long)-1,NumMyRows2,0LL,Comm);
529  Epetra_Map ColMap((long long)NumMyElements2,NumMyElements2,0LL,Comm);
530  // The DomainMap needs to be different from the ColMap for the test to
531  // be meaningful.
532  Epetra_Map DomainMap((long long)NumMyElements2,0LL,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((long long)-1,NumMyRangeElements2,0LL,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  long long RowLoopLength = RowMap.MaxMyGID64()-RowMap.MinMyGID64()+1;
557  for (long long i=0; i<RowLoopLength; i++)
558  RowLeftScaleValues[i] = (double) ((i + RowMap.MinMyGID64() ) % 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  long long RangeLoopLength = RangeMap.MaxMyGID64()-RangeMap.MinMyGID64()+1;
564  double * RangeLeftScaleValues = new double[(std::size_t) RangeLoopLength];
565  long long DomainLoopLength = DomainMap.MaxMyGID64()-DomainMap.MinMyGID64()+1;
566  double * DomainRightScaleValues = new double[(std::size_t) DomainLoopLength];
567  for (long long i=0; i<RangeLoopLength; i++)
568  RangeLeftScaleValues[i] = 1.0/((i + RangeMap.MinMyGID64() ) % 2 + 1);
569  for (long long i=0; i<DomainLoopLength;i++)
570  DomainRightScaleValues[i] = 1.0/((i + DomainMap.MinMyGID64() ) % 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((long long) 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  long long * myGlobalElements = new long long[NumMyRows3];
752  for (int i=0; i<NumMyRows3; i++) myGlobalElements[i] = MyPID+i;
753  Epetra_Map RowMap3((long long)NumProc*2, NumMyRows3, myGlobalElements, 0LL, Comm);
754  int NumMyElements3 = 5;
755  Epetra_CrsMatrix A3(Copy, RowMap3, NumMyElements3);
756  double * Values3 = new double[NumMyElements3];
757  long long * Indices3 = new long long[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((long long)NumProc+1, 0LL, Comm);
766  Epetra_Map DomainMap3((long long)NumMyElements3, 0LL, 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  long long * Indices3cm = new long long[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  long long nnz_A3cm = A3cm.Graph().NumGlobalNonzeros64();
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((long long) -1,NumMyRows, (long long) 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  long long GID = Map1.GID64(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  delete A1;
938  }
939 
940 
941 
942  delete [] Values2;
943  delete [] Indices2;
944  delete [] myGlobalElements;
945  delete [] Values3;
946  delete [] Indices3;
947  delete [] Values3cm;
948  delete [] Indices3cm;
949  delete [] RangeLeftScaleValues;
950  delete [] RowLeftScaleValues;
951 #ifdef EPETRA_MPI
952  MPI_Finalize() ;
953 #endif
954 
955 /* end main
956 */
957 return ierr ;
958 }
959 
960 int power_method(bool TransA, Epetra_CrsMatrix& A, Epetra_Vector& q, Epetra_Vector& z,
961  Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose)
962 {
963 
964  // Fill z with random Numbers
965  z.Random();
966 
967  // variable needed for iteration
968  double normz, residual;
969 
970  int ierr = 1;
971 
972  for(int iter = 0; iter < niters; iter++) {
973  z.Norm2(&normz); // Compute 2-norm of z
974  q.Scale(1.0/normz, z);
975  A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
976  q.Dot(z, lambda); // Approximate maximum eigenvaluE
977  if(iter%100==0 || iter+1==niters) {
978  resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
979  resid.Norm2(&residual);
980  if(verbose) cout << "Iter = " << iter << " Lambda = " << *lambda
981  << " Residual of A*q - lambda*q = " << residual << endl;
982  }
983  if(residual < tolerance) {
984  ierr = 0;
985  break;
986  }
987  }
988  return(ierr);
989 }
990 
991 int check(Epetra_CrsMatrix& A, int NumMyRows1, long long NumGlobalRows1, int NumMyNonzeros1,
992  long long NumGlobalNonzeros1, long long* MyGlobalElements, bool verbose)
993 {
994  (void)MyGlobalElements;
995  int ierr = 0, forierr = 0;
996  int NumGlobalIndices;
997  int NumMyIndices;
998  int* MyViewIndices = 0;
999  long long* GlobalViewIndices = 0;
1000  double* MyViewValues = 0;
1001  double* GlobalViewValues = 0;
1002  int MaxNumIndices = A.Graph().MaxNumIndices();
1003  int* MyCopyIndices = new int[MaxNumIndices];
1004  long long* GlobalCopyIndices = new long long[MaxNumIndices];
1005  double* MyCopyValues = new double[MaxNumIndices];
1006  double* GlobalCopyValues = new double[MaxNumIndices];
1007 
1008  // Test query functions
1009 
1010  int NumMyRows = A.NumMyRows();
1011  if (verbose) cout << "\n\nNumber of local Rows = " << NumMyRows << endl<< endl;
1012 
1013  EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);
1014 
1015  int NumMyNonzeros = A.NumMyNonzeros();
1016  if (verbose) cout << "\n\nNumber of local Nonzero entries = " << NumMyNonzeros << endl<< endl;
1017 
1018  EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);
1019 
1020  long long NumGlobalRows = A.NumGlobalRows64();
1021  if (verbose) cout << "\n\nNumber of global Rows = " << NumGlobalRows << endl<< endl;
1022 
1023  EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);
1024 
1025  long long NumGlobalNonzeros = A.NumGlobalNonzeros64();
1026  if (verbose) cout << "\n\nNumber of global Nonzero entries = " << NumGlobalNonzeros << endl<< endl;
1027 
1028  EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);
1029 
1030  // GlobalRowView should be illegal (since we have local indices)
1031 
1032  EPETRA_TEST_ERR(!(A.ExtractGlobalRowView(A.RowMap().MaxMyGID64(), NumGlobalIndices, GlobalViewValues, GlobalViewIndices)==-2),ierr);
1033 
1034  // Other binary tests
1035 
1036  EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
1037  EPETRA_TEST_ERR(!(A.Filled()),ierr);
1038  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID64())),ierr);
1039  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID64())),ierr);
1040  EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID64()),ierr);
1041  EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID64()),ierr);
1042  EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
1043  EPETRA_TEST_ERR(!(A.MyLRID(NumMyRows-1)),ierr);
1044  EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
1045  EPETRA_TEST_ERR(A.MyLRID(NumMyRows),ierr);
1046 
1047  forierr = 0;
1048  for (int i = 0; i < NumMyRows; i++) {
1049  long long Row = A.GRID64(i);
1050  A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyValues, GlobalCopyIndices);
1051  A.ExtractMyRowView(i, NumMyIndices, MyViewValues, MyViewIndices); // this is where the problem comes from
1052  forierr += !(NumGlobalIndices == NumMyIndices);
1053  for(int j = 1; j < NumMyIndices; j++) {
1054  forierr += !(MyViewIndices[j-1] < MyViewIndices[j]); // this is where the test fails
1055  }
1056  for(int j = 0; j < NumGlobalIndices; j++) {
1057  forierr += !(GlobalCopyIndices[j] == A.GCID64(MyViewIndices[j]));
1058  forierr += !(A.LCID(GlobalCopyIndices[j]) == MyViewIndices[j]);
1059  forierr += !(GlobalCopyValues[j] == MyViewValues[j]);
1060  }
1061  }
1062  EPETRA_TEST_ERR(forierr,ierr);
1063 
1064  forierr = 0;
1065  for (int i = 0; i < NumMyRows; i++) {
1066  long long Row = A.GRID64(i);
1067  A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyValues, GlobalCopyIndices);
1068  A.ExtractMyRowCopy(i, MaxNumIndices, NumMyIndices, MyCopyValues, MyCopyIndices);
1069  forierr += !(NumGlobalIndices == NumMyIndices);
1070  for (int j = 1; j < NumMyIndices; j++)
1071  forierr += !(MyCopyIndices[j-1] < MyCopyIndices[j]);
1072  for (int j = 0; j < NumGlobalIndices; j++) {
1073  forierr += !(GlobalCopyIndices[j] == A.GCID64(MyCopyIndices[j]));
1074  forierr += !(A.LCID(GlobalCopyIndices[j]) == MyCopyIndices[j]);
1075  forierr += !(GlobalCopyValues[j] == MyCopyValues[j]);
1076  }
1077 
1078  }
1079  EPETRA_TEST_ERR(forierr,ierr);
1080 
1081  delete [] MyCopyIndices;
1082  delete [] GlobalCopyIndices;
1083  delete [] MyCopyValues;
1084  delete [] GlobalCopyValues;
1085 
1086  if (verbose) cout << "\n\nRows sorted check OK" << endl<< endl;
1087 
1088  return (ierr);
1089 }
1090 
1092 {
1093  int numLocalElems = 5;
1094  int localProc = Comm.MyPID();
1095  long long firstElem = localProc*numLocalElems;
1096  int err;
1097  Epetra_Map map((long long)-1, numLocalElems, 0LL, Comm);
1098 
1099  Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, map, 1);
1100 
1101  for (int i=0; i<numLocalElems; ++i) {
1102  long long row = firstElem+i;
1103  long long col = row;
1104  double val = 1.0;
1105 
1106  err = A->InsertGlobalValues(row, 1, &val, &col);
1107  if (err != 0) {
1108  cerr << "A->InsertGlobalValues("<<row<<") returned err="<<err<<endl;
1109  return(err);
1110  }
1111  }
1112 
1113  A->FillComplete(false);
1114 
1115  Epetra_CrsMatrix B(Copy, A->Graph());
1116 
1117  delete A;
1118 
1119  for (int i=0; i<numLocalElems; ++i) {
1120  long long row = firstElem+i;
1121  long long col = row;
1122  double val = 1.0;
1123 
1124  err = B.ReplaceGlobalValues(row, 1, &val, &col);
1125  if (err != 0) {
1126  cerr << "B.ReplaceGlobalValues("<<row<<") returned err="<<err<<endl;
1127  return(err);
1128  }
1129  }
1130 
1131  return(0);
1132 }
1133 
1134 
long long MinMyGID64() const
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...
long long NumGlobalRows64() const
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.
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.
long long GRID64(int LRID_in) const
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.
long long NumGlobalNonzeros64() const
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.
long long GID64(int LID) const
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...
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 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.
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 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...
long long GCID64(int LCID_in) const
long long MaxMyGID64() const
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 NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor&#39;s portion of the matrix.