Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
B.cc
Go to the documentation of this file.
1/*
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42*/
43
44#include <stdio.h>
45#include <stdlib.h>
46#include <ctype.h>
47#include <assert.h>
48#include <string.h>
49#include <math.h>
50#ifdef PETRA_MPI
51#include "mpi.h"
52#endif
53#ifndef __cplusplus
54#define __cplusplus
55#endif
56#include "Petra_Comm.h"
57#include "Petra_Map.h"
58#include "Petra_Time.h"
59#include "Petra_RDP_CRS_Matrix.h"
60
61int main(int argc, char *argv[])
62{
63 int ierr = 0, i, j;
64 bool debug = false;
65
66#ifdef PETRA_MPI
67
68 // Initialize MPI
69
70 MPI_Init(&argc,&argv);
71 int size, rank; // Number of MPI processes, My process ID
72
73 MPI_Comm_size(MPI_COMM_WORLD, &size);
74 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
75
76#else
77
78 int size = 1; // Serial case (not using MPI)
79 int rank = 0;
80
81#endif
82
83 bool verbose = false;
84
85 // Check if we should print results to standard out
86 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
87
88
89
90#ifdef PETRA_MPI
91 Petra_Comm & Comm = *new Petra_Comm( MPI_COMM_WORLD );
92#else
93 Petra_Comm & Comm = *new Petra_Comm();
94#endif
95
96
97 //char tmp;
98 //if (rank==0) cout << "Press any key to continue..."<< endl;
99 //if (rank==0) cin >> tmp;
100 //Comm.Barrier();
101
102 int MyPID = Comm.MyPID();
103 int NumProc = Comm.NumProc();
104 if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
105 << " is alive."<<endl;
106
107 bool verbose1 = verbose;
108
109 // Redefine verbose to only print on PE 0
110 if (verbose && rank!=0) verbose = false;
111
112 int NumMyPoints = 10000;
113 int NumGlobalPoints = NumMyPoints*NumProc+minfn(NumProc,3);
114 if (MyPID < 3) NumMyPoints++;
115 int IndexBase = 0;
116 bool DistributedGlobal = (NumGlobalPoints>NumMyPoints);
117
118 // Construct a Source Map that puts approximately the same Number of equations on each processor in
119 // uniform global ordering
120
121 Petra_Map& SourceMap = *new Petra_Map(NumGlobalPoints, NumMyPoints, 0, Comm);
122
123 // Get update list and number of local equations from newly created Map
124 int NumMyElements = SourceMap.NumMyElements();
125 int * SourceMyGlobalElements = new int[NumMyElements];
126 SourceMap.MyGlobalElements(SourceMyGlobalElements);
127
128
129 // Construct a Target Map that will contain:
130 // some unchanged elements (relative to the soure map),
131 // some permuted elements
132 // some off-processor elements
133 Petra_RDP_Vector & RandVec = *new Petra_RDP_Vector(SourceMap);
134 RandVec.Random(); // This creates a vector of random numbers between negative one and one.
135
136 int *TargetMyGlobalElements = new int[NumMyElements];
137
138 for (i=0; i< NumMyPoints/2; i++) TargetMyGlobalElements[i] = i; // Half will be the same...
139 for (i=NumMyPoints/2; i<NumMyPoints; i++) {
140 int index = abs((int)(((double) (NumGlobalPoints-1) ) * RandVec[i]));
141 TargetMyGlobalElements[i] = minfn(NumGlobalPoints-1,maxfn(0,index));
142 }
143
144 int NumSameIDs = 0;
145 int NumPermutedIDs = 0;
146 int NumRemoteIDs = 0;
147 bool StillContiguous = true;
148 for (i=0; i < NumMyPoints; i++) {
149 if (SourceMyGlobalElements[i]==TargetMyGlobalElements[i] && StillContiguous)
150 NumSameIDs++;
151 else if (SourceMap.MyGID(TargetMyGlobalElements[i])) {
152 StillContiguous = false;
153 NumPermutedIDs++;
154 }
155 else {
156 StillContiguous = false;
157 NumRemoteIDs++;
158 }
159 }
160 assert(NumMyPoints==NumSameIDs+NumPermutedIDs+NumRemoteIDs);
161
162 Petra_Map & TargetMap = *new Petra_Map(-1, NumMyElements, TargetMyGlobalElements, 0, Comm);
163
164 // Create a multivector whose elements are GlobalID * (column number +1)
165
166 int NumVectors = 3;
167 Petra_RDP_MultiVector & SourceMultiVector = *new Petra_RDP_MultiVector(SourceMap, NumVectors);
168 for (j=0; j < NumVectors; j++)
169 for (i=0; i < NumMyElements; i++)
170 SourceMultiVector[j][i] = (double) SourceMyGlobalElements[i]*(j+1);
171
172 // Create a target multivector that we will fill using an Import
173
174 Petra_RDP_MultiVector & TargetMultiVector = *new Petra_RDP_MultiVector(TargetMap, NumVectors);
175
176 Petra_Import & Importer = *new Petra_Import(TargetMap, SourceMap);
177
178 assert(TargetMultiVector.Import(SourceMultiVector, Importer, Insert)==0);
179
180 // Test Target against expected values
181
182 for (j=0; j < NumVectors; j++)
183 for (i=0; i < NumMyElements; i++) {
184 if (TargetMultiVector[j][i]!= (double) TargetMyGlobalElements[i]*(j+1))
185 cout << "TargetMultiVector["<<i<<"]["<<j<<"] = " << TargetMultiVector[j][i]
186 << " TargetMyGlobalElements[i]*(j+1) = " << TargetMyGlobalElements[i]*(j+1) << endl;
187 assert(TargetMultiVector[j][i]== (double) TargetMyGlobalElements[i]*(j+1));
188 }
189
190 if (verbose) cout << "MultiVector Import using Importer Check OK" << endl << endl;
191
192
194
195 // Now use Importer to do an export
196
197 Petra_RDP_Vector & TargetVector = *new Petra_RDP_Vector(SourceMap);
198 Petra_RDP_Vector & ExpectedTarget = *new Petra_RDP_Vector(SourceMap);
199 Petra_RDP_Vector & SourceVector = *new Petra_RDP_Vector(TargetMap);
200
201 NumSameIDs = Importer.NumSameIDs();
202 int NumPermuteIDs = Importer.NumPermuteIDs();
203 int NumExportIDs = Importer.NumExportIDs();
204 int *PermuteToLIDs = Importer.PermuteToLIDs();
205 int *PermuteFromLIDs = Importer.PermuteFromLIDs();
206 int *ExportLIDs = Importer.ExportLIDs();
207 int *ExportPIDs = Importer.ExportPIDs();
208
209 for (i=0; i < NumSameIDs; i++) ExpectedTarget[i] = (double) (MyPID+1);
210 for (i=0; i < NumPermuteIDs; i++) ExpectedTarget[PermuteFromLIDs[i]] =
211 (double) (MyPID+1);
212 for (i=0; i < NumExportIDs; i++) ExpectedTarget[ExportLIDs[i]] +=
213 (double) (ExportPIDs[i]+1);
214
215 for (i=0; i < NumMyElements; i++) SourceVector[i] = (double) (MyPID+1);
216
217 assert(TargetVector.Export(SourceVector, Importer, Add)==0);
218
219 for (i=0; i < NumMyElements; i++) {
220 if (TargetVector[i]!= ExpectedTarget[i])
221 cout << " TargetVector["<<i<<"] = " << TargetVector[i]
222 << " ExpectedTarget["<<i<<"] = " << ExpectedTarget[i] << " on PE " << MyPID << endl;
223 assert(TargetVector[i]== ExpectedTarget[i]);
224 }
225
226 if (verbose) cout << "Vector Export using Importer Check OK" << endl << endl;
227
228
229
231 // Build a tridiagonal system two ways:
232 // 1) From "standard" matrix view where equations are uniquely owned.
233 // 2) From 1D PDE view where nodes (equations) between processors are shared and partial contributions are done
234 // in parallel, then merged together at the end of the construction process.
235 //
237
238
239
240 // Construct a Standard Map that puts approximately the same number of equations on each processor in
241 // uniform global ordering
242
243 Petra_Map& StandardMap = *new Petra_Map(NumGlobalPoints, NumMyPoints, 0, Comm);
244
245 // Get update list and number of local equations from newly created Map
246 NumMyElements = StandardMap.NumMyElements();
247 int * StandardMyGlobalElements = new int[NumMyElements];
248 StandardMap.MyGlobalElements(StandardMyGlobalElements);
249
250
251 // Create a standard Petra_CRS_Graph
252
253 Petra_CRS_Graph& StandardGraph = *new Petra_CRS_Graph(Copy, StandardMap, 3);
254 assert(!StandardGraph.IndicesAreGlobal());
255 assert(!StandardGraph.IndicesAreLocal());
256
257 // Add rows one-at-a-time
258 // Need some vectors to help
259 // Off diagonal Values will always be -1
260
261
262 int *Indices = new int[2];
263 int NumEntries;
264
265 for (i=0; i<NumMyPoints; i++)
266 {
267 if (StandardMyGlobalElements[i]==0)
268 {
269 Indices[0] = 1;
270 NumEntries = 1;
271 }
272 else if (StandardMyGlobalElements[i] == NumGlobalPoints-1)
273 {
274 Indices[0] = NumGlobalPoints-2;
275 NumEntries = 1;
276 }
277 else
278 {
279 Indices[0] = StandardMyGlobalElements[i]-1;
280 Indices[1] = StandardMyGlobalElements[i]+1;
281 NumEntries = 2;
282 }
283 assert(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], NumEntries, Indices)==0);
284 assert(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], 1, StandardMyGlobalElements+i)==0); // Put in the diagonal entry
285 }
286
287 // Finish up
288 assert(StandardGraph.IndicesAreGlobal());
289 assert(StandardGraph.FillComplete()==0);
290 assert(StandardGraph.IndicesAreLocal());
291 assert(!StandardGraph.StorageOptimized());
292 StandardGraph.OptimizeStorage();
293 assert(StandardGraph.StorageOptimized());
294 assert(!StandardGraph.UpperTriangular());
295 assert(!StandardGraph.LowerTriangular());
296
297
298 // Create Petra_RDP_CRS_Matrix using the just-built graph
299
300 Petra_RDP_CRS_Matrix& StandardMatrix = *new Petra_RDP_CRS_Matrix(Copy, StandardGraph);
301 assert(!StandardMatrix.IndicesAreGlobal());
302 assert(StandardMatrix.IndicesAreLocal());
303
304 // Add rows one-at-a-time
305 // Need some vectors to help
306 // Off diagonal Values will always be -1
307
308
309 double *Values = new double[2];
310 Values[0] = -1.0; Values[1] = -1.0;
311 double two = 2.0;
312
313 for (i=0; i<NumMyPoints; i++)
314 {
315 if (StandardMyGlobalElements[i]==0)
316 {
317 Indices[0] = 1;
318 NumEntries = 1;
319 }
320 else if (StandardMyGlobalElements[i] == NumGlobalPoints-1)
321 {
322 Indices[0] = NumGlobalPoints-2;
323 NumEntries = 1;
324 }
325 else
326 {
327 Indices[0] = StandardMyGlobalElements[i]-1;
328 Indices[1] = StandardMyGlobalElements[i]+1;
329 NumEntries = 2;
330 }
331 assert(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], NumEntries, Values, Indices)==0);
332 assert(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], 1, &two, StandardMyGlobalElements+i)==0); // Put in the diagonal entry
333 }
334
335 // Finish up
336 assert(StandardMatrix.IndicesAreLocal());
337 assert(StandardMatrix.FillComplete()==0);
338 assert(StandardMatrix.IndicesAreLocal());
339 assert(StandardMatrix.StorageOptimized());
340 assert(!StandardMatrix.UpperTriangular());
341 assert(!StandardMatrix.LowerTriangular());
342
343 // Construct an Overlapped Map of StandardMap that include the endpoints from two neighboring processors.
344
345 int OverlapNumMyElements;
346 int OverlapMinMyGID;
347
348 OverlapNumMyElements = NumMyElements + 1;
349 if (MyPID==0) OverlapNumMyElements--;
350
351 if (MyPID==0) OverlapMinMyGID = StandardMap.MinMyGID();
352 else OverlapMinMyGID = StandardMap.MinMyGID()-1;
353
354 int * OverlapMyGlobalElements = new int[OverlapNumMyElements];
355
356 for (i=0; i< OverlapNumMyElements; i++) OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
357
358 Petra_Map& OverlapMap = *new Petra_Map(-1, OverlapNumMyElements, OverlapMyGlobalElements, 0, Comm);
359
360 // Create the Overlap Petra_Matrix
361
362 Petra_RDP_CRS_Matrix& OverlapMatrix = *new Petra_RDP_CRS_Matrix(Copy, OverlapMap, 4);
363 assert(!OverlapMatrix.IndicesAreGlobal());
364 assert(!OverlapMatrix.IndicesAreLocal());
365
366 // Add matrix element one cell at a time.
367 // Each cell does an incoming and outgoing flux calculation
368
369
370 double pos_one = 1.0;
371 double neg_one = -1.0;
372
373 for (i=0; i<OverlapNumMyElements; i++)
374 {
375 int node_left = OverlapMyGlobalElements[i]-1;
376 int node_center = node_left + 1;
377 int node_right = node_left + 2;
378 if (i>0) {
379 if (node_left>-1)
380 assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_left)==0);
381 assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
382 }
383 if (i<OverlapNumMyElements-1) {
384 assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
385 if (node_right<NumGlobalPoints)
386 assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_right)==0);
387 }
388 }
389
390 // Handle endpoints
391 if (MyPID==0) {
392 int node_center = 0;
393 assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
394 }
395 if (MyPID==NumProc-1) {
396 int node_center = OverlapMyGlobalElements[OverlapNumMyElements-1];
397 assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
398 }
399
400 assert(OverlapMatrix.FillComplete()==0);
401
402 // Make a gathered matrix from OverlapMatrix. It should be identical to StandardMatrix
403
404 Petra_RDP_CRS_Matrix& GatheredMatrix = *new Petra_RDP_CRS_Matrix(Copy, StandardGraph);
405 Petra_Export & Exporter = *new Petra_Export(OverlapMap, StandardMap);
406 assert(GatheredMatrix.Export(OverlapMatrix, Exporter, Add)==0);
407 assert(GatheredMatrix.FillComplete()==0);
408
409 // Check if entries of StandardMatrix and GatheredMatrix are identical
410
411 int StandardNumEntries, GatheredNumEntries;
412 int * StandardIndices, * GatheredIndices;
413 double * StandardValues, * GatheredValues;
414
415 int StandardNumMyNonzeros = StandardMatrix.NumMyNonzeros();
416 int GatheredNumMyNonzeros = GatheredMatrix.NumMyNonzeros();
417 assert(StandardNumMyNonzeros==GatheredNumMyNonzeros);
418
419 int StandardNumMyRows = StandardMatrix.NumMyRows();
420 int GatheredNumMyRows = GatheredMatrix.NumMyRows();
421 assert(StandardNumMyRows==GatheredNumMyRows);
422
423 for (i=0; i< StandardNumMyRows; i++)
424 {
425 assert(StandardMatrix.ExtractMyRowView(i, StandardNumEntries, StandardValues, StandardIndices)==0);
426 assert(GatheredMatrix.ExtractMyRowView(i, GatheredNumEntries, GatheredValues, GatheredIndices)==0);
427 assert(StandardNumEntries==GatheredNumEntries);
428 for (j=0; j < StandardNumEntries; j++) {
429 //if (StandardIndices[j]!=GatheredIndices[j])
430 // cout << "MyPID = " << MyPID << " i = " << i << " StandardIndices[" << j << "] = " << StandardIndices[j]
431 // << " GatheredIndices[" << j << "] = " << GatheredIndices[j] << endl;
432 //if (StandardValues[j]!=GatheredValues[j])
433 //cout << "MyPID = " << MyPID << " i = " << i << " StandardValues[" << j << "] = " << StandardValues[j]
434 // << " GatheredValues[" << j << "] = " << GatheredValues[j] << endl;
435 assert(StandardIndices[j]==GatheredIndices[j]);
436 assert(StandardValues[j]==GatheredValues[j]);
437 }
438 }
439
440 if (verbose) cout << "Matrix Export Check OK" << endl;
441 // Release all objects
442
443 delete &SourceVector;
444 delete &TargetVector;
445 delete &ExpectedTarget;
446
447
448 delete &Importer;
449 delete &SourceMap;
450 delete &TargetMap;
451
452 delete [] SourceMyGlobalElements;
453 delete [] TargetMyGlobalElements;
454
455 delete &SourceMultiVector;
456 delete &TargetMultiVector;
457 delete &RandVec;
458
459 delete &Exporter;
460 delete &GatheredMatrix;
461 delete &OverlapMatrix;
462 delete &OverlapMap;
463 delete [] OverlapMyGlobalElements;
464
465 delete &StandardMatrix;
466 delete &StandardGraph;
467 delete &StandardMap;
468 delete [] StandardMyGlobalElements;
469
470 delete [] Values;
471 delete [] Indices;
472 delete &Comm;
473
474#ifdef PETRA_MPI
475 MPI_Finalize() ;
476#endif
477
478/* end main
479*/
480return ierr ;
481}
int main(int argc, char *argv[])
Definition: B.cc:61
@ Insert
@ Copy