Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
59int check(Epetra_CrsMatrix& A, int NumMyRows1, long long NumGlobalRows1, int NumMyNonzeros1,
60 long long NumGlobalNonzeros1, long long * MyGlobalElements, bool verbose);
61
62int power_method(bool TransA, Epetra_CrsMatrix& A,
65 Epetra_Vector& resid,
66 double * lambda, int niters, double tolerance,
67 bool verbose);
68
70
71int 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;
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
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();
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);
736cout << xCol;
737 EPETRA_TEST_ERR(A2.RightScale(xCol),ierr);
738 float A2oneNormFloat = A2.NormOne();
739cout << 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
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*/
957return ierr ;
958}
959
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
991int 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
#define EPETRA_MIN(x, y)
@ View
@ Copy
std::string Epetra_Version()
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long MaxMyGID64() const
long long GID64(int LID) const
long long MinMyGID64() const
int NumMyElements() const
Number of elements on the calling processor.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual int MyPID() const =0
Return my process ID.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
double Flops() const
Returns the number of floating point operations with this multi-vector.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int MaxNumIndices() const
Returns the maximum number of nonzero entries across all rows on this processor.
long long NumGlobalNonzeros64() const
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
bool MyLRID(int LRID_in) const
Returns true if the LRID passed in belongs to the calling processor in this map, otherwise returns fa...
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
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 InvRowSums(Epetra_Vector &x) const
Computes the inverse of the sum of absolute values of the rows of the Epetra_CrsMatrix,...
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the left with a Epetra_Vector x.
int RightScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the right with a Epetra_Vector x.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
bool MyGRID(int GRID_in) const
Returns true if the GRID 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.
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
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.
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
long long NumGlobalRows64() const
long long NumGlobalNonzeros64() const
int InvColMaxs(Epetra_Vector &x) const
Computes the max of absolute values of the columns of the Epetra_CrsMatrix, results returned in x.
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...
double NormOne() const
Returns the one norm of the global matrix.
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
int InvRowMaxs(Epetra_Vector &x) const
Computes the inverse of the max of absolute values of the rows of the Epetra_CrsMatrix,...
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 NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor's portion of the matrix.
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A <- ScalarConstant * A).
int InvColSums(Epetra_Vector &x) const
Computes the inverse of the sum of absolute values of the columns of the Epetra_CrsMatrix,...
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Add this list of entries to existing values for a given global row of the matrix.
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
Returns internal data pointers associated with Crs matrix format.
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor'...
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false.
bool IndicesAreLocal() const
If matrix indices has been transformed to local, this query returns true, otherwise it returns false.
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.
int RemoveEmptyProcessesInPlace(const Epetra_BlockMap *NewMap)
Remove processes owning zero rows from the Maps and their communicator.
int NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
bool LowerTriangular() const
If matrix is lower triangular in local index space, this query returns true, otherwise it returns fal...
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
bool UpperTriangular() const
If matrix is upper triangular in local index space, this query returns true, otherwise it returns fal...
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
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.
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.
long long GRID64(int LRID_in) const
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.
bool NoDiagonal() const
If matrix has no diagonal entries in global index space, this query returns true, otherwise it return...
bool IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
double NormInf() const
Returns the infinity norm of the global matrix.
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
long long GCID64(int LCID_in) const
Epetra_Flops: The Epetra Floating Point Operations Class.
Definition: Epetra_Flops.h:58
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Definition: Epetra_Flops.h:77
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_Map * RemoveEmptyProcesses() const
Return a new BlockMap with processes with zero elements removed.
Definition: Epetra_Map.cpp:179
Epetra_MpiComm: The Epetra MPI Communication Class.
void Barrier() const
Epetra_MpiComm Barrier function.
int Broadcast(double *MyVals, int Count, int Root) const
Epetra_MpiComm Broadcast function.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
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 Random()
Set multi-vector values to random numbers.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Time: The Epetra Timing Class.
Definition: Epetra_Time.h:75
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
#define EPETRA_TEST_ERR(a, b)
int main(int argc, char *argv[])
int check(Epetra_CrsMatrix &A, int NumMyRows1, long long NumGlobalRows1, int NumMyNonzeros1, long long NumGlobalNonzeros1, long long *MyGlobalElements, bool verbose)
int check_graph_sharing(Epetra_Comm &Comm)
int power_method(bool TransA, Epetra_CrsMatrix &A, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, double *lambda, int niters, double tolerance, bool verbose)