74int 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();
114 if (
verbose) cout <<
"Processor "<<MyPID<<
" of "<< NumProc
115 <<
" is alive."<<endl;
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;
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;
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()];
258 Map1.MyGlobalElements(MyGlobalElements1);
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;
308 ierr = A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1);
310 ierr = A1.InsertGlobalValues(MyGlobalElements1[i], 1, &two1, MyGlobalElements1+i);
317 if (
verbose) cout <<
"\n\nPrint out tridiagonal matrix, each part on each processor.\n\n" << endl;
324 delete [] MyGlobalElements1;
int main(int argc, char *argv[])
int power_method(Epetra_CrsMatrix &A, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, double *lambda, int niters, double tolerance, bool verbose)