42 #include "Epetra_Map.h" 43 #include "Epetra_Time.h" 44 #include "Epetra_CrsMatrix.h" 46 #include "HYPRE_IJ_mv.h" 48 #include "Epetra_Vector.h" 49 #include "Epetra_Flops.h" 51 #include "Epetra_MpiComm.h" 54 #include "Epetra_SerialComm.h" 56 #include "../epetra_test_err.h" 57 #include "Epetra_Version.h" 66 int checkValues(
double x,
double y,
string message =
"",
bool verbose =
false) {
67 if (fabs((x-y)/x) > 0.01) {
69 if (verbose) cout <<
"********** " << message <<
" check failed.********** " << endl;
72 if (verbose) cout << message <<
" check OK." << endl;
77 int checkMultiVectors( Epetra_MultiVector & X, Epetra_MultiVector & Y,
string message =
"",
bool verbose =
false) {
78 int numVectors = X.NumVectors();
79 int length = Y.MyLength();
81 int globalbadvalue = 0;
82 for (
int j=0; j<numVectors; j++)
83 for (
int i=0; i< length; i++)
85 X.Map().Comm().MaxAll(&badvalue, &globalbadvalue, 1);
88 if (globalbadvalue==0) cout << message <<
" check OK." << endl;
89 else cout <<
"********* " << message <<
" check failed.********** " << endl;
91 return(globalbadvalue);
94 int check(Epetra_RowMatrix & A, Epetra_RowMatrix & B,
bool verbose);
96 int powerMethodTests(Epetra_RowMatrix & A, Epetra_RowMatrix & JadA, Epetra_Map & Map,
97 Epetra_Vector & q, Epetra_Vector & z, Epetra_Vector & resid,
bool verbose);
102 Epetra_Vector& resid,
103 double * lambda,
int niters,
double tolerance,
106 int main(
int argc,
char *argv[])
108 int ierr = 0, i, forierr = 0;
113 MPI_Init(&argc,&argv);
116 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
117 Epetra_MpiComm Comm( MPI_COMM_WORLD );
122 Epetra_SerialComm Comm;
126 bool verbose =
false;
129 if (argc>1)
if (argv[1][0]==
'-' && argv[1][1]==
'v') verbose =
true;
131 int verbose_int = verbose ? 1 : 0;
132 Comm.Broadcast(&verbose_int, 1, 0);
133 verbose = verbose_int==1 ? true :
false;
141 Comm.SetTracebackMode(0);
142 int MyPID = Comm.MyPID();
143 int NumProc = Comm.NumProc();
145 if(verbose && MyPID==0)
146 cout << Epetra_Version() << endl << endl;
148 if (verbose) cout <<
"Processor "<<MyPID<<
" of "<< NumProc
149 <<
" is alive."<<endl;
152 if(verbose && rank!=0)
155 int NumMyEquations = 10000;
156 int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);
162 Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
165 vector<int> MyGlobalElements(Map.NumMyElements());
166 Map.MyGlobalElements(&MyGlobalElements[0]);
171 vector<int> NumNz(NumMyEquations);
176 for(i = 0; i < NumMyEquations; i++)
177 if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
184 Epetra_CrsMatrix A(Copy, Map, &NumNz[0]);
193 vector<double> Values(2);
196 vector<int> Indices(2);
201 for(i = 0; i < NumMyEquations; i++) {
202 if(MyGlobalElements[i] == 0) {
206 else if (MyGlobalElements[i] == NumGlobalEquations-1) {
207 Indices[0] = NumGlobalEquations-2;
211 Indices[0] = MyGlobalElements[i]-1;
212 Indices[1] = MyGlobalElements[i]+1;
215 forierr += !(A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0])==0);
216 forierr += !(A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i])>0);
224 HYPRE_IJMatrix Matrix;
225 int ilower = Map.MinMyGID();
226 int iupper = Map.MaxMyGID();
229 HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &Matrix);
230 HYPRE_IJMatrixSetObjectType(Matrix, HYPRE_PARCSR);
231 HYPRE_IJMatrixInitialize(Matrix);
233 for(i = 0; i < A.NumMyRows(); i++){
235 A.NumMyRowEntries(i, numElements);
236 vector<int> my_indices; my_indices.resize(numElements);
237 vector<double> my_values; my_values.resize(numElements);
239 A.ExtractMyRowCopy(i, numElements, numEntries, &my_values[0], &my_indices[0]);
240 for(
int j = 0; j < numEntries; j++) {
241 my_indices[j] = A.GCID(my_indices[j]);
244 GlobalRow[0] = A.GRID(i);
245 HYPRE_IJMatrixSetValues(Matrix, 1, &numEntries, GlobalRow, &my_indices[0], &my_values[0]);
247 HYPRE_IJMatrixAssemble(Matrix);
251 JadA.SetMaps(JadA.RowMatrixRowMap(), A.RowMatrixColMap());
254 Epetra_Vector q(Map);
255 Epetra_Vector z(Map); z.Random();
256 Epetra_Vector resid(Map);
258 Epetra_Flops flopcounter;
259 A.SetFlopCounter(flopcounter);
262 resid.SetFlopCounter(A);
263 JadA.SetFlopCounter(A);
265 if (verbose) cout <<
"=======================================" << endl
266 <<
"Testing Jad using CrsMatrix as input..." << endl
267 <<
"=======================================" << endl;
280 Epetra_Vector & q, Epetra_Vector & z, Epetra_Vector & resid,
bool verbose) {
286 double tolerance = 1.0e-2;
293 Epetra_Time timer(Map.Comm());
295 double startTime = timer.ElapsedTime();
297 double elapsed_time = timer.ElapsedTime() - startTime;
298 double total_flops = q.Flops();
299 double MFLOPs = total_flops/elapsed_time/1000000.0;
300 double lambdaref = lambda;
301 double flopsref = total_flops;
304 cout <<
"\n\nTotal MFLOPs for reference first solve = " << MFLOPs << endl
305 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
308 startTime = timer.ElapsedTime();
311 elapsed_time = timer.ElapsedTime() - startTime;
312 total_flops = q.Flops();
313 MFLOPs = total_flops/elapsed_time/1000000.0;
316 cout <<
"\n\nTotal MFLOPs for candidate first solve = " << MFLOPs << endl
317 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
326 if (verbose) cout <<
"\n\nUsing transpose of matrix and solving again (should give same result).\n\n" 330 startTime = timer.ElapsedTime();
332 elapsed_time = timer.ElapsedTime() - startTime;
333 total_flops = q.Flops();
334 MFLOPs = total_flops/elapsed_time/1000000.0;
336 flopsref = total_flops;
339 cout <<
"\n\nTotal MFLOPs for reference transpose solve = " << MFLOPs << endl
340 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
343 startTime = timer.ElapsedTime();
345 elapsed_time = timer.ElapsedTime() - startTime;
346 total_flops = q.Flops();
347 MFLOPs = total_flops/elapsed_time/1000000.0;
350 cout <<
"\n\nTotal MFLOPs for candidate transpose solve = " << MFLOPs << endl
351 <<
"Total FLOPS = " <<total_flops <<endl<<endl;
360 int power_method(
bool TransA, Epetra_RowMatrix& A, Epetra_Vector& q, Epetra_Vector& z0,
361 Epetra_Vector& resid,
double* lambda,
int niters,
double tolerance,
bool verbose)
368 double normz, residual;
372 for(
int iter = 0; iter < niters; iter++) {
374 q.Scale(1.0/normz, z);
375 A.Multiply(TransA, q, z);
377 if(iter%100==0 || iter+1==niters) {
378 resid.Update(1.0, z, -(*lambda), q, 0.0);
379 resid.Norm2(&residual);
380 if(verbose) cout <<
"Iter = " << iter <<
" Lambda = " << *lambda
381 <<
" Residual of A*q - lambda*q = " << residual << endl;
383 if(residual < tolerance) {
391 int check(Epetra_RowMatrix& A, Epetra_RowMatrix & B,
bool verbose) {
408 for (
int i=0; i<A.NumMyRows(); i++) {
410 A.NumMyRowEntries(i,nA); B.NumMyRowEntries(i,nB);
414 EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
415 EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
417 EPETRA_TEST_ERR(!A.RowMatrixColMap().SameAs(B.RowMatrixColMap()),ierr);
418 EPETRA_TEST_ERR(!A.RowMatrixRowMap().SameAs(B.RowMatrixRowMap()),ierr);
424 Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
425 Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
426 Epetra_MultiVector YA2(YA1);
427 Epetra_MultiVector YB1(YA1);
428 Epetra_MultiVector YB2(YA1);
432 A.SetUseTranspose(transA);
433 B.SetUseTranspose(transA);
435 A.Multiply(transA, X, YA2);
439 B.Multiply(transA, X, YB2);
444 Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
445 Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
446 Epetra_MultiVector YA2(YA1);
447 Epetra_MultiVector YB1(YA1);
448 Epetra_MultiVector YB2(YA1);
452 A.SetUseTranspose(transA);
453 B.SetUseTranspose(transA);
455 A.Multiply(transA, X, YA2);
459 B.Multiply(transA, X,YB2);
464 Epetra_Vector diagA(A.RowMatrixRowMap());
466 Epetra_Vector diagB(B.RowMatrixRowMap());
470 Epetra_Vector rowA(A.RowMatrixRowMap());
472 Epetra_Vector rowB(B.RowMatrixRowMap());
476 Epetra_Vector colA(A.RowMatrixColMap());
478 Epetra_Vector colB(B.RowMatrixColMap());
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;
#define EPETRA_TEST_ERR(a, b)
int powerMethodTests(Epetra_RowMatrix &A, Epetra_RowMatrix &JadA, Epetra_Map &Map, Epetra_Vector &q, Epetra_Vector &z, Epetra_Vector &resid, bool verbose)
int checkValues(double x, double y, string message="", bool verbose=false)
int main(int argc, char *argv[])
int check(Epetra_RowMatrix &A, Epetra_RowMatrix &B, bool verbose)
int checkMultiVectors(Epetra_MultiVector &X, Epetra_MultiVector &Y, string message="", bool verbose=false)
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)