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);
89 Epetra_MpiComm Comm( MPI_COMM_WORLD );
94 Epetra_SerialComm Comm;
101 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v')
verbose =
true;
108 int MyPID = Comm.MyPID();
109 int NumProc = Comm.NumProc();
114 if (
verbose) cout <<
"Processor "<<MyPID<<
" of "<< NumProc
115 <<
" is alive."<<endl;
122 int NumMyEquations = 10000;
123 int NumGlobalEquations = NumMyEquations*NumProc+EPETRA_MIN(NumProc,3);
124 if (MyPID < 3) NumMyEquations++;
128 Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
131 int * MyGlobalElements =
new int[Map.NumMyElements()];
132 Map.MyGlobalElements(MyGlobalElements);
137 int * NumNz =
new int[NumMyEquations];
142 for (i=0; i<NumMyEquations; i++)
143 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == NumGlobalEquations-1)
150 Epetra_CrsMatrix A(Copy, Map, NumNz);
157 double *Values =
new double[2];
158 Values[0] = -1.0; Values[1] = -1.0;
159 int *Indices =
new int[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;
182 ierr2 = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
184 ierr2 = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i);
193 Epetra_Vector q(Map);
194 Epetra_Vector z(Map);
195 Epetra_Vector resid(Map);
200 double tolerance = 1.0e-3;
203 Epetra_Time timer(Comm);
205 double elapsed_time = timer.ElapsedTime();
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" 216 if (A.MyGlobalRow(0)) {
217 int numvals = A.NumGlobalEntries(0);
218 double * Rowvals =
new double [numvals];
219 int * Rowinds =
new int [numvals];
220 A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds);
222 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
224 A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
228 timer.ResetStartTime();
229 A.ResetFlops(); q.ResetFlops(); z.ResetFlops(); resid.ResetFlops();
231 elapsed_time = timer.ElapsedTime();
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 int NumGlobalEquations1 = NumMyEquations1*NumProc;
254 Epetra_Map Map1(-1, NumMyElements1, 0, Comm);
257 int * MyGlobalElements1 =
new int[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)
276 Epetra_CrsMatrix A1(Copy, Map1, NumNz1);
283 double *Values1 =
new double[2];
284 Values1[0] = -1.0; Values1[1] = -1.0;
285 int *Indices1 =
new int[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 ierr2 = A1.InsertGlobalValues(MyGlobalElements1[i], NumEntries1, Values1, Indices1);
310 ierr2 = 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;
340 Epetra_Vector& resid,
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);
358 A.Multiply(
false, q, 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) {
int power_method(Epetra_CrsMatrix &A, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, double *lambda, int niters, double tolerance, bool verbose)
int main(int argc, char *argv[])
std::string Ifpack_Version()
#define IFPACK_CHK_ERR(ifpack_err)