Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
lesson05_redistribution.cpp
Go to the documentation of this file.
1
8// This defines useful macros like HAVE_MPI, which is defined if and
9// only if Epetra was built with MPI enabled.
10#include <Epetra_config.h>
11
12#ifdef HAVE_MPI
13# include <mpi.h>
14// Epetra's wrapper for MPI_Comm. This header file only exists if
15// Epetra was built with MPI enabled.
16# include <Epetra_MpiComm.h>
17#else
18# include <Epetra_SerialComm.h>
19#endif // HAVE_MPI
20
21#include <Epetra_CrsMatrix.h>
22#include <Epetra_Export.h>
23#include <Epetra_Map.h>
24#include <Epetra_Vector.h>
25#include <Epetra_Version.h>
26
27#include <sstream>
28#include <stdexcept>
29
30
31// The type of global indices. You could just set this to int,
32// but we want the example to work for Epetra64 as well.
33#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
34// Epetra was compiled only with 64-bit global index support,
35// so use 64-bit global indices.
36# define EXAMPLE_USES_64BIT_GLOBAL_INDICES 1
37typedef long long global_ordinal_type;
38#else
39# ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
40// Epetra was compiled only with 32-bit global index support,
41// so use 32-bit global indices.
42typedef int global_ordinal_type;
43# else
44# define EXAMPLE_USES_64BIT_GLOBAL_INDICES 1
45// Epetra was compiled with both 64-bit and 32-bit global index
46// support. Use 64-bit global indices for maximum generality.
47typedef long long global_ordinal_type;
48# endif // EPETRA_NO_64BIT_GLOBAL_INDICES
49#endif // EPETRA_NO_32BIT_GLOBAL_INDICES
50
51
52// Create and return a pointer to an example CrsMatrix, with row
53// distribution over the given Map. The caller is responsible for
54// freeing the result.
57{
58 const Epetra_Comm& comm = map.Comm ();
59
60 // Create an Epetra_CrsMatrix using the Map, with dynamic allocation.
61 Epetra_CrsMatrix* A = new Epetra_CrsMatrix (Copy, map, 3);
62
63 // The list of global indices owned by this MPI process.
64
65 const global_ordinal_type* myGblElts = NULL;
66 global_ordinal_type numGblElts = 0;
67#ifdef EXAMPLE_USES_64BIT_GLOBAL_INDICES
68 myGblElts = map.MyGlobalElements64 ();
69 numGblElts = map.NumGlobalElements64 ();
70#else
71 myGblElts = map.MyGlobalElements ();
72 numGblElts = map.NumGlobalElements ();
73#endif // EXAMPLE_USES_64BIT_GLOBAL_INDICES
74
75 // The number of global indices owned by this MPI process.
76 const int numMyElts = map.NumMyElements ();
77
78 // In general, tests like this really should synchronize across all
79 // processes. However, the likely cause for this case is a
80 // misconfiguration of Epetra, so we expect it to happen on all
81 // processes, if it happens at all.
82 if (numMyElts > 0 && myGblElts == NULL) {
83 throw std::logic_error ("Failed to get the list of global indices");
84 }
85
86 // Local error code for use below.
87 int lclerr = 0;
88
89 // Fill the sparse matrix, one row at a time.
90 double tempVals[3];
91 global_ordinal_type tempGblInds[3];
92 for (int i = 0; i < numMyElts; ++i) {
93 // A(0, 0:1) = [2, -1]
94 if (myGblElts[i] == 0) {
95 tempVals[0] = 2.0;
96 tempVals[1] = -1.0;
97 tempGblInds[0] = myGblElts[i];
98 tempGblInds[1] = myGblElts[i] + 1;
99 if (lclerr == 0) {
100 lclerr = A->InsertGlobalValues (myGblElts[i], 2, tempVals, tempGblInds);
101 }
102 if (lclerr != 0) {
103 break;
104 }
105 }
106 // A(N-1, N-2:N-1) = [-1, 2]
107 else if (myGblElts[i] == numGblElts - 1) {
108 tempVals[0] = -1.0;
109 tempVals[1] = 2.0;
110 tempGblInds[0] = myGblElts[i] - 1;
111 tempGblInds[1] = myGblElts[i];
112 if (lclerr == 0) {
113 lclerr = A->InsertGlobalValues (myGblElts[i], 2, tempVals, tempGblInds);
114 }
115 if (lclerr != 0) {
116 break;
117 }
118 }
119 // A(i, i-1:i+1) = [-1, 2, -1]
120 else {
121 tempVals[0] = -1.0;
122 tempVals[1] = 2.0;
123 tempVals[2] = -1.0;
124 tempGblInds[0] = myGblElts[i] - 1;
125 tempGblInds[1] = myGblElts[i];
126 tempGblInds[2] = myGblElts[i] + 1;
127 if (lclerr == 0) {
128 lclerr = A->InsertGlobalValues (myGblElts[i], 3, tempVals, tempGblInds);
129 }
130 if (lclerr != 0) {
131 break;
132 }
133 }
134 }
135
136 // If any process failed to insert at least one entry, throw.
137 int gblerr = 0;
138 (void) comm.MaxAll (&lclerr, &gblerr, 1);
139 if (gblerr != 0) {
140 if (A != NULL) {
141 delete A;
142 }
143 throw std::runtime_error ("Some process failed to insert an entry.");
144 }
145
146 // Tell the sparse matrix that we are done adding entries to it.
147 gblerr = A->FillComplete ();
148 if (gblerr != 0) {
149 if (A != NULL) {
150 delete A;
151 }
152 std::ostringstream os;
153 os << "A->FillComplete() failed with error code " << gblerr << ".";
154 throw std::runtime_error (os.str ());
155 }
156
157 return A;
158}
159
160
161void
162example (const Epetra_Comm& comm)
163{
164 // The global number of rows in the matrix A to create. We scale
165 // this relative to the number of (MPI) processes, so that no matter
166 // how many MPI processes you run, every process will have 10 rows.
167 const global_ordinal_type numGblElts = 10 * comm.NumProc ();
168 // The global min global index in all the Maps here.
169 const global_ordinal_type indexBase = 0;
170
171 // Local error code for use below.
172 //
173 // In the ideal case, we would use this to emulate behavior like
174 // that of Haskell's Maybe in the context of MPI. That is, if one
175 // process experiences an error, we don't want to abort early and
176 // cause the other processes to deadlock on MPI communication
177 // operators. Rather, we want to chain along the local error state,
178 // until we reach a point where it's natural to pass along that
179 // state with other processes. For example, if one is doing an
180 // MPI_Allreduce anyway, it makes sense to pass along one more bit
181 // of information: whether the calling process is in a local error
182 // state. Epetra's interface doesn't let one chain the local error
183 // state in this way, so we use extra collectives below to propagate
184 // that state. The code below uses very conservative error checks;
185 // typical user code would not need to be so conservative and could
186 // therefore avoid all the all-reduces.
187 int lclerr = 0;
188
189 // Construct a Map that is global (not locally replicated), but puts
190 // all the equations on MPI Proc 0.
191 const int procZeroMapNumLclElts = (comm.MyPID () == 0) ?
192 numGblElts :
193 static_cast<global_ordinal_type> (0);
194 Epetra_Map procZeroMap (numGblElts, procZeroMapNumLclElts, indexBase, comm);
195
196 // Construct a Map that puts approximately the same number of
197 // equations on each processor.
198 Epetra_Map globalMap (numGblElts, indexBase, comm);
199
200 // Create a sparse matrix using procZeroMap.
201 Epetra_CrsMatrix* A = createCrsMatrix (procZeroMap);
202 if (A == NULL) {
203 lclerr = 1;
204 }
205
206 // Make sure that sparse matrix creation succeeded. Normally you
207 // don't have to check this; we are being extra conservative because
208 // this example is also a test. Even though the matrix's rows live
209 // entirely on Process 0, the matrix is nonnull on all processes in
210 // its Map's communicator.
211 int gblerr = 0;
212 (void) comm.MaxAll (&lclerr, &gblerr, 1);
213 if (gblerr != 0) {
214 throw std::runtime_error ("createCrsMatrix returned NULL on at least one "
215 "process.");
216 }
217
218 //
219 // We've created a sparse matrix whose rows live entirely on MPI
220 // Process 0. Now we want to distribute it over all the processes.
221 //
222
223 // Redistribute the matrix. Since both the source and target Maps
224 // are one-to-one, we could use either an Import or an Export. If
225 // only the source Map were one-to-one, we would have to use an
226 // Import; if only the target Map were one-to-one, we would have to
227 // use an Export. We do not allow redistribution using Import or
228 // Export if neither source nor target Map is one-to-one.
229
230 // Make an export object with procZeroMap as the source Map, and
231 // globalMap as the target Map. The Export type has the same
232 // template parameters as a Map. Note that Export does not depend
233 // on the Scalar template parameter of the objects it
234 // redistributes. You can reuse the same Export for different
235 // Tpetra object types, or for Tpetra objects of the same type but
236 // different Scalar template parameters (e.g., Scalar=float or
237 // Scalar=double).
238 Epetra_Export exporter (procZeroMap, globalMap);
239
240 // Make a new sparse matrix whose row map is the global Map.
241 Epetra_CrsMatrix B (Copy, globalMap, 0);
242
243 // Redistribute the data, NOT in place, from matrix A (which lives
244 // entirely on Proc 0) to matrix B (which is distributed evenly over
245 // the processes).
246 //
247 // Export() has collective semantics, so we must always call it on
248 // all processes collectively. This is why we don't select on
249 // lclerr, as we do for the local operations above.
250 lclerr = B.Export (*A, exporter, Insert);
251
252 // Make sure that the Export succeeded. Normally you don't have to
253 // check this; we are being extra conservative because this example
254 // example is also a test. We test both min and max, since lclerr
255 // may be negative, zero, or positive.
256 gblerr = 0;
257 (void) comm.MinAll (&lclerr, &gblerr, 1);
258 if (gblerr != 0) {
259 throw std::runtime_error ("Export() failed on at least one process.");
260 }
261 (void) comm.MaxAll (&lclerr, &gblerr, 1);
262 if (gblerr != 0) {
263 throw std::runtime_error ("Export() failed on at least one process.");
264 }
265
266 // FillComplete has collective semantics, so we must always call it
267 // on all processes collectively. This is why we don't select on
268 // lclerr, as we do for the local operations above.
269 lclerr = B.FillComplete ();
270
271 // Make sure that FillComplete succeeded. Normally you don't have
272 // to check this; we are being extra conservative because this
273 // example is also a test. We test both min and max, since lclerr
274 // may be negative, zero, or positive.
275 gblerr = 0;
276 (void) comm.MinAll (&lclerr, &gblerr, 1);
277 if (gblerr != 0) {
278 throw std::runtime_error ("B.FillComplete() failed on at least one process.");
279 }
280 (void) comm.MaxAll (&lclerr, &gblerr, 1);
281 if (gblerr != 0) {
282 throw std::runtime_error ("B.FillComplete() failed on at least one process.");
283 }
284
285 if (A != NULL) {
286 delete A;
287 }
288}
289
290
291int
292main (int argc, char *argv[])
293{
294 using std::cout;
295 using std::endl;
296
297#ifdef HAVE_MPI
298 MPI_Init (&argc, &argv);
299 Epetra_MpiComm comm (MPI_COMM_WORLD);
300#else
302#endif // HAVE_MPI
303
304 const int myRank = comm.MyPID ();
305 const int numProcs = comm.NumProc ();
306
307 if (myRank == 0) {
308 // Print out the Epetra software version.
309 cout << Epetra_Version () << endl << endl
310 << "Total number of processes: " << numProcs << endl;
311 }
312
313 example (comm); // Run the whole example.
314
315 // This tells the Trilinos test framework that the test passed.
316 if (myRank == 0) {
317 cout << "End Result: TEST PASSED" << endl;
318 }
319
320#ifdef HAVE_MPI
321 (void) MPI_Finalize ();
322#endif // HAVE_MPI
323
324 return 0;
325}
@ Insert
@ Copy
std::string Epetra_Version()
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
long long * MyGlobalElements64() const
long long NumGlobalElements64() const
int NumGlobalElements() const
Number of elements across all processors.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
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 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 MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
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.
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.
Definition: Epetra_Export.h:62
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_SerialComm: The Epetra Serial Communication Class.
int main(int argc, char *argv[])
void example(const Epetra_Comm &comm)
long long global_ordinal_type
Epetra_CrsMatrix * createCrsMatrix(const Epetra_Map &map)