70 int main(
int argc,
char *argv[])
78 MPI_Init(&argc,&argv);
88 int MyPID = Comm.
MyPID();
90 bool verbose = (MyPID==0);
95 std::cout << Comm << std::endl;
101 std::cout <<
"Usage: " << argv[0] <<
" number_of_equations" << std::endl;
104 long long NumGlobalElements = std::atoi(argv[1]);
106 if (NumGlobalElements < NumProc)
109 std::cout <<
"numGlobalBlocks = " << NumGlobalElements
110 <<
" cannot be < number of processors = " << NumProc << std::endl;
123 std::vector<long long> MyGlobalElements(NumMyElements);
130 std::vector<int> NumNz(NumMyElements);
135 for (i=0; i<NumMyElements; i++)
136 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
150 std::vector<double> Values(2);
151 Values[0] = -1.0; Values[1] = -1.0;
152 std::vector<long long> Indices(2);
156 for (i=0; i<NumMyElements; i++)
158 if (MyGlobalElements[i]==0)
163 else if (MyGlobalElements[i] == NumGlobalElements-1)
165 Indices[0] = NumGlobalElements-2;
170 Indices[0] = MyGlobalElements[i]-1;
171 Indices[1] = MyGlobalElements[i]+1;
174 ierr = A.
InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
190 int niters = (int) NumGlobalElements*10;
191 double tolerance = 1.0e-2;
197 ierr +=
power_method(A, lambda, niters, tolerance, verbose);
199 double total_flops =counter.
Flops();
200 double MFLOPs = total_flops/elapsed_time/1000000.0;
203 std::cout <<
"\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
207 std::cout <<
"\nIncreasing magnitude of first diagonal term, solving again\n\n"
212 std::vector<double> Rowvals(numvals);
213 std::vector<long long> Rowinds(numvals);
215 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
224 ierr +=
power_method(A, lambda, niters, tolerance, verbose);
226 total_flops = counter.
Flops();
227 MFLOPs = total_flops/elapsed_time/1000000.0;
230 std::cout <<
"\n\nTotal MFLOPs for second solve = " << MFLOPs << std::endl<< std::endl;
244 double tolerance,
bool verbose) {
254 resid.SetFlopCounter(A);
261 double normz, residual;
265 for (
int iter = 0; iter < niters; iter++)
268 q.Scale(1.0/normz, z);
271 if (iter%100==0 || iter+1==niters)
273 resid.Update(1.0, z, -lambda, q, 0.0);
274 resid.Norm2(&residual);
275 if (verbose) std::cout <<
"Iter = " << iter <<
" Lambda = " << lambda
276 <<
" Residual of A*q - lambda*q = "
277 << residual << std::endl;
279 if (residual < tolerance) {
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
Epetra_Map: A class for partitioning vectors and matrices.
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
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 power_method(Epetra_CrsMatrix &A, double &lambda, int niters, double tolerance, bool verbose)
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int MyPID() const
Return my process ID.
Epetra_MpiComm: The Epetra MPI Communication Class.
std::string Epetra_Version()
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
Epetra_Time: The Epetra Timing Class.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
int NumMyElements() const
Number of elements on the calling processor.
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 NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
Epetra_SerialComm: The Epetra Serial Communication Class.
void ResetFlops()
Resets the number of floating point operations to zero for this multi-vector.
Epetra_Flops: The Epetra Floating Point Operations Class.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int main(int argc, char *argv[])
void ResetStartTime(void)
Epetra_Time function to reset the start time for a timer object.
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...
Epetra_Flops * GetFlopCounter() const
Get the pointer to the Epetra_Flops() object associated with this object, returns 0 if none...