Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/ImportExport/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_IntVector.h"
48#include "Epetra_Import.h"
49#include "Epetra_Export.h"
50#include "Epetra_OffsetIndex.h"
51#ifdef EPETRA_MPI
52#include "Epetra_MpiComm.h"
53#include "mpi.h"
54#else
55#include "Epetra_SerialComm.h"
56#endif
57#ifndef __cplusplus
58#define __cplusplus
59#endif
60#include "../epetra_test_err.h"
61#include "Epetra_Version.h"
62
66
67int main(int argc, char *argv[])
68{
69 int ierr = 0, i, j, forierr = 0;
70
71#ifdef EPETRA_MPI
72 // Initialize MPI
73 MPI_Init(&argc,&argv);
74 Epetra_MpiComm Comm( MPI_COMM_WORLD );
75#else
77#endif
78
79 bool verbose = false;
80
81 // Check if we should print results to standard out
82 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
83
84
85
86
87 //char tmp;
88 //if (Comm.MyPID()==0) cout << "Press any key to continue..."<< endl;
89 //if (Comm.MyPID()==0) cin >> tmp;
90 //Comm.Barrier();
91
92 Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
93 int MyPID = Comm.MyPID();
94 int NumProc = Comm.NumProc();
95
96 if (verbose && MyPID==0)
97 cout << Epetra_Version() << endl << endl;
98
99 if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
100 << " is alive."<<endl;
101
102 // Redefine verbose to only print on PE 0
103 if (verbose && Comm.MyPID()!=0) verbose = false;
104
105 int NumMyEquations = 20;
106 int NumGlobalEquations = NumMyEquations*NumProc+EPETRA_MIN(NumProc,3);
107 if (MyPID < 3) NumMyEquations++;
108 // Construct a Source Map that puts approximately the same Number of equations on each processor in
109 // uniform global ordering
110
111 Epetra_Map SourceMap(NumGlobalEquations, NumMyEquations, 0, Comm);
112
113 // Get update list and number of local equations from newly created Map
114 int NumMyElements = SourceMap.NumMyElements();
115 int * SourceMyGlobalElements = new int[NumMyElements];
116 SourceMap.MyGlobalElements(SourceMyGlobalElements);
117
118 // Construct a Target Map that will contain:
119 // some unchanged elements (relative to the soure map),
120 // some permuted elements
121 // some off-processor elements
122 Epetra_Vector RandVec(SourceMap);
123 RandVec.Random(); // This creates a vector of random numbers between negative one and one.
124
125 int *TargetMyGlobalElements = new int[NumMyElements];
126
127 int MinGID = SourceMap.MinMyGID();
128 for (i=0; i< NumMyEquations/2; i++) TargetMyGlobalElements[i] = i + MinGID; // Half will be the same...
129 for (i=NumMyEquations/2; i<NumMyEquations; i++) {
130 int index = abs((int)(((double) (NumGlobalEquations-1) ) * RandVec[i]));
131 TargetMyGlobalElements[i] = EPETRA_MIN(NumGlobalEquations-1,EPETRA_MAX(0,index));
132 }
133
134 int NumSameIDs = 0;
135 int NumPermutedIDs = 0;
136 int NumRemoteIDs = 0;
137 bool StillContiguous = true;
138 for (i=0; i < NumMyEquations; i++) {
139 if (SourceMyGlobalElements[i]==TargetMyGlobalElements[i] && StillContiguous)
140 NumSameIDs++;
141 else if (SourceMap.MyGID(TargetMyGlobalElements[i])) {
142 StillContiguous = false;
143 NumPermutedIDs++;
144 }
145 else {
146 StillContiguous = false;
147 NumRemoteIDs++;
148 }
149 }
150 EPETRA_TEST_ERR(!(NumMyEquations==NumSameIDs+NumPermutedIDs+NumRemoteIDs),ierr);
151
152 Epetra_Map TargetMap(-1, NumMyElements, TargetMyGlobalElements, 0, Comm);
153
154 // Create a multivector whose elements are GlobalID * (column number +1)
155
156 int NumVectors = 3;
157 Epetra_MultiVector SourceMultiVector(SourceMap, NumVectors);
158 for (j=0; j < NumVectors; j++)
159 for (i=0; i < NumMyElements; i++)
160 SourceMultiVector[j][i] = (double) SourceMyGlobalElements[i]*(j+1);
161
162 // Create a target multivector that we will fill using an Import
163
164 Epetra_MultiVector TargetMultiVector(TargetMap, NumVectors);
165
166 Epetra_Import Importer(TargetMap, SourceMap);
167
168 EPETRA_TEST_ERR(!(TargetMultiVector.Import(SourceMultiVector, Importer, Insert)==0),ierr);
169
170 // Test Target against expected values
171 forierr = 0;
172 for (j=0; j < NumVectors; j++)
173 for (i=0; i < NumMyElements; i++) {
174 if (TargetMultiVector[j][i]!= (double) TargetMyGlobalElements[i]*(j+1))
175 cout << "TargetMultiVector["<<i<<"]["<<j<<"] = " << TargetMultiVector[j][i]
176 << " TargetMyGlobalElements[i]*(j+1) = " << TargetMyGlobalElements[i]*(j+1) << endl;
177 forierr += !(TargetMultiVector[j][i]== (double) TargetMyGlobalElements[i]*(j+1));
178 }
179 EPETRA_TEST_ERR(forierr,ierr);
180
181 if (verbose) cout << "MultiVector Import using Importer Check OK" << endl << endl;
182
183
185
186 // Now use Importer to do an export
187
188 Epetra_Vector TargetVector(SourceMap);
189 Epetra_Vector ExpectedTarget(SourceMap);
190 Epetra_Vector SourceVector(TargetMap);
191
192 NumSameIDs = Importer.NumSameIDs();
193 int NumPermuteIDs = Importer.NumPermuteIDs();
194 int NumExportIDs = Importer.NumExportIDs();
195 int *PermuteFromLIDs = Importer.PermuteFromLIDs();
196 int *ExportLIDs = Importer.ExportLIDs();
197 int *ExportPIDs = Importer.ExportPIDs();
198
199 for (i=0; i < NumSameIDs; i++) ExpectedTarget[i] = (double) (MyPID+1);
200 for (i=0; i < NumPermuteIDs; i++) ExpectedTarget[PermuteFromLIDs[i]] =
201 (double) (MyPID+1);
202 for (i=0; i < NumExportIDs; i++) ExpectedTarget[ExportLIDs[i]] +=
203 (double) (ExportPIDs[i]+1);
204
205 for (i=0; i < NumMyElements; i++) SourceVector[i] = (double) (MyPID+1);
206
207 EPETRA_TEST_ERR(!(TargetVector.Export(SourceVector, Importer, Add)==0),ierr);
208
209 forierr = 0;
210 for (i=0; i < NumMyElements; i++) {
211 if (TargetVector[i]!= ExpectedTarget[i])
212 cout << " TargetVector["<<i<<"] = " << TargetVector[i]
213 << " ExpectedTarget["<<i<<"] = " << ExpectedTarget[i] << " on PE " << MyPID << endl;
214 forierr += !(TargetVector[i]== ExpectedTarget[i]);
215 }
216 EPETRA_TEST_ERR(forierr,ierr);
217
218 if (verbose) cout << "Vector Export using Importer Check OK" << endl << endl;
219
221 // Now use Importer to create a reverse exporter
222 TargetVector.PutScalar(0.0);
223 Epetra_Export ReversedImport(Importer);
224
225 EPETRA_TEST_ERR(!(TargetVector.Export(SourceVector, ReversedImport, Add)==0),ierr);
226
227 forierr = 0;
228 for (i=0; i < NumMyElements; i++) {
229 if (TargetVector[i]!= ExpectedTarget[i])
230 cout << " TargetVector["<<i<<"] = " << TargetVector[i]
231 << " ExpectedTarget["<<i<<"] = " << ExpectedTarget[i] << " on PE " << MyPID << endl;
232 forierr += !(TargetVector[i]== ExpectedTarget[i]);
233 }
234 EPETRA_TEST_ERR(forierr,ierr);
235
236 if (verbose) cout << "Vector Export using Reversed Importer Check OK" << endl << endl;
237
239 // Now use Exporter to create a reverse importer
240 TargetVector.PutScalar(0.0);
241 Epetra_Import ReversedExport(ReversedImport);
242
243 EPETRA_TEST_ERR(!(TargetVector.Export(SourceVector, ReversedExport, Add)==0),ierr);
244
245 forierr = 0;
246 for (i=0; i < NumMyElements; i++) {
247 if (TargetVector[i]!= ExpectedTarget[i])
248 cout << " TargetVector["<<i<<"] = " << TargetVector[i]
249 << " ExpectedTarget["<<i<<"] = " << ExpectedTarget[i] << " on PE " << MyPID << endl;
250 forierr += !(TargetVector[i]== ExpectedTarget[i]);
251 }
252 EPETRA_TEST_ERR(forierr,ierr);
253
254 if (verbose) cout << "Vector Export using Reversed Exporter Check OK" << endl << endl;
255
257 // Build a tridiagonal system two ways:
258 // 1) From "standard" matrix view where equations are uniquely owned.
259 // 2) From 1D PDE view where nodes (equations) between processors are shared and partial contributions are done
260 // in parallel, then merged together at the end of the construction process.
261 //
263
264
265
266 // Construct a Standard Map that puts approximately the same number of equations on each processor in
267 // uniform global ordering
268
269 Epetra_Map StandardMap(NumGlobalEquations, NumMyEquations, 0, Comm);
270
271 // Get update list and number of local equations from newly created Map
272 NumMyElements = StandardMap.NumMyElements();
273 int * StandardMyGlobalElements = new int[NumMyElements];
274 StandardMap.MyGlobalElements(StandardMyGlobalElements);
275
276
277 // Create a standard Epetra_CrsGraph
278
279 Epetra_CrsGraph StandardGraph(Copy, StandardMap, 3);
280 EPETRA_TEST_ERR(StandardGraph.IndicesAreGlobal(),ierr);
281 EPETRA_TEST_ERR(StandardGraph.IndicesAreLocal(),ierr);
282
283 // Add rows one-at-a-time
284 // Need some vectors to help
285 // Off diagonal Values will always be -1
286
287
288 int *Indices = new int[2];
289 int NumEntries;
290
291 forierr = 0;
292 for (i=0; i<NumMyEquations; i++)
293 {
294 if (StandardMyGlobalElements[i]==0)
295 {
296 Indices[0] = 1;
297 NumEntries = 1;
298 }
299 else if (StandardMyGlobalElements[i] == NumGlobalEquations-1)
300 {
301 Indices[0] = NumGlobalEquations-2;
302 NumEntries = 1;
303 }
304 else
305 {
306 Indices[0] = StandardMyGlobalElements[i]-1;
307 Indices[1] = StandardMyGlobalElements[i]+1;
308 NumEntries = 2;
309 }
310 forierr += !(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], NumEntries, Indices)==0);
311 forierr += !(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], 1, StandardMyGlobalElements+i)==0); // Put in the diagonal entry
312 }
313 EPETRA_TEST_ERR(forierr,ierr);
314
315 // Finish up
316 EPETRA_TEST_ERR(!(StandardGraph.IndicesAreGlobal()),ierr);
317 EPETRA_TEST_ERR(!(StandardGraph.FillComplete()==0),ierr);
318 EPETRA_TEST_ERR(!(StandardGraph.IndicesAreLocal()),ierr);
319 EPETRA_TEST_ERR(StandardGraph.StorageOptimized(),ierr);
320 StandardGraph.OptimizeStorage();
321 EPETRA_TEST_ERR(!(StandardGraph.StorageOptimized()),ierr);
322 EPETRA_TEST_ERR(StandardGraph.UpperTriangular(),ierr);
323 EPETRA_TEST_ERR(StandardGraph.LowerTriangular(),ierr);
324
325 // Create Epetra_CrsMatrix using the just-built graph
326
327 Epetra_CrsMatrix StandardMatrix(Copy, StandardGraph);
328 EPETRA_TEST_ERR(StandardMatrix.IndicesAreGlobal(),ierr);
329 EPETRA_TEST_ERR(!(StandardMatrix.IndicesAreLocal()),ierr);
330
331 // Add rows one-at-a-time
332 // Need some vectors to help
333 // Off diagonal Values will always be -1
334
335
336 double *Values = new double[2];
337 Values[0] = -1.0; Values[1] = -1.0;
338 double two = 2.0;
339
340 forierr = 0;
341 for (i=0; i<NumMyEquations; i++)
342 {
343 if (StandardMyGlobalElements[i]==0)
344 {
345 Indices[0] = 1;
346 NumEntries = 1;
347 }
348 else if (StandardMyGlobalElements[i] == NumGlobalEquations-1)
349 {
350 Indices[0] = NumGlobalEquations-2;
351 NumEntries = 1;
352 }
353 else
354 {
355 Indices[0] = StandardMyGlobalElements[i]-1;
356 Indices[1] = StandardMyGlobalElements[i]+1;
357 NumEntries = 2;
358 }
359 forierr += !(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], NumEntries, Values, Indices)==0);
360 // Put in the diagonal entry
361 forierr += !(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], 1, &two, StandardMyGlobalElements+i)==0);
362 }
363 EPETRA_TEST_ERR(forierr,ierr);
364
365 // Finish up
366 EPETRA_TEST_ERR(!(StandardMatrix.IndicesAreLocal()),ierr);
367 EPETRA_TEST_ERR(!(StandardMatrix.FillComplete()==0),ierr);
368 EPETRA_TEST_ERR(!(StandardMatrix.IndicesAreLocal()),ierr);
369 // EPETRA_TEST_ERR((StandardMatrix.StorageOptimized()),ierr);
370 EPETRA_TEST_ERR((StandardMatrix.OptimizeStorage()),ierr);
371 EPETRA_TEST_ERR(!(StandardMatrix.StorageOptimized()),ierr);
372 EPETRA_TEST_ERR(StandardMatrix.UpperTriangular(),ierr);
373 EPETRA_TEST_ERR(StandardMatrix.LowerTriangular(),ierr);
374
375 // Construct an Overlapped Map of StandardMap that include the endpoints from two neighboring processors.
376
377 int OverlapNumMyElements;
378 int OverlapMinMyGID;
379
380 OverlapNumMyElements = NumMyElements + 1;
381 if (MyPID==0) OverlapNumMyElements--;
382
383 if (MyPID==0) OverlapMinMyGID = StandardMap.MinMyGID();
384 else OverlapMinMyGID = StandardMap.MinMyGID()-1;
385
386 int * OverlapMyGlobalElements = new int[OverlapNumMyElements];
387
388 for (i=0; i< OverlapNumMyElements; i++) OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
389
390 Epetra_Map OverlapMap(-1, OverlapNumMyElements, OverlapMyGlobalElements, 0, Comm);
391
392 // Create the Overlap Epetra_Matrix
393
394 Epetra_CrsMatrix OverlapMatrix(Copy, OverlapMap, 4);
395 EPETRA_TEST_ERR(OverlapMatrix.IndicesAreGlobal(),ierr);
396 EPETRA_TEST_ERR(OverlapMatrix.IndicesAreLocal(),ierr);
397
398 // Add matrix element one cell at a time.
399 // Each cell does an incoming and outgoing flux calculation
400
401
402 double pos_one = 1.0;
403 double neg_one = -1.0;
404
405 forierr = 0;
406 for (i=0; i<OverlapNumMyElements; i++)
407 {
408 int node_left = OverlapMyGlobalElements[i]-1;
409 int node_center = node_left + 1;
410 int node_right = node_left + 2;
411 if (i>0) {
412 if (node_left>-1)
413 forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_left)==0);
414 forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
415 }
416 if (i<OverlapNumMyElements-1) {
417 forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
418 if (node_right<NumGlobalEquations)
419 forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_right)==0);
420 }
421 }
422 EPETRA_TEST_ERR(forierr,ierr);
423
424 // Handle endpoints
425 if (MyPID==0) {
426 int node_center = 0;
427 EPETRA_TEST_ERR(!(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0),ierr);
428 }
429 if (MyPID==NumProc-1) {
430 int node_center = OverlapMyGlobalElements[OverlapNumMyElements-1];
431 EPETRA_TEST_ERR(!(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0),ierr);
432 }
433
434 EPETRA_TEST_ERR(!(OverlapMatrix.FillComplete()==0),ierr);
435
436 // Make a gathered matrix from OverlapMatrix. It should be identical to StandardMatrix
437
438 Epetra_CrsMatrix GatheredMatrix(Copy, StandardGraph);
439 Epetra_Export Exporter(OverlapMap, StandardMap);
440 EPETRA_TEST_ERR(!(GatheredMatrix.Export(OverlapMatrix, Exporter, Add)==0),ierr);
441 EPETRA_TEST_ERR(!(GatheredMatrix.FillComplete()==0),ierr);
442
443 // Check if entries of StandardMatrix and GatheredMatrix are identical
444
445 int StandardNumEntries, GatheredNumEntries;
446 int * StandardIndices, * GatheredIndices;
447 double * StandardValues, * GatheredValues;
448
449 int StandardNumMyNonzeros = StandardMatrix.NumMyNonzeros();
450 int GatheredNumMyNonzeros = GatheredMatrix.NumMyNonzeros();
451 EPETRA_TEST_ERR(!(StandardNumMyNonzeros==GatheredNumMyNonzeros),ierr);
452
453 int StandardNumMyRows = StandardMatrix.NumMyRows();
454 int GatheredNumMyRows = GatheredMatrix.NumMyRows();
455 EPETRA_TEST_ERR(!(StandardNumMyRows==GatheredNumMyRows),ierr);
456
457 forierr = 0;
458 for (i=0; i< StandardNumMyRows; i++)
459 {
460 forierr += !(StandardMatrix.ExtractMyRowView(i, StandardNumEntries, StandardValues, StandardIndices)==0);
461 forierr += !(GatheredMatrix.ExtractMyRowView(i, GatheredNumEntries, GatheredValues, GatheredIndices)==0);
462 forierr += !(StandardNumEntries==GatheredNumEntries);
463 for (j=0; j < StandardNumEntries; j++) {
464 //if (StandardIndices[j]!=GatheredIndices[j])
465 // cout << "MyPID = " << MyPID << " i = " << i << " StandardIndices[" << j << "] = " << StandardIndices[j]
466 // << " GatheredIndices[" << j << "] = " << GatheredIndices[j] << endl;
467 //if (StandardValues[j]!=GatheredValues[j])
468 //cout << "MyPID = " << MyPID << " i = " << i << " StandardValues[" << j << "] = " << StandardValues[j]
469 // << " GatheredValues[" << j << "] = " << GatheredValues[j] << endl;
470 forierr += !(StandardIndices[j]==GatheredIndices[j]);
471 forierr += !(StandardValues[j]==GatheredValues[j]);
472 }
473 }
474 EPETRA_TEST_ERR(forierr,ierr);
475
476 if (verbose) cout << "Matrix Export Check OK" << endl << endl;
477
478 //Do Again with use of Epetra_OffsetIndex object for speed
479 Epetra_OffsetIndex OffsetIndex( OverlapMatrix.Graph(), GatheredMatrix.Graph(), Exporter );
480 EPETRA_TEST_ERR(!(GatheredMatrix.Export(OverlapMatrix, Exporter, Add)==0),ierr);
481
482 if (verbose) cout << "Optimized Matrix Export Check OK" << endl << endl;
483
484 bool passed;
485 Epetra_IntVector v1(StandardMap); v1.PutValue(2);
486 Epetra_IntVector v2(StandardMap); v2.PutValue(3);
487
488 Epetra_Export identExporter(StandardMap,StandardMap); // Identity exporter
489 EPETRA_TEST_ERR(!(v2.Export(v1, identExporter, Insert)==0),ierr);
490 passed = (v2.MinValue()==2);
491 EPETRA_TEST_ERR(!passed,ierr);
492
493 v1.PutValue(1);
494 Epetra_Import identImporter(StandardMap,StandardMap); // Identity importer
495 EPETRA_TEST_ERR(!(v2.Import(v1, identExporter, Insert)==0),ierr);
496 passed = passed && (v2.MaxValue()==1);
497 EPETRA_TEST_ERR(!passed,ierr);
498
499 if (verbose) {
500 if (passed) cout << "Identity Import/Export Check OK" << endl << endl;
501 else cout << "Identity Import/Export Check Failed" << endl << endl;
502 }
503
504 int NumSubMapElements = StandardMap.NumMyElements()/2;
505 int SubStart = Comm.MyPID();
506 NumSubMapElements = EPETRA_MIN(NumSubMapElements,StandardMap.NumMyElements()-SubStart);
507 Epetra_Map SubMap(-1, NumSubMapElements, StandardMyGlobalElements+SubStart, 0, Comm);
508
509 Epetra_IntVector v3(View, SubMap, SubMap.MyGlobalElements()); // Fill v3 with GID values for variety
510 Epetra_Export subExporter(SubMap, StandardMap); // Export to a subset of indices of standard map
511 EPETRA_TEST_ERR(!(v2.Export(v3,subExporter,Insert)==0),ierr);
512
513 forierr = 0;
514 for (i=0; i<SubMap.NumMyElements(); i++) {
515 int i1 = StandardMap.LID(SubMap.GID(i));
516 forierr += !(v3[i]==v2[i1]);
517 }
518 EPETRA_TEST_ERR(forierr,ierr);
519
520 Epetra_Import subImporter(StandardMap, SubMap); // Import to a subset of indices of standard map
521 EPETRA_TEST_ERR(!(v1.Import(v3,subImporter,Insert)==0),ierr);
522
523 for (i=0; i<SubMap.NumMyElements(); i++) {
524 int i1 = StandardMap.LID(SubMap.GID(i));
525 forierr += !(v3[i]==v1[i1]);
526 }
527 EPETRA_TEST_ERR(forierr,ierr);
528
529 if (verbose) {
530 if (forierr==0) cout << "SubMap Import/Export Check OK" << endl << endl;
531 else cout << "SubMap Import/Export Check Failed" << endl << endl;
532 }
533
534
535#ifdef DOESNT_WORK_IN_PARALLEL
536 forierr = special_submap_import_test(Comm);
537 EPETRA_TEST_ERR(forierr, ierr);
538
539 if (verbose) {
540 if (forierr==0) cout << "Special SubMap Import Check OK" << endl << endl;
541 else cout << "Special SubMap Import Check Failed" << endl << endl;
542 }
543#endif
544
545
546 forierr = alternate_import_constructor_test(Comm);
547 EPETRA_TEST_ERR(forierr, ierr);
548
549 if (verbose) {
550 if (forierr==0) cout << "Alternative Import Constructor Check OK" << endl << endl;
551 else cout << "Alternative Import Constructor Check Failed" << endl << endl;
552 }
553
554
555 // Now let's test Importer replacement: Coalesce to 1 proc...
556 Epetra_CrsMatrix FunMatrix(StandardMatrix);
557 if(Comm.NumProc()!=1) {
558 forierr=0;
559 long long num_global_elements1 = FunMatrix.DomainMap().NumGlobalElements64();
560 long long num_global_elements2 = FunMatrix.DomainMap().MaxAllGID64()- FunMatrix.DomainMap().MinAllGID64()+1;
561 if(num_global_elements1 == num_global_elements2) {
562 // The original domain map is linear. Let's have fun
563 int NumMyElements = Comm.MyPID()==0 ? num_global_elements1 : 0;
564 if(FunMatrix.DomainMap().GlobalIndicesLongLong()) {
565 Epetra_Map NewMap((long long)-1,NumMyElements,(long long)0,Comm);
566 Epetra_Import NewImport(FunMatrix.ColMap(),NewMap);
567 FunMatrix.ReplaceDomainMapAndImporter(NewMap,&NewImport);
568 }
569 else {
570 Epetra_Map NewMap(-1,NumMyElements,0,Comm);
571 Epetra_Import NewImport(FunMatrix.ColMap(),NewMap);
572 FunMatrix.ReplaceDomainMapAndImporter(NewMap,&NewImport);
573 }
574
575 // Now let's test the new importer...
576
577 // Fill a random vector on the original map
578 Epetra_Vector OriginalVec(StandardMatrix.DomainMap());
579 Epetra_Vector OriginalY(FunMatrix.RangeMap(),true);
580 OriginalVec.SetSeed(24601);
581 OriginalVec.Random();
582
583 // Move said random vector to a single proc
584 Epetra_Vector NewVec(FunMatrix.DomainMap(),true);
585 Epetra_Vector NewY(FunMatrix.RangeMap(),true);
586 Epetra_Import ImportOld2New(FunMatrix.DomainMap(),StandardMatrix.DomainMap());
587 NewVec.Import(OriginalVec,ImportOld2New,Add);
588
589 // Test the Importer Copy Constructor
590 Epetra_Vector ColVec1(FunMatrix.ColMap(),true);
591 Epetra_Import ColImport(FunMatrix.ColMap(),FunMatrix.DomainMap());
592 ColVec1.Import(NewVec,ColImport,Add);
593
594 Epetra_Vector ColVec2(FunMatrix.ColMap(),true);
595 Epetra_Import ColImport2(ColImport);
596 ColVec2.Import(NewVec,ColImport2,Add);
597
598 double norm;
599 ColVec1.Update(-1.0,ColVec2,1.0);
600 NewY.Norm2(&norm);
601 if(norm > 1e-12) forierr=-1;
602 if (verbose) {
603 if (forierr==0) cout << "Import Copy Constructor Check OK" << endl << endl;
604 else cout << "Import Copy Constructor Check Failed" << endl << endl;
605 }
606
607 // Test replaceDomainMapAndImporter
608 // Now do two multiplies and compare
609 StandardMatrix.Apply(OriginalVec,OriginalY);
610 FunMatrix.Apply(NewVec,NewY);
611 NewY.Update(-1.0,OriginalY,1.0);
612 NewY.Norm2(&norm);
613 if(norm > 1e-12) forierr=-1;
614 if (verbose) {
615 if (forierr==0) cout << "ReplaceDomainMapAndImporter Check OK" << endl << endl;
616 else cout << "ReplaceDomainMapAndImporter Check Failed" << endl << endl;
617 }
618 }
619 }
620
621
622
623 // Release all objects
624
625 delete [] SourceMyGlobalElements;
626 delete [] TargetMyGlobalElements;
627 delete [] OverlapMyGlobalElements;
628 delete [] StandardMyGlobalElements;
629
630 delete [] Values;
631 delete [] Indices;
632
633#ifdef EPETRA_MPI
634 MPI_Finalize() ;
635#endif
636
637/* end main
638*/
639return ierr ;
640}
641
643{
644 int localProc = Comm.MyPID();
645
646 //set up ids_source and ids_target such that ids_source are only
647 //a subset of ids_target, and furthermore that ids_target are ordered
648 //such that the LIDs don't match up. In other words, even if gid 2 does
649 //exist in both ids_source and ids_target, it will correspond to different
650 //LIDs on at least 1 proc.
651 //
652 //This is to test a certain bug-fix in Epetra_Import where the 'RemoteLIDs'
653 //array wasn't being calculated correctly on all procs.
654
655 int ids_source[1];
656 ids_source[0] = localProc*2+2;
657
658 int ids_target[3];
659 ids_target[0] = localProc*2+2;
660 ids_target[1] = localProc*2+1;
661 ids_target[2] = localProc*2+0;
662
663 Epetra_Map map_source(-1, 1, &ids_source[0], 0, Comm);
664 Epetra_Map map_target(-1, 3, &ids_target[0], 0, Comm);
665
666 Epetra_Import importer(map_target, map_source);
667
668 Epetra_IntVector vec_source(map_source);
669 Epetra_IntVector vec_target(map_target);
670
671 vec_target.PutValue(0);
672
673 //set vec_source's contents so that entry[i] == GID[i].
674 int* GIDs = map_source.MyGlobalElements();
675 for(int i=0; i<map_source.NumMyElements(); ++i) {
676 vec_source[i] = GIDs[i];
677 }
678
679 //Import vec_source into vec_target. This should result in the contents
680 //of vec_target remaining 0 for the entries that don't exist in vec_source,
681 //and other entries should be equal to the corresponding GID in the map.
682
683 vec_target.Import(vec_source, importer, Insert);
684
685 GIDs = map_target.MyGlobalElements();
686 int test_failed = 0;
687
688 //the test passes if the i-th entry in vec_target equals either 0 or
689 //GIDs[i].
690 for(int i=0; i<vec_target.MyLength(); ++i) {
691 if (vec_target[i] != GIDs[i] && vec_target[i] != 0) test_failed = 1;
692 }
693
694 int global_result;
695 Comm.MaxAll(&test_failed, &global_result, 1);
696
697 //If test didn't fail on any procs, global_result should be 0.
698 //If test failed on any proc, global_result should be 1.
699 return global_result;
700}
701
703{
704 int localProc = Comm.MyPID();
705
706
707 int ids_source[1];
708 ids_source[0] = localProc*2+2;
709
710 int ids_target[3];
711 ids_target[0] = localProc*2+2;
712 ids_target[1] = localProc*2+1;
713 ids_target[2] = localProc*2+0;
714
715 Epetra_Map map_source(-1, 1, &ids_source[0], 0, Comm);
716 Epetra_Map map_target(-1, 3, &ids_target[0], 0, Comm);
717
718 Epetra_Import importer(map_target, map_source);
719
720 Epetra_IntVector vec_source(map_source);
721 Epetra_IntVector vec_target(map_target);
722
723 vec_target.PutValue(0);
724
725 //set vec_source's contents so that entry[i] == GID[i].
726 int* GIDs = map_source.MyGlobalElements();
727 for(int i=0; i<map_source.NumMyElements(); ++i) {
728 vec_source[i] = GIDs[i];
729 }
730
731 //Import vec_source into vec_target. This should result in the contents
732 //of vec_target remaining 0 for the entries that don't exist in vec_source,
733 //and other entries should be equal to the corresponding GID in the map.
734
735 vec_target.Import(vec_source, importer, Insert);
736
737 GIDs = map_target.MyGlobalElements();
738 int test_failed = 0;
739
740 //the test passes if the i-th entry in vec_target equals either 0 or
741 //GIDs[i].
742 for(int i=0; i<vec_target.MyLength(); ++i) {
743 if (vec_target[i] != GIDs[i] && vec_target[i] != 0) test_failed = 1;
744 }
745
746 int global_result;
747 Comm.MaxAll(&test_failed, &global_result, 1);
748
749 //If test didn't fail on any procs, global_result should be 0.
750 //If test failed on any proc, global_result should be 1.
751 return global_result;
752}
753
754
755int test_import_gid(const char * name,Epetra_IntVector & Source, Epetra_IntVector & Target, const Epetra_Import & Import){
756 int i;
757 bool test_passed=true;
758
759 // Setup
760 for(i=0; i<Source.MyLength(); i++)
761 Source[i] = Source.Map().GID(i);
762 Target.PutValue(0);
763
764 // Import
765 Target.Import(Source,Import,Add);
766
767 // Test
768 for(i=0; i<Target.MyLength(); i++){
769 if(Target[i] != Target.Map().GID(i)) test_passed=false;
770 }
771
772 if(!test_passed){
773 printf("[%d] test_import_gid %s failed: ",Source.Map().Comm().MyPID(),name);
774 for(i=0; i<Target.MyLength(); i++)
775 printf("%2d(%2d) ",Target[i],Target.Map().GID(i));
776 printf("\n");
777 fflush(stdout);
778 }
779
780 return !test_passed;
781}
782
783
784
786 int rv=0;
787 int nodes_per_proc=10;
788 int numprocs = Comm.NumProc();
789 int mypid = Comm.MyPID();
790
791 // Only run if we have multiple procs & MPI
792 if(numprocs==0) return 0;
793#ifndef HAVE_MPI
794 return 0;
795#endif
796
797 // Build Map 1 - linear
798 Epetra_Map Map1(-1,nodes_per_proc,0,Comm);
799
800 // Build Map 2 - mod striped
801 std::vector<int> MyGIDs(nodes_per_proc);
802 for(int i=0; i<nodes_per_proc; i++)
803 MyGIDs[i] = (mypid*nodes_per_proc + i) % numprocs;
804 Epetra_Map Map2(-1,nodes_per_proc,&MyGIDs[0],0,Comm);
805
806 // For testing
807 Epetra_IntVector Source(Map1), Target(Map2);
808
809
810 // Build Import 1 - normal
811 Epetra_Import Import1(Map2,Map1);
812 rv = rv|| test_import_gid("Alt test: 2 map constructor",Source,Target, Import1);
813
814 // Build Import 2 - no-comm constructor
815 int Nremote=Import1.NumRemoteIDs();
816 const int * RemoteLIDs = Import1.RemoteLIDs();
817 std::vector<int> RemotePIDs(Nremote+1); // I hate you, stl vector....
818 std::vector<int> AllPIDs;
819 Epetra_Util::GetPids(Import1,AllPIDs,true);
820
821 for(int i=0; i<Nremote; i++) {
822 RemotePIDs[i]=AllPIDs[RemoteLIDs[i]];
823 }
824 Epetra_Import Import2(Import1.TargetMap(),Import1.SourceMap(),Nremote,&RemotePIDs[0],Import1.NumExportIDs(),Import1.ExportLIDs(),Import1.ExportPIDs());
825
826 rv = rv || test_import_gid("Alt test: no comm constructor",Source,Target,Import2);
827
828
829 // Build Import 3 - Remotes only
830 Epetra_Import Import3(Import1.TargetMap(),Import1.SourceMap(),Nremote,&RemotePIDs[0]);
831 rv = rv || test_import_gid("Alt test: remote only constructor",Source,Target, Import3);
832
833
834 return rv;
835}
#define EPETRA_MIN(x, y)
#define EPETRA_MAX(x, y)
std::string Epetra_Version()
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long MaxAllGID64() const
int MinMyGID() const
Returns the minimum global ID owned by this processor.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
int GID(int LID) const
Returns global ID of local ID, return IndexBase-1 if not found on this processor.
long long NumGlobalElements64() const
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
long long MinAllGID64() const
int NumMyElements() const
Number of elements on the calling processor.
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition Epetra_Comm.h:73
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
virtual int NumProc() const =0
Returns total number of processes.
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int InsertGlobalIndices(int GlobalRow, int NumIndices, int *Indices)
Enter a list of elements in a specified global row of the graph.
bool IndicesAreLocal() const
If column indices are in local range, this query returns true, otherwise it returns false.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false.
bool LowerTriangular() const
If graph is lower triangular in local index space, this query returns true, otherwise it returns fals...
bool IndicesAreGlobal() const
If column indices are in global range, this query returns true, otherwise it returns false.
bool UpperTriangular() const
If graph is upper triangular in local index space, this query returns true, otherwise it returns fals...
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int ReplaceDomainMapAndImporter(const Epetra_Map &NewDomainMap, const Epetra_Import *NewImporter)
Replaces the current DomainMap & Importer with the user-specified map object.
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.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
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 NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor's portion of the matrix.
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
const Epetra_Map & RangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
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.
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...
const Epetra_Map & DomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
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 IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements.
int NumPermuteIDs() const
Returns the number of elements that are local to the calling processor, but not part of the first Num...
int * ExportPIDs() const
List of processors to which elements will be sent, ExportLIDs() [i] will be sent to processor ExportP...
int * PermuteFromLIDs() const
List of elements in the source map that are permuted.
int NumExportIDs() const
Returns the number of elements that must be sent by the calling processor to other processors.
int NumSameIDs() const
Returns the number of elements that are identical between the source and target maps,...
int * RemoteLIDs() const
List of elements in the target map that are coming from other processors.
int NumRemoteIDs() const
Returns the number of elements that are not on the calling processor.
int * ExportLIDs() const
List of elements that will be sent to other processors.
const Epetra_BlockMap & SourceMap() const
Returns the SourceMap used to construct this importer.
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this importer.
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer.
int MaxValue()
Find maximum value.
int PutValue(int Value)
Set all elements of the vector to Value.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int MinValue()
Find minimum value.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
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.
int SetSeed(unsigned int Seed_in)
Set seed for Random function.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
Epetra_SerialComm: The Epetra Serial Communication Class.
static int GetPids(const Epetra_Import &Importer, std::vector< int > &pids, bool use_minus_one_for_local)
Epetra_Util GetPids 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 special_submap_import_test(Epetra_Comm &Comm)
int combine_mode_test(Epetra_Comm &Comm)
int test_import_gid(const char *name, Epetra_IntVector &Source, Epetra_IntVector &Target, const Epetra_Import &Import)
int alternate_import_constructor_test(Epetra_Comm &Comm)