10#include <Epetra_config.h>
43 const double tolerance);
46main (
int argc,
char *argv[])
52 MPI_Init (&argc, &argv);
58 const int myRank = comm.
MyPID ();
59 const int numProcs = comm.
NumProc ();
64 <<
"Total number of processes: " << numProcs << endl;
69#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
86 Epetra_Map map (numGlobalElements, indexBase, comm);
99#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
109 if (numMyElements > 0 && myGlobalElements == NULL) {
110 throw std::logic_error (
"Failed to get the list of global indices");
114 cout << endl <<
"Creating the sparse matrix" << endl;
129 for (
int i = 0; i < numMyElements; ++i) {
131 if (myGlobalElements[i] == 0) {
134 tempGblInds[0] = myGlobalElements[i];
135 tempGblInds[1] = myGlobalElements[i] + 1;
144 else if (myGlobalElements[i] == numGlobalElements - 1) {
147 tempGblInds[0] = myGlobalElements[i] - 1;
148 tempGblInds[1] = myGlobalElements[i];
161 tempGblInds[0] = myGlobalElements[i] - 1;
162 tempGblInds[1] = myGlobalElements[i];
163 tempGblInds[2] = myGlobalElements[i] + 1;
175 (void) comm.
MaxAll (&lclerr, &gblerr, 1);
177 throw std::runtime_error (
"Some process failed to insert an entry.");
183 std::ostringstream os;
184 os <<
"A.FillComplete() failed with error code " << gblerr <<
".";
185 throw std::runtime_error (os.str ());
189 const int niters = 500;
191 const double tolerance = 1.0e-2;
194 double lambda =
powerMethod (A, niters, tolerance);
196 cout << endl <<
"Estimated max eigenvalue: " << lambda << endl;
208 cout << endl <<
"Increasing magnitude of A(0,0), solving again" << endl;
218 const int lidOfFirstRow = A.
RowMap ().
LID (gidOfFirstRow);
220 double* rowvals =
new double [numEntriesInRow];
239 numEntriesInRow, numEntriesInRow,
243 for (
int i = 0; i < numEntriesInRow; ++i) {
244 if (rowinds[i] == gidOfFirstRow) {
261 if (rowvals != NULL) {
264 if (rowinds != NULL) {
272 (void) comm.
MaxAll (&lclerr, &gblerr, 1);
274 throw std::runtime_error (
"One of the owning process(es) of global "
275 "row 0 failed to replace an entry.");
281 cout << endl <<
"Estimated max eigenvalue: " << lambda << endl;
286 cout <<
"End Result: TEST PASSED" << endl;
290 (void) MPI_Finalize ();
300 const double tolerance)
307 const int myRank = comm.
MyPID ();
336 double residual = 0.0;
338 const double zero = 0.0;
339 const double one = 1.0;
345 const int reportFrequency = 10;
349 for (
int iter = 0; iter < niters; ++iter) {
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;
365 if (iter % reportFrequency == 0 || iter + 1 == niters) {
373 lclerr = (lclerr == 0) ? resid.
Update (one, z, -lambda, q, zero) : lclerr;
374 lclerr = (lclerr == 0) ? resid.
Norm2 (&residual) : lclerr;
377 cout <<
"Iteration " << iter <<
":" << endl
378 <<
"- lambda = " << lambda << endl
379 <<
"- ||A*q - lambda*q||_2 = " << residual << endl;
382 if (residual < tolerance) {
384 cout <<
"Converged after " << iter <<
" iterations" << endl;
387 }
else if (iter + 1 == niters) {
389 cout <<
"Failed to converge after " << niters <<
" iterations" << endl;
397 (void) comm.
MaxAll (&lclerr, &gblerr, 1);
399 throw std::runtime_error (
"The power method failed in some way.");
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.
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.
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