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) {
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
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.
Epetra_Flops * GetFlopCounter() const
Get the pointer to the Epetra_Flops() object associated with this object, returns 0 if none...
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.
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
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.
double Flops() const
Returns the number of floating point operations with this object and resets the count.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
int NumMyElements() const
Number of elements on the calling processor.
int power_method(Epetra_CrsMatrix &A, double &lambda, int niters, double tolerance, bool verbose)
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
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.
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 main(int argc, char *argv[])
int MyPID() const
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
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...
double ElapsedTime(void) const
Epetra_Time elapsed time function.