63int checkValues(
double x,
double y,
string message =
"",
bool verbose =
false) {
64 if (fabs((x-y)/x) > 0.01) {
66 if (verbose) cout <<
"********** " << message <<
" check failed.********** " << endl;
69 if (verbose) cout << message <<
" check OK." << endl;
78 int globalbadvalue = 0;
79 for (
int j=0; j<numVectors; j++)
80 for (
int i=0; i< length; i++)
85 if (globalbadvalue==0) cout << message <<
" check OK." << endl;
86 else cout <<
"********* " << message <<
" check failed.********** " << endl;
88 return(globalbadvalue);
100 double * lambda,
int niters,
double tolerance,
103int main(
int argc,
char *argv[])
105 int ierr = 0, i, forierr = 0;
110 MPI_Init(&argc,&argv);
113 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
123 bool verbose =
false;
126 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
128 int verbose_int = verbose ? 1 : 0;
130 verbose = verbose_int==1 ? true :
false;
139 int MyPID = Comm.
MyPID();
142 if(verbose && MyPID==0)
145 if (verbose) cout <<
"Processor "<<MyPID<<
" of "<< NumProc
146 <<
" is alive."<<endl;
149 if(verbose && rank!=0)
152 int NumMyEquations = 10000;
153 int NumGlobalEquations = (NumMyEquations * NumProc) +
EPETRA_MIN(NumProc,3);
159 Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
168 vector<int> NumNz(NumMyEquations);
173 for(i = 0; i < NumMyEquations; i++)
174 if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
190 vector<double> Values(2);
193 vector<int> Indices(2);
198 for(i = 0; i < NumMyEquations; i++) {
199 if(MyGlobalElements[i] == 0) {
203 else if (MyGlobalElements[i] == NumGlobalEquations-1) {
204 Indices[0] = NumGlobalEquations-2;
208 Indices[0] = MyGlobalElements[i]-1;
209 Indices[1] = MyGlobalElements[i]+1;
212 forierr += !(A.
InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0])==0);
213 forierr += !(A.
InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i])>0);
241 if (verbose) cout <<
"=======================================" << endl
242 <<
"Testing Jad using CrsMatrix as input..." << endl
243 <<
"=======================================" << endl;
250 if (verbose) cout <<
"\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
256 vector<double> Rowvals(numvals);
257 vector<int> Rowinds(numvals);
260 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
268 if (verbose) cout <<
"================================================================" << endl
269 <<
"Testing Jad using Jad matrix as input matrix for construction..." << endl
270 <<
"================================================================" << endl;
288 double tolerance = 1.0e-2;
299 double elapsed_time = timer.
ElapsedTime() - startTime;
300 double total_flops = q.
Flops();
301 double MFLOPs = total_flops/elapsed_time/1000000.0;
302 double lambdaref = lambda;
303 double flopsref = total_flops;
306 cout <<
"\n\nTotal MFLOPs for reference first solve = " << MFLOPs << endl
307 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
313 total_flops = q.
Flops();
314 MFLOPs = total_flops/elapsed_time/1000000.0;
317 cout <<
"\n\nTotal MFLOPs for candidate first solve = " << MFLOPs << endl
318 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
327 if (verbose) cout <<
"\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
334 total_flops = q.
Flops();
335 MFLOPs = total_flops/elapsed_time/1000000.0;
337 flopsref = total_flops;
340 cout <<
"\n\nTotal MFLOPs for reference transpose solve = " << MFLOPs << endl
341 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
347 total_flops = q.
Flops();
348 MFLOPs = total_flops/elapsed_time/1000000.0;
351 cout <<
"\n\nTotal MFLOPs for candidate transpose solve = " << MFLOPs << endl
352 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
362 Epetra_Vector& resid,
double* lambda,
int niters,
double tolerance,
bool verbose)
369 double normz, residual;
373 for(
int iter = 0; iter < niters; iter++) {
375 q.
Scale(1.0/normz, z);
378 if(iter%100==0 || iter+1==niters) {
379 resid.
Update(1.0, z, -(*lambda), q, 0.0);
380 resid.
Norm2(&residual);
381 if(verbose) cout <<
"Iter = " << iter <<
" Lambda = " << *lambda
382 <<
" Residual of A*q - lambda*q = " << residual << endl;
384 if(residual < tolerance) {
506 for (
int j=0; j<nA; j++) {
507 double curVal = valuesA[j];
508 int curIndex = indicesA[j];
509 bool notfound =
true;
511 while (notfound && jj< nB) {
512 if (!
checkValues(curVal, valuesB[jj])) notfound =
false;
516 vector<int>::iterator p = find(indicesB.begin(),indicesB.end(),curIndex);
521 if (verbose) cout <<
"RowMatrix Methods check OK" << endl;
std::string Epetra_Version()
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
virtual int NumProc() const =0
Returns total number of processes.
virtual int MyPID() const =0
Return my process ID.
void ResetFlops() const
Resets the number of floating point operations to zero for this multi-vector.
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
double Flops() const
Returns the number of floating point operations with this multi-vector.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int NumGlobalEntries(long long Row) const
Returns the current number of nonzero entries in specified global row on this processor.
bool MyGlobalRow(int GID) const
Returns true of GID is owned by the calling processor, otherwise it returns false.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Replace specified existing values with this list of entries for a given global row of the matrix.
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified global row in user-provided arrays.
bool IndicesAreLocal() const
If matrix indices has been transformed to local, this query returns true, otherwise it returns false.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
bool IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
Epetra_Flops: The Epetra Floating Point Operations Class.
Epetra_JadMatrix: A class for constructing matrix objects optimized for common kernels.
int UpdateValues(const Epetra_RowMatrix &Matrix, bool CheckStructure=false)
Update values using a matrix with identical structure.
Epetra_Map: A class for partitioning vectors and matrices.
Epetra_MpiComm: The Epetra MPI Communication Class.
int Broadcast(double *MyVals, int Count, int Root) const
Epetra_MpiComm Broadcast function.
int NumProc() const
Returns total number of processes.
int MyPID() const
Return my process ID.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int Random()
Set multi-vector values to random numbers.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, transpose of this operator will be applied.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector X in Y.
virtual const Epetra_Comm & Comm() const =0
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual const Epetra_Map & OperatorDomainMap() const =0
Returns the Epetra_Map object associated with the domain of this operator.
virtual bool HasNormInf() const =0
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const Epetra_Map & OperatorRangeMap() const =0
Returns the Epetra_Map object associated with the range of this operator.
virtual bool UseTranspose() const =0
Returns the current UseTranspose setting.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
virtual int NumMyRows() const =0
Returns the number of matrix rows owned by the calling processor.
virtual int NumMyCols() const =0
Returns the number of matrix columns owned by the calling processor.
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const =0
Returns a copy of the main diagonal in a user-provided vector.
virtual int NumGlobalCols() const =0
Returns the number of global matrix columns.
virtual int NumGlobalNonzeros() const =0
Returns the number of nonzero entries in the global matrix.
virtual const Epetra_Map & RowMatrixColMap() const =0
Returns the Epetra_Map object associated with the columns of this matrix.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
virtual int NumMyNonzeros() const =0
Returns the number of nonzero entries in the calling processor's portion of the matrix.
virtual int NumMyDiagonals() const =0
Returns the number of local nonzero diagonal entries, based on global row/column index comparisons.
virtual int NumMyRowEntries(int MyRow, int &NumEntries) const =0
Returns the number of nonzero entries in MyRow.
virtual int NumGlobalRows() const =0
Returns the number of global matrix rows.
virtual bool LowerTriangular() const =0
If matrix is lower triangular in local index space, this query returns true, otherwise it returns fal...
virtual int NumGlobalDiagonals() const =0
Returns the number of global nonzero diagonal entries, based on global row/column index comparisons.
virtual int MaxNumEntries() const =0
Returns the maximum of NumMyRowEntries() over all rows.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
Returns a copy of the specified local row in user-provided arrays.
virtual int InvRowSums(Epetra_Vector &x) const =0
Computes the sum of absolute values of the rows of the Epetra_RowMatrix, results returned in x.
virtual bool UpperTriangular() const =0
If matrix is upper triangular in local index space, this query returns true, otherwise it returns fal...
virtual int RightScale(const Epetra_Vector &x)=0
Scales the Epetra_RowMatrix on the right with a Epetra_Vector x.
virtual double NormOne() const =0
Returns the one norm of the global matrix.
virtual double NormInf() const =0
Returns the infinity norm of the global matrix.
virtual bool Filled() const =0
If FillComplete() has been called, this query returns true, otherwise it returns false.
virtual int InvColSums(Epetra_Vector &x) const =0
Computes the sum of absolute values of the columns of the Epetra_RowMatrix, results returned in x.
virtual int LeftScale(const Epetra_Vector &x)=0
Scales the Epetra_RowMatrix on the left with a Epetra_Vector x.
Epetra_SerialComm: The Epetra Serial Communication Class.
virtual const Epetra_BlockMap & Map() const =0
Returns a reference to the Epetra_BlockMap for this object.
Epetra_Time: The Epetra Timing Class.
double ElapsedTime(void) const
Epetra_Time elapsed time function.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
#define EPETRA_TEST_ERR(a, b)
int main(int argc, char *argv[])
int powerMethodTests(Epetra_RowMatrix &A, Epetra_RowMatrix &JadA, Epetra_Map &Map, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, bool verbose)
int power_method(bool TransA, Epetra_RowMatrix &A, Epetra_Vector &q, Epetra_Vector &z0, Epetra_Vector &resid, double *lambda, int niters, double tolerance, bool verbose)
int check(Epetra_RowMatrix &A, Epetra_RowMatrix &B, bool verbose)
int checkValues(double x, double y, string message="", bool verbose=false)
int checkMultiVectors(Epetra_MultiVector &X, Epetra_MultiVector &Y, string message="", bool verbose=false)