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();
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;
182 ierr =
A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
184 ierr =
A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i);
200 double tolerance = 1.0e-3;
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 long long * Rowinds =
new long long [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);
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);
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 Norm2(double *Result) const
int MyGlobalElements(int *MyGlobalElementList) const
int Dot(const Epetra_MultiVector &A, double *Result) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int main(int argc, char *argv[])
int FillComplete(bool OptimizeDataStorage=true)
int Scale(double ScalarValue)
int NumMyElements() const
std::string Ifpack_Version()
int power_method(Epetra_CrsMatrix &A, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, double *lambda, int niters, double tolerance, bool verbose)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
#define IFPACK_CHK_ERR(ifpack_err)
void ResetStartTime(void)
double ElapsedTime(void) const