44 #include "Epetra_ConfigDefs.h" 48 #include "Epetra_MpiComm.h" 50 #include "Epetra_SerialComm.h" 52 #include "Trilinos_Util.h" 53 #include "Epetra_Comm.h" 54 #include "Epetra_Map.h" 55 #include "Epetra_Time.h" 56 #include "Epetra_BlockMap.h" 57 #include "Epetra_MultiVector.h" 58 #include "Epetra_Vector.h" 59 #include "Epetra_Export.h" 61 #include "Epetra_VbrMatrix.h" 62 #include "Epetra_CrsMatrix.h" 85 int checkValues(
double x,
double y, std::string message =
"",
bool verbose =
false) {
86 if (fabs((x-y)/x) > 0.01 && x > 1.0e-12) {
87 if (verbose) std::cout <<
"********** " << message <<
" check failed.********** " << std::endl;
91 if (verbose) std::cout << message <<
" check OK." << std::endl;
95 int runTests(Epetra_Map & map, Epetra_CrsMatrix & A, Epetra_Vector & x, Epetra_Vector & b, Epetra_Vector & xexact,
bool verbose);
100 int main(
int argc,
char *argv[]) {
103 MPI_Init(&argc,&argv);
104 Epetra_MpiComm comm (MPI_COMM_WORLD);
106 Epetra_SerialComm comm;
109 int MyPID = comm.MyPID();
111 bool verbose =
false;
112 bool verbose1 =
false;
115 if ((argv[1][0] ==
'-') && (argv[1][1] ==
'v')) {
117 if (MyPID==0) verbose =
true;
123 if (verbose1) std::cout << comm << std::endl;
132 Epetra_CrsMatrix * A;
135 Epetra_Vector * xexact;
137 int nx = 20*comm.NumProc();
140 int xoff[] = {-1, 0, 1, -1, 0, 1, 0};
141 int yoff[] = {-1, -1, -1, 0, 0, 0, 1};
146 Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
148 ierr +=
runTests(*map, *A, *x, *b, *xexact, verbose);
157 Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, 1);
159 ierr +=
runTests(*map, *A, *x, *b, *xexact, verbose);
168 Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, -1);
170 ierr +=
runTests(*map, *A, *x, *b, *xexact, verbose);
199 Epetra_Vector X(A.RowMap());
200 EPETRA_CHK_ERR(X.Random());
201 Epetra_Vector Y(A.RowMap());
202 EPETRA_CHK_ERR(A.Multiply(
false, X, Y));
207 int runTests(Epetra_Map & map, Epetra_CrsMatrix & A, Epetra_Vector & x, Epetra_Vector & b, Epetra_Vector & xexact,
bool verbose) {
212 Epetra_MultiVector X( map, 2,
false );
213 Epetra_MultiVector B( map, 2,
false );
214 Epetra_MultiVector Xexact( map, 2,
false );
216 for (
int i=0; i<X.NumVectors(); ++i) {
222 std::vector<double> residualmv(2);
223 residual = A.NormInf();
double rAInf = residual;
224 if (verbose) std::cout <<
"Inf Norm of A = " << residual << std::endl;
225 residual = A.NormOne();
double rAOne = residual;
226 if (verbose) std::cout <<
"One Norm of A = " << residual << std::endl;
227 xexact.Norm2(&residual);
double rxx = residual;
228 Xexact.Norm2(&residualmv[0]); std::vector<double> rXX( residualmv );
229 if (verbose) std::cout <<
"Norm of xexact = " << residual << std::endl;
230 if (verbose) std::cout <<
"Norm of Xexact = (" << residualmv[0] <<
", " <<residualmv[1] <<
")"<< std::endl;
231 Epetra_Vector tmp1(map);
232 Epetra_MultiVector tmp1mv(map,2,
false);
233 A.Multiply(
false, xexact, tmp1);
234 A.Multiply(
false, Xexact, tmp1mv);
235 tmp1.Norm2(&residual);
double rAx = residual;
236 tmp1mv.Norm2(&residualmv[0]); std::vector<double> rAX( residualmv );
237 if (verbose) std::cout <<
"Norm of Ax = " << residual << std::endl;
238 if (verbose) std::cout <<
"Norm of AX = (" << residualmv[0] <<
", " << residualmv[1] <<
")"<< std::endl;
239 b.Norm2(&residual);
double rb = residual;
240 B.Norm2(&residualmv[0]); std::vector<double> rB( residualmv );
241 if (verbose) std::cout <<
"Norm of b (should equal norm of Ax) = " << residual << std::endl;
242 if (verbose) std::cout <<
"Norm of B (should equal norm of AX) = (" << residualmv[0] <<
", " << residualmv[1] <<
")"<< std::endl;
243 tmp1.Update(1.0, b, -1.0);
244 tmp1mv.Update(1.0, B, -1.0);
245 tmp1.Norm2(&residual);
246 tmp1mv.Norm2(&residualmv[0]);
247 if (verbose) std::cout <<
"Norm of difference between compute Ax and Ax from file = " << residual << std::endl;
248 if (verbose) std::cout <<
"Norm of difference between compute AX and AX from file = (" << residualmv[0] <<
", " << residualmv[1] <<
")"<< std::endl;
249 map.Comm().Barrier();
252 "This is the official EpetraExt test map generated by the EpetraExt regression tests"));
255 "This is the official EpetraExt test matrix generated by the EpetraExt regression tests"));
258 "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));
261 "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));
264 "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));
267 "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));
270 "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));
273 "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));
280 Epetra_CrsMatrix * A1;
281 Epetra_CrsMatrix * A2;
282 Epetra_CrsMatrix * A3;
285 Epetra_Vector * xexact1;
286 Epetra_MultiVector * X1;
287 Epetra_MultiVector * B1;
288 Epetra_MultiVector * Xexact1;
292 if (map.SameAs(*map1)) {
293 if (verbose) std::cout <<
"Maps are equal. In/Out works." << std::endl;
296 if (verbose) std::cout <<
"Maps are not equal. In/Out fails." << std::endl;
310 residual = A1->NormInf();
double rA1Inf = residual;
311 if (verbose) std::cout <<
"Inf Norm of A1 = " << residual << std::endl;
312 ierr +=
checkValues(rA1Inf,rAInf,
"Inf Norm of A", verbose);
314 residual = A1->NormOne();
double rA1One = residual;
315 if (verbose) std::cout <<
"One Norm of A1 = " << residual << std::endl;
316 ierr +=
checkValues(rA1One,rAOne,
"One Norm of A", verbose);
318 xexact1->Norm2(&residual);
double rxx1 = residual;
319 if (verbose) std::cout <<
"Norm of xexact1 = " << residual << std::endl;
320 ierr +=
checkValues(rxx1,rxx,
"Norm of xexact", verbose);
322 Xexact1->Norm2(&residualmv[0]); std::vector<double> rXX1(residualmv);
323 if (verbose) std::cout <<
"Norm of Xexact1 = (" << residualmv[0] <<
", " <<residualmv[1]<<
")"<< std::endl;
324 ierr +=
checkValues(rXX1[0],rXX[0],
"Norm of Xexact", verbose);
325 ierr +=
checkValues(rXX1[1],rXX[1],
"Norm of Xexact", verbose);
327 Epetra_Vector tmp11(*map1);
328 A1->Multiply(
false, *xexact1, tmp11);
330 Epetra_MultiVector tmp11mv(*map1,2,
false);
331 A1->Multiply(
false, *Xexact1, tmp11mv);
333 tmp11.Norm2(&residual);
double rAx1 = residual;
334 if (verbose) std::cout <<
"Norm of A1*x1 = " << residual << std::endl;
335 ierr +=
checkValues(rAx1,rAx,
"Norm of A1*x", verbose);
337 tmp11mv.Norm2(&residualmv[0]); std::vector<double> rAX1(residualmv);
338 if (verbose) std::cout <<
"Norm of A1*X1 = (" << residualmv[0] <<
", "<<residualmv[1]<<
")"<< std::endl;
339 ierr +=
checkValues(rAX1[0],rAX[0],
"Norm of A1*X", verbose);
340 ierr +=
checkValues(rAX1[1],rAX[1],
"Norm of A1*X", verbose);
342 if (map1->IndexBase()==0) {
343 Epetra_Vector tmp12(*map1);
344 A2->Multiply(
false, *xexact1, tmp12);
346 tmp12.Norm2(&residual);
double rAx2 = residual;
347 if (verbose) std::cout <<
"Norm of A2*x1 = " << residual << std::endl;
348 ierr +=
checkValues(rAx2,rAx,
"Norm of A2*x", verbose);
350 Epetra_Vector tmp13(*map1);
351 A3->Multiply(
false, *xexact1, tmp13);
353 tmp13.Norm2(&residual);
double rAx3 = residual;
354 if (verbose) std::cout <<
"Norm of A3*x1 = " << residual << std::endl;
355 ierr +=
checkValues(rAx3,rAx,
"Norm of A3*x", verbose);
357 b1->Norm2(&residual);
double rb1 = residual;
358 if (verbose) std::cout <<
"Norm of b1 (should equal norm of Ax) = " << residual << std::endl;
361 B1->Norm2(&residualmv[0]); std::vector<double> rB1(residualmv);
362 if (verbose) std::cout <<
"Norm of B1 (should equal norm of AX) = (" << residualmv[0] <<
", "<<residualmv[1]<<
")"<< std::endl;
363 ierr +=
checkValues(rB1[0],rB[0],
"Norm of B", verbose);
364 ierr +=
checkValues(rB1[1],rB[1],
"Norm of B", verbose);
366 tmp11.Update(1.0, *b1, -1.0);
367 tmp11.Norm2(&residual);
368 if (verbose) std::cout <<
"Norm of difference between computed A1x1 and A1x1 from file = " << residual << std::endl;
369 ierr +=
checkValues(residual,0.0,
"Norm of difference between computed A1x1 and A1x1 from file", verbose);
371 tmp11mv.Update(1.0, *B1, -1.0);
372 tmp11mv.Norm2(&residualmv[0]);
373 if (verbose) std::cout <<
"Norm of difference between computed A1X1 and A1X1 from file = (" << residualmv[0] <<
", "<<residualmv[1]<<
")"<< std::endl;
374 ierr +=
checkValues(residualmv[0],0.0,
"Norm of difference between computed A1X1 and A1X1 from file", verbose);
375 ierr +=
checkValues(residualmv[1],0.0,
"Norm of difference between computed A1X1 and A1X1 from file", verbose);
377 if (map1->IndexBase()==0) {
delete A2;
delete A3;}
397 "This is the official EpetraExt test operator generated by the EpetraExt regression tests"));
400 A.OperatorRangeMap().Comm().Barrier();
401 A.OperatorRangeMap().Comm().Barrier();
402 Epetra_CrsMatrix * A1;
403 Epetra_CrsMatrix * A2;
408 residual = A.NormInf();
double rAInf = residual;
409 if (verbose) std::cout <<
"Inf Norm of Operator A = " << residual << std::endl;
410 residual = A1->NormInf();
double rA1Inf = residual;
411 if (verbose) std::cout <<
"Inf Norm of Matrix A1 = " << residual << std::endl;
412 ierr +=
checkValues(rA1Inf,rAInf,
"Inf Norm of A", verbose);
415 Epetra_Vector x(A.OperatorDomainMap()); x.Random();
416 Epetra_Vector y1(A.OperatorRangeMap());
417 Epetra_Vector y2(A.OperatorRangeMap());
418 Epetra_Vector y3(A.OperatorRangeMap());
420 A1->Multiply(
false, x, y2);
421 A2->Multiply(
false, x, y3);
423 y1.Norm2(&residual);
double rAx1 = residual;
424 if (verbose) std::cout <<
"Norm of A*x = " << residual << std::endl;
426 y2.Norm2(&residual);
double rAx2 = residual;
427 if (verbose) std::cout <<
"Norm of A1*x = " << residual << std::endl;
428 ierr +=
checkValues(rAx1,rAx2,
"Norm of A1*x", verbose);
430 y3.Norm2(&residual);
double rAx3 = residual;
431 if (verbose) std::cout <<
"Norm of A2*x = " << residual << std::endl;
432 ierr +=
checkValues(rAx1,rAx3,
"Norm of A2*x", verbose);
441 int MyPID = comm.MyPID();
442 int NumProc = comm.NumProc();
445 int ilower = MyPID *
N;
446 int iupper = (MyPID+1)*
N-1;
448 double filePID = (double)MyPID/(
double)100000;
449 std::ostringstream stream;
451 stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
453 std::string fileName(filename);
454 fileName += stream.str().substr(1,7);
456 std::ofstream myfile(fileName.c_str());
458 if(myfile.is_open()){
459 myfile << ilower <<
" " << iupper <<
" " << ilower <<
" " << iupper << std::endl;
460 for(
int i = ilower; i <= iupper; i++){
461 for(
int j=i-5; j <= i+5; j++){
462 if(j >= 0 && j <
N*NumProc)
463 myfile << i <<
" " << j <<
" " << (double)rand()/(double)RAND_MAX << std::endl;
469 std::cout <<
"\nERROR:\nCouldn't open file.\n";
int OperatorToMatrixMarketFile(const char *filename, const Epetra_Operator &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_Operator object to a Matrix Market format file, forming the coefficients by applying...
int MatrixMarketFileToMultiVector(const char *filename, const Epetra_BlockMap &map, Epetra_MultiVector *&A)
Constructs an Epetra_MultiVector object from a Matrix Market format file.
int MultiVectorToMatlabFile(const char *filename, const Epetra_MultiVector &A)
Writes an Epetra_MultiVector object to a file that is compatible with Matlab.
int BlockMapToMatrixMarketFile(const char *filename, const Epetra_BlockMap &map, const char *mapName, const char *mapDescription, bool writeHeader)
Writes an Epetra_BlockMap or Epetra_Map object to a Matrix Market format file.
int runHypreTest(Epetra_CrsMatrix &A)
std::string EpetraExt_Version()
int runTests(Epetra_Map &map, Epetra_CrsMatrix &A, Epetra_Vector &x, Epetra_Vector &b, Epetra_Vector &xexact, bool verbose)
int RowMatrixToMatlabFile(const char *filename, const Epetra_RowMatrix &A)
Writes an Epetra_RowMatrix object to a file that is compatible with Matlab.
int MatrixMarketFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
int MatrixMarketFileToMap(const char *filename, const Epetra_Comm &comm, Epetra_Map *&map)
Constructs an Epetra_BlockMap object from a Matrix Market format file.
int generateHyprePrintOut(const char *filename, const Epetra_Comm &comm)
int MultiVectorToMatrixMarketFile(const char *filename, const Epetra_MultiVector &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_MultiVector object to a Matrix Market format file.
Poisson2dOperator: A sample implementation of the Epetra_Operator class.
int runOperatorTests(Epetra_Operator &A, bool verbose)
int MatlabFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&A)
Constructs an Epetra_CrsMatrix object from a Matlab format file, distributes rows evenly across proce...
int RowMatrixToMatrixMarketFile(const char *filename, const Epetra_RowMatrix &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_RowMatrix object to a Matrix Market format file.
int VectorToMatrixMarketFile(const char *filename, const Epetra_Vector &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_Vector object to a Matrix Market format file.
int OperatorToMatlabFile(const char *filename, const Epetra_Operator &A)
Writes an Epetra_Operator object to a file that is compatible with Matlab.
int MatrixMarketFileToVector(const char *filename, const Epetra_BlockMap &map, Epetra_Vector *&A)
Constructs an Epetra_Vector object from a Matrix Market format file.
int HypreFileToCrsMatrix(const char *filename, const Epetra_Comm &comm, Epetra_CrsMatrix *&Matrix)
Constructs an Epetra_CrsMatrix object from a Hypre Matrix Print command, the row map is specified...
int main(int argc, char *argv[])
int checkValues(double x, double y, std::string message="", bool verbose=false)