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]);
177 ierr =
A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
182 ierr =
A.FillComplete();
190 int niters = (int) NumGlobalElements*10;
191 double tolerance = 1.0e-2;
195 A.SetFlopCounter(counter);
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" 210 if (
A.MyGlobalRow(0)) {
211 int numvals =
A.NumGlobalEntries(0);
212 std::vector<double> Rowvals(numvals);
213 std::vector<long long> Rowinds(numvals);
214 A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]);
215 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
217 A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[0]);
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);
269 A.Multiply(
false, q, 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) {
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.
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()
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)
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 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.
double ElapsedTime(void) const
Epetra_Time elapsed time function.