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,
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 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);
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 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;
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"
218 double * Rowvals =
new double [numvals];
219 int * Rowinds =
new int [numvals];
222 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
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;
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 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;
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) {
#define IFPACK_CHK_ERR(ifpack_err)
std::string Ifpack_Version()
int MyGlobalElements(int *MyGlobalElementList) const
int NumMyElements() const
int NumGlobalEntries(long long Row) const
bool MyGlobalRow(int GID) const
int FillComplete(bool OptimizeDataStorage=true)
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
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
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int Scale(double ScalarValue)
int Dot(const Epetra_MultiVector &A, double *Result) const
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Norm2(double *Result) const
void ResetStartTime(void)
double ElapsedTime(void) const
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)