55 #include "../epetra_test_err.h" 63 int 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,
103 int 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);
232 A.SetFlopCounter(flopcounter);
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" 254 if (
A.MyGlobalRow(0)) {
255 int numvals =
A.NumGlobalEntries(0);
256 vector<double> Rowvals(numvals);
257 vector<int> Rowinds(numvals);
258 A.ExtractGlobalRowCopy(0, numvals, numvals, &Rowvals[0], &Rowinds[0]);
260 for (i=0; i<numvals; i++)
if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
262 A.ReplaceGlobalValues(0, numvals, &Rowvals[0], &Rowinds[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;
310 startTime = timer.ElapsedTime();
312 elapsed_time = timer.ElapsedTime() - startTime;
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" 331 startTime = timer.ElapsedTime();
333 elapsed_time = timer.ElapsedTime() - startTime;
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;
344 startTime = timer.ElapsedTime();
346 elapsed_time = timer.ElapsedTime() - startTime;
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);
376 A.Multiply(TransA, q, 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) {
409 for (
int i=0; i<
A.NumMyRows(); i++) {
411 A.NumMyRowEntries(i,nA);
B.NumMyRowEntries(i,nB);
432 A.SetUseTranspose(transA);
433 B.SetUseTranspose(transA);
435 A.Multiply(transA, X, YA2);
439 B.Multiply(transA, X, YB2);
452 A.SetUseTranspose(transA);
453 B.SetUseTranspose(transA);
455 A.Multiply(transA, X, YA2);
459 B.Multiply(transA, X,YB2);
496 vector<double> valuesA(
A.MaxNumEntries());
497 vector<int> indicesA(
A.MaxNumEntries());
498 vector<double> valuesB(
B.MaxNumEntries());
499 vector<int> indicesB(
B.MaxNumEntries());
501 for (
int i=0; i<
A.NumMyRows(); i++) {
503 EPETRA_TEST_ERR(
A.ExtractMyRowCopy(i,
A.MaxNumEntries(), nA, &valuesA[0], &indicesA[0]),ierr);
504 EPETRA_TEST_ERR(
B.ExtractMyRowCopy(i,
B.MaxNumEntries(), nB, &valuesB[0], &indicesB[0]),ierr);
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;
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
void SetFlopCounter(const Epetra_Flops &FlopCounter_in)
Set the internal Epetra_Flops() pointer.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
Epetra_Map: A class for partitioning vectors and matrices.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int Random()
Set multi-vector values to random numbers.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
#define EPETRA_TEST_ERR(a, b)
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
Epetra_MpiComm: The Epetra MPI Communication Class.
std::string Epetra_Version()
Epetra_JadMatrix: A class for constructing matrix objects optimized for common kernels.
Epetra_Time: The Epetra Timing Class.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
int NumProc() const
Returns total number of processes (always returns 1 for SerialComm).
int powerMethodTests(Epetra_RowMatrix &A, Epetra_RowMatrix &JadA, Epetra_Map &Map, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, bool verbose)
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
double Flops() const
Returns the number of floating point operations with this multi-vector.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int Broadcast(double *MyVals, int Count, int Root) const
Epetra_SerialComm Broadcast function.
int NumMyElements() const
Number of elements on the calling processor.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
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 UpdateValues(const Epetra_RowMatrix &Matrix, bool CheckStructure=false)
Update values using a matrix with identical structure.
int main(int argc, char *argv[])
Epetra_SerialComm: The Epetra Serial Communication Class.
Epetra_Flops: The Epetra Floating Point Operations Class.
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 check(Epetra_RowMatrix &A, Epetra_RowMatrix &B, bool verbose)
int MyPID() const
Return my process ID.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices...
int checkValues(double x, double y, string message="", bool verbose=false)
int checkMultiVectors(Epetra_MultiVector &X, Epetra_MultiVector &Y, string message="", bool verbose=false)
double ElapsedTime(void) const
Epetra_Time elapsed time function.
void ResetFlops() const
Resets the number of floating point operations to zero for this multi-vector.