50 #include "Epetra_Comm.h"
51 #include "Epetra_Map.h"
52 #include "Epetra_Time.h"
53 #include "Epetra_CrsMatrix.h"
54 #include "Epetra_Vector.h"
55 #include "Epetra_Object.h"
59 #include "Epetra_MpiComm.h"
61 #include "Epetra_SerialComm.h"
70 double * lambda,
int niters,
double tolerance,
74 int main(
int argc,
char *argv[])
85 MPI_Init(&argc,&argv);
88 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
101 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
108 int MyPID = Comm.
MyPID();
111 if (verbose && MyPID==0)
114 if (verbose) cout <<
"Processor "<<MyPID<<
" of "<< NumProc
115 <<
" is alive."<<endl;
120 if (verbose && rank!=0) verbose =
false;
122 long long NumMyEquations = 10000;
123 long long NumGlobalEquations = ((
long long)NumMyEquations)*NumProc+
EPETRA_MIN(NumProc,3);
124 if (MyPID < 3) NumMyEquations++;
128 Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
131 long long * MyGlobalElements =
new long long[Map.
NumMyElements()];
137 int * NumNz =
new int[NumMyEquations];
142 for (i=0; i<NumMyEquations; i++)
143 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalEquations-1)
157 double *Values =
new double[2];
158 Values[0] = -1.0; Values[1] = -1.0;
159 long long *Indices =
new long long[2];
163 for (i=0; i<NumMyEquations; i++)
165 if (MyGlobalElements[i]==0)
170 else if (MyGlobalElements[i] == NumGlobalEquations-1)
172 Indices[0] = NumGlobalEquations-2;
177 Indices[0] = MyGlobalElements[i]-1;
178 Indices[1] = MyGlobalElements[i]+1;
200 double tolerance = 1.0e-3;
204 ierr +=
power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
206 double total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
207 double MFLOPs = total_flops/elapsed_time/1000000.0;
209 if (verbose) cout <<
"\n\nTotal MFLOPs for first solve = " << MFLOPs << endl<< endl;
212 if (verbose) cout <<
"\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
218 double * Rowvals =
new double [numvals];
219 long long * Rowinds =
new long long [numvals];
222 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
229 A.ResetFlops(); q.ResetFlops(); z.ResetFlops(); resid.ResetFlops();
230 ierr +=
power_method(A, q, z, resid, &lambda, niters, tolerance, verbose);
232 total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
233 MFLOPs = total_flops/elapsed_time/1000000.0;
235 if (verbose) cout <<
"\n\nTotal MFLOPs for second solve = " << MFLOPs << endl<< endl;
242 delete [] MyGlobalElements;
250 int NumMyElements1 = 2;
251 int NumMyEquations1 = NumMyElements1;
252 long long NumGlobalEquations1 = ((
long long)NumMyEquations1)*NumProc;
254 Epetra_Map Map1(-1LL, NumMyElements1, 0, Comm);
257 long long * MyGlobalElements1 =
new long long[Map1.
NumMyElements()];
263 int * NumNz1 =
new int[NumMyEquations1];
268 for (i=0; i<NumMyEquations1; i++)
269 if (MyGlobalElements1[i]==0 || MyGlobalElements1[i] == NumGlobalEquations1-1)
283 double *Values1 =
new double[2];
284 Values1[0] = -1.0; Values1[1] = -1.0;
285 long long *Indices1 =
new long long[2];
289 for (i=0; i<NumMyEquations1; i++)
291 if (MyGlobalElements1[i]==0)
296 else if (MyGlobalElements1[i] == NumGlobalEquations1-1)
298 Indices1[0] = NumGlobalEquations1-2;
303 Indices1[0] = MyGlobalElements1[i]-1;
304 Indices1[1] = MyGlobalElements1[i]+1;
317 if (verbose) cout <<
"\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
324 delete [] MyGlobalElements1;
341 double * lambda,
int niters,
double tolerance,
350 double normz, residual;
354 for (
int iter = 0;
iter < niters;
iter++)
357 q.Scale(1.0/normz, z);
360 if (iter%100==0 || iter+1==niters)
362 resid.Update(1.0, z, -(*lambda), q, 0.0);
363 resid.Norm2(&residual);
364 if (verbose) cout <<
"Iter = " <<
iter <<
" Lambda = " << *lambda
365 <<
" Residual of A*q - lambda*q = " << residual << endl;
367 if (residual < tolerance) {
bool MyGlobalRow(int GID) const
int NumGlobalEntries(long long Row) const
int power_method(Epetra_CrsMatrix &A, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, double *lambda, int niters, double tolerance, bool verbose)
int ExtractGlobalRowCopy(int_type Row, int Length, int &NumEntries, double *values, int_type *Indices) const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
int MyGlobalElements(int *MyGlobalElementList) const
double ElapsedTime(void) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int FillComplete(bool OptimizeDataStorage=true)
int NumMyElements() const
std::string Ifpack_Version()
int main(int argc, char *argv[])
#define IFPACK_CHK_ERR(ifpack_err)
void ResetStartTime(void)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)