71 int main(
int argc,
char *argv[])
79 MPI_Init(&argc,&argv);
89 int MyPID = Comm.
MyPID();
91 bool verbose = (MyPID==0);
96 std::cout << Comm << std::endl;
102 std::cout <<
"Usage: " << argv[0] <<
" number_of_equations" << std::endl;
105 int NumGlobalElements = std::atoi(argv[1]);
107 if (NumGlobalElements < NumProc)
110 std::cout <<
"numGlobalBlocks = " << NumGlobalElements
111 <<
" cannot be < number of processors = " << NumProc << std::endl;
124 std::vector<int> MyGlobalElements(NumMyElements);
131 std::vector<int> NumNz(NumMyElements);
136 for (i=0; i<NumMyElements; i++)
137 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalElements-1)
151 std::vector<double> Values(2);
152 Values[0] = -1.0; Values[1] = -1.0;
153 std::vector<int> Indices(2);
157 for (i=0; i<NumMyElements; i++)
159 if (MyGlobalElements[i]==0)
164 else if (MyGlobalElements[i] == NumGlobalElements-1)
166 Indices[0] = NumGlobalElements-2;
171 Indices[0] = MyGlobalElements[i]-1;
172 Indices[1] = MyGlobalElements[i]+1;
175 ierr =
A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
178 ierr =
A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
183 ierr =
A.FillComplete();
191 int niters = NumGlobalElements*10;
192 double tolerance = 1.0e-2;
196 A.SetFlopCounter(counter);
200 double total_flops =counter.
Flops();
201 double MFLOPs = total_flops/elapsed_time/1000000.0;
204 std::cout <<
"\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;
208 std::cout <<
"\nIncreasing magnitude of first diagonal term, solving again\n\n" 211 if (
A.MyGlobalRow(0)) {
212 int numvals =
A.NumGlobalEntries(0);
213 std::vector<double> Rowvals(numvals);
214 std::vector<int> Rowinds(numvals);
216 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
227 total_flops = counter.
Flops();
228 MFLOPs = total_flops/elapsed_time/1000000.0;
231 std::cout <<
"\n\nTotal MFLOPs for second solve = " << MFLOPs << std::endl<< std::endl;
245 double tolerance,
bool verbose) {
255 resid.SetFlopCounter(
A);
262 double normz, residual;
266 for (
int iter = 0; iter < niters; iter++)
269 q.Scale(1.0/normz, z);
270 A.Multiply(
false, q, z);
272 if (iter%100==0 || iter+1==niters)
274 resid.Update(1.0, z, -lambda, q, 0.0);
275 resid.Norm2(&residual);
276 if (verbose) std::cout <<
"Iter = " << iter <<
" Lambda = " << lambda
277 <<
" Residual of A*q - lambda*q = " 278 << residual << std::endl;
280 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.
int main(int argc, char *argv[])
T * Epetra_Util_data_ptr(std::vector< T > &vec)
Function that returns either a pointer to the first entry in the vector or, if the vector is empty...
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.
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.
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 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.