Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
lesson03_power_method.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_Map.h>
23#include <Epetra_Vector.h>
24#include <Epetra_Version.h>
25
26#include <sstream>
27#include <stdexcept>
28
29// Implementation of the power method for estimating the eigenvalue of
30// maximum magnitude of a matrix. This function returns the
31// eigenvalue estimate.
32//
33// A: The sparse matrix or operator, as an Epetra_Operator.
34// niters: Maximum number of iterations of the power method.
35// tolerance: If the 2-norm of the residual A*x-lambda*x (for the
36// current eigenvalue estimate lambda) is less than this, stop
37// iterating.
38// out: output stream to which to print the current status of the
39// power method.
40double
42 const int niters,
43 const double tolerance);
44
45int
46main (int argc, char *argv[])
47{
48 using std::cout;
49 using std::endl;
50
51#ifdef HAVE_MPI
52 MPI_Init (&argc, &argv);
53 Epetra_MpiComm comm (MPI_COMM_WORLD);
54#else
56#endif // HAVE_MPI
57
58 const int myRank = comm.MyPID ();
59 const int numProcs = comm.NumProc ();
60
61 if (myRank == 0) {
62 // Print out the Epetra software version.
63 cout << Epetra_Version () << endl << endl
64 << "Total number of processes: " << numProcs << endl;
65 }
66
67 // The type of global indices. You could just set this to int,
68 // but we want the example to work for Epetra64 as well.
69#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
70 // Epetra was compiled only with 64-bit global index support, so use
71 // 64-bit global indices.
72 typedef long long global_ordinal_type;
73#else
74 // Epetra was compiled with 32-bit global index support. If
75 // EPETRA_NO_64BIT_GLOBAL_INDICES is defined, it does not also
76 // support 64-bit indices.
77 typedef int global_ordinal_type;
78#endif // EPETRA_NO_32BIT_GLOBAL_INDICES
79
80 // The number of rows and columns in the matrix.
81 const global_ordinal_type numGlobalElements = 50;
82
83 // Construct a Map that puts approximately the same number of
84 // equations on each processor.
85 const global_ordinal_type indexBase = 0;
86 Epetra_Map map (numGlobalElements, indexBase, comm);
87
88 // Get the list of global indices that this process owns. In this
89 // example, this is unnecessary, because we know that we created a
90 // contiguous Map (see above). (Thus, we really only need the min
91 // and max global index on this process.) However, in general, we
92 // don't know what global indices the Map owns, so if we plan to add
93 // entries into the sparse matrix using global indices, we have to
94 // get the list of global indices this process owns.
95 const int numMyElements = map.NumMyElements ();
96
97 global_ordinal_type* myGlobalElements = NULL;
98
99#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
100 myGlobalElements = map.MyGlobalElements64 ();
101#else
102 myGlobalElements = map.MyGlobalElements ();
103#endif // EPETRA_NO_32BIT_GLOBAL_INDICES
104
105 // In general, tests like this really should synchronize across all
106 // processes. However, the likely cause for this case is a
107 // misconfiguration of Epetra, so we expect it to happen on all
108 // processes, if it happens at all.
109 if (numMyElements > 0 && myGlobalElements == NULL) {
110 throw std::logic_error ("Failed to get the list of global indices");
111 }
112
113 if (myRank == 0) {
114 cout << endl << "Creating the sparse matrix" << endl;
115 }
116
117 // Create a Epetra sparse matrix whose rows have distribution given
118 // by the Map. The max number of entries per row is 3.
119 Epetra_CrsMatrix A (Copy, map, 3);
120
121 // Local error code for use below.
122 int lclerr = 0;
123
124 // Fill the sparse matrix, one row at a time. InsertGlobalValues
125 // adds entries to the sparse matrix, using global column indices.
126 // It changes both the graph structure and the values.
127 double tempVals[3];
128 global_ordinal_type tempGblInds[3];
129 for (int i = 0; i < numMyElements; ++i) {
130 // A(0, 0:1) = [2, -1]
131 if (myGlobalElements[i] == 0) {
132 tempVals[0] = 2.0;
133 tempVals[1] = -1.0;
134 tempGblInds[0] = myGlobalElements[i];
135 tempGblInds[1] = myGlobalElements[i] + 1;
136 if (lclerr == 0) {
137 lclerr = A.InsertGlobalValues (myGlobalElements[i], 2, tempVals, tempGblInds);
138 }
139 if (lclerr != 0) {
140 break;
141 }
142 }
143 // A(N-1, N-2:N-1) = [-1, 2]
144 else if (myGlobalElements[i] == numGlobalElements - 1) {
145 tempVals[0] = -1.0;
146 tempVals[1] = 2.0;
147 tempGblInds[0] = myGlobalElements[i] - 1;
148 tempGblInds[1] = myGlobalElements[i];
149 if (lclerr == 0) {
150 lclerr = A.InsertGlobalValues (myGlobalElements[i], 2, tempVals, tempGblInds);
151 }
152 if (lclerr != 0) {
153 break;
154 }
155 }
156 // A(i, i-1:i+1) = [-1, 2, -1]
157 else {
158 tempVals[0] = -1.0;
159 tempVals[1] = 2.0;
160 tempVals[2] = -1.0;
161 tempGblInds[0] = myGlobalElements[i] - 1;
162 tempGblInds[1] = myGlobalElements[i];
163 tempGblInds[2] = myGlobalElements[i] + 1;
164 if (lclerr == 0) {
165 lclerr = A.InsertGlobalValues (myGlobalElements[i], 3, tempVals, tempGblInds);
166 }
167 if (lclerr != 0) {
168 break;
169 }
170 }
171 }
172
173 // If any process failed to insert at least one entry, throw.
174 int gblerr = 0;
175 (void) comm.MaxAll (&lclerr, &gblerr, 1);
176 if (gblerr != 0) {
177 throw std::runtime_error ("Some process failed to insert an entry.");
178 }
179
180 // Tell the sparse matrix that we are done adding entries to it.
181 gblerr = A.FillComplete ();
182 if (gblerr != 0) {
183 std::ostringstream os;
184 os << "A.FillComplete() failed with error code " << gblerr << ".";
185 throw std::runtime_error (os.str ());
186 }
187
188 // Number of iterations
189 const int niters = 500;
190 // Desired (absolute) residual tolerance
191 const double tolerance = 1.0e-2;
192
193 // Run the power method and report the result.
194 double lambda = powerMethod (A, niters, tolerance);
195 if (myRank == 0) {
196 cout << endl << "Estimated max eigenvalue: " << lambda << endl;
197 }
198
199 //
200 // Now we're going to change values in the sparse matrix and run the
201 // power method again.
202 //
203
204 //
205 // Increase diagonal dominance
206 //
207 if (myRank == 0) {
208 cout << endl << "Increasing magnitude of A(0,0), solving again" << endl;
209 }
210
211 if (A.RowMap ().MyGID (0)) {
212 // Get a copy of the row with with global index 0. Modify the
213 // diagonal entry of that row. Submit the modified values to the
214 // matrix.
215 const global_ordinal_type gidOfFirstRow = 0;
216 // Since 0 is a GID in the row Map on the calling process,
217 // lidOfFirstRow is a valid LID of that GID in the row Map.
218 const int lidOfFirstRow = A.RowMap ().LID (gidOfFirstRow);
219 int numEntriesInRow = A.NumMyEntries (lidOfFirstRow);
220 double* rowvals = new double [numEntriesInRow];
221 global_ordinal_type* rowinds = new global_ordinal_type [numEntriesInRow];
222
223 // Get a copy of the entries and column indices of global row 0.
224 // Get global column indices, so that we can figure out which
225 // entry corresponds to the diagonal entry in this row. (The row
226 // Map and column Map of the matrix may differ, so their local
227 // indices may not be the same.)
228 //
229 // Note that it's legal (though we don't exercise it in this
230 // example) for the row Map of the sparse matrix not to be one to
231 // one. This means that more than one process might own entries
232 // in the first row. In general, multiple processes might own the
233 // (0,0) entry, so that the global A(0,0) value is really the sum
234 // of all processes' values for that entry. However, scaling the
235 // entry by a constant factor distributes across that sum, so it's
236 // OK to do so.
237 if (lclerr == 0) {
238 lclerr = A.ExtractGlobalRowCopy (gidOfFirstRow,
239 numEntriesInRow, numEntriesInRow,
240 rowvals, rowinds);
241 }
242 if (lclerr == 0) { // no error
243 for (int i = 0; i < numEntriesInRow; ++i) {
244 if (rowinds[i] == gidOfFirstRow) {
245 // We have found the diagonal entry; modify it.
246 rowvals[i] *= 10.0;
247 }
248 }
249 // "Replace global values" means modify the values, but not the
250 // structure of the sparse matrix. If the specified columns
251 // aren't already populated in this row on this process, then this
252 // method fails and returns nonzero. Since we have already called
253 // FillComplete() on this matrix, we may not modify its graph
254 // structure any more.
255 if (lclerr == 0) {
256 lclerr = A.ReplaceGlobalValues (gidOfFirstRow, numEntriesInRow,
257 rowvals, rowinds);
258 }
259 }
260
261 if (rowvals != NULL) {
262 delete [] rowvals;
263 }
264 if (rowinds != NULL) {
265 delete [] rowinds;
266 }
267 }
268
269 // If the owning process(es) of global row 0 failed to replace the
270 // one entry, throw.
271 gblerr = 0;
272 (void) comm.MaxAll (&lclerr, &gblerr, 1);
273 if (gblerr != 0) {
274 throw std::runtime_error ("One of the owning process(es) of global "
275 "row 0 failed to replace an entry.");
276 }
277
278 // Run the power method again.
279 lambda = powerMethod (A, niters, tolerance);
280 if (myRank == 0) {
281 cout << endl << "Estimated max eigenvalue: " << lambda << endl;
282 }
283
284 // This tells the Trilinos test framework that the test passed.
285 if (myRank == 0) {
286 cout << "End Result: TEST PASSED" << endl;
287 }
288
289#ifdef HAVE_MPI
290 (void) MPI_Finalize ();
291#endif // HAVE_MPI
292
293 return 0;
294}
295
296
297double
299 const int niters,
300 const double tolerance)
301{
302 using std::cout;
303 using std::endl;
304
305 // An Operator doesn't have a Comm, but its domain Map does.
306 const Epetra_Comm& comm = A.OperatorDomainMap ().Comm ();
307 const int myRank = comm.MyPID ();
308
309 // Create three vectors for iterating the power method. Since the
310 // power method computes z = A*q, q should be in the domain of A and
311 // z should be in the range. (Obviously the power method requires
312 // that the domain and the range are equal, but it's a good idea to
313 // get into the habit of thinking whether a particular vector
314 // "belongs" in the domain or range of the matrix.) The residual
315 // vector "resid" is of course in the range of A.
318 Epetra_Vector resid (A.OperatorRangeMap ());
319
320 // Local error code for use below.
321 int lclerr = 0;
322
323 // Fill the iteration vector z with random numbers to start. Don't
324 // have grand expectations about the quality of our pseudorandom
325 // number generator; this is usually good enough for eigensolvers.
326 //
327 // If this call fails, just let the code below finish before trying
328 // to catch the error globally. Ditto for other calls below.
329 lclerr = z.Random ();
330
331 // lambda: the current approximation of the eigenvalue of maximum magnitude.
332 // normz: the 2-norm of the current iteration vector z.
333 // residual: the 2-norm of the current residual vector "resid"
334 double lambda = 0.0;
335 double normz = 0.0;
336 double residual = 0.0;
337
338 const double zero = 0.0;
339 const double one = 1.0;
340
341 // How often to report progress in the power method. Reporting
342 // progress requires computing a residual which can be expensive.
343 // However, if you don't compute the residual often enough, you
344 // might keep iterating even after you've converged.
345 const int reportFrequency = 10;
346
347 // Do the power method, until the method has converged or the
348 // maximum iteration count has been reached.
349 for (int iter = 0; iter < niters; ++iter) {
350 // If you feel confident that your code is correct, you may omit
351 // everything having to do with lclerr, and just write the following:
352 //
353 // z.Norm2 (&normz); // Compute the 2-norm of z
354 // q.Scale (one / normz, z); // q := z / normz
355 // A.Apply (q, z); // z := A * q
356 // q.Dot (z, &lambda); // Approx. max eigenvalue
357
358 lclerr = (lclerr == 0) ? z.Norm2 (&normz) : lclerr;
359 lclerr = (lclerr == 0) ? q.Scale (one / normz, z) : lclerr;
360 lclerr = (lclerr == 0) ? A.Apply (q, z) : lclerr;
361 lclerr = (lclerr == 0) ? q.Dot (z, &lambda) : lclerr;
362
363 // Compute and report the residual norm every reportFrequency
364 // iterations, or if we've reached the maximum iteration count.
365 if (iter % reportFrequency == 0 || iter + 1 == niters) {
366 // If you feel confident that your code is correct, you may omit
367 // everything having to do with lclerr, and just write the
368 // following:
369 //
370 // resid.Update (one, z, -lambda, q, zero); // z := A*q - lambda*q
371 // resid.Norm2 (&residual); // 2-norm of the residual vector
372
373 lclerr = (lclerr == 0) ? resid.Update (one, z, -lambda, q, zero) : lclerr;
374 lclerr = (lclerr == 0) ? resid.Norm2 (&residual) : lclerr;
375
376 if (myRank == 0) {
377 cout << "Iteration " << iter << ":" << endl
378 << "- lambda = " << lambda << endl
379 << "- ||A*q - lambda*q||_2 = " << residual << endl;
380 }
381 }
382 if (residual < tolerance) {
383 if (myRank == 0) {
384 cout << "Converged after " << iter << " iterations" << endl;
385 }
386 break;
387 } else if (iter + 1 == niters) {
388 if (myRank == 0) {
389 cout << "Failed to converge after " << niters << " iterations" << endl;
390 }
391 break;
392 }
393 }
394
395 // If any process failed to insert at least one entry, throw.
396 int gblerr = 0;
397 (void) comm.MaxAll (&lclerr, &gblerr, 1);
398 if (gblerr != 0) {
399 throw std::runtime_error ("The power method failed in some way.");
400 }
401
402 return lambda;
403}
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
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
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...
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 MyPID() const =0
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
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 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 NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this 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.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_MpiComm: The Epetra MPI Communication Class.
int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const
Epetra_MpiComm Global Max 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.
Epetra_Operator: A pure virtual class for using real-valued double-precision operators.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
virtual const Epetra_Map & OperatorDomainMap() const =0
Returns the Epetra_Map object associated with the domain of this operator.
virtual const Epetra_Map & OperatorRangeMap() const =0
Returns the Epetra_Map object associated with the range of this operator.
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int main(int argc, char *argv[])
double powerMethod(const Epetra_Operator &A, const int niters, const double tolerance)
long long global_ordinal_type