EpetraExt Package Browser (Single Doxygen Collection) Development
test/inout/cxx_main.cpp
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Epetra_ConfigDefs.h"
00031 #include "EpetraExt_Version.h"
00032 #ifdef EPETRA_MPI
00033 #include "mpi.h"
00034 #include "Epetra_MpiComm.h"
00035 #else
00036 #include "Epetra_SerialComm.h"
00037 #endif
00038 #include "Trilinos_Util.h"
00039 #include "Epetra_Comm.h"
00040 #include "Epetra_Map.h"
00041 #include "Epetra_Time.h"
00042 #include "Epetra_BlockMap.h"
00043 #include "Epetra_MultiVector.h"
00044 #include "Epetra_Vector.h"
00045 #include "Epetra_Export.h"
00046 
00047 #include "Epetra_VbrMatrix.h"
00048 #include "Epetra_CrsMatrix.h"
00049 #include "EpetraExt_RowMatrixOut.h"
00050 #include "EpetraExt_OperatorOut.h"
00051 #include "EpetraExt_MultiVectorOut.h"
00052 #include "EpetraExt_VectorOut.h"
00053 #include "EpetraExt_BlockMapOut.h"
00054 #include "EpetraExt_BlockMapIn.h"
00055 #include "EpetraExt_CrsMatrixIn.h"
00056 #include "EpetraExt_MultiVectorIn.h"
00057 #include "EpetraExt_VectorIn.h"
00058 // include the header files again to check if guards are in place
00059 #include "EpetraExt_RowMatrixOut.h"
00060 #include "EpetraExt_VectorOut.h"
00061 #include "EpetraExt_BlockMapOut.h"
00062 #include "EpetraExt_BlockMapIn.h"
00063 #include "EpetraExt_CrsMatrixIn.h"
00064 #include "EpetraExt_MultiVectorIn.h"
00065 #include "EpetraExt_VectorIn.h"
00066 #include <string>
00067 #include <vector>
00068 #include "Poisson2dOperator.h"
00069 // prototypes
00070 
00071 int checkValues( double x, double y, string message = "", bool verbose = false) { 
00072   if (fabs((x-y)/x) > 0.01) {
00073     if (verbose) cout << "********** " << message << " check failed.********** " << endl;
00074     return(1); 
00075   }
00076   else {
00077     if (verbose) cout << message << " check OK." << endl;    
00078     return(0);
00079   }
00080 }
00081 int runTests(Epetra_Map & map, Epetra_CrsMatrix & A, Epetra_Vector & x, Epetra_Vector & b, Epetra_Vector & xexact, bool verbose);
00082 int runOperatorTests(Epetra_Operator & A, bool verbose);
00083 int generateHyprePrintOut(const char* filename, const Epetra_Comm &comm);
00084 int runHypreTest(Epetra_CrsMatrix &A);
00085 
00086 int main(int argc, char *argv[]) {
00087 
00088 #ifdef EPETRA_MPI
00089   MPI_Init(&argc,&argv);
00090   Epetra_MpiComm comm (MPI_COMM_WORLD);
00091 #else
00092   Epetra_SerialComm comm;
00093 #endif
00094 
00095   int MyPID = comm.MyPID();
00096 
00097   bool verbose = false;
00098   bool verbose1 = false; 
00099   // Check if we should print results to standard out
00100   if (argc > 1) {
00101     if ((argv[1][0] == '-') && (argv[1][1] == 'v')) {
00102       verbose1 = true;
00103       if (MyPID==0) verbose = true;
00104     }
00105   }
00106   if (verbose)
00107     cout << EpetraExt::EpetraExt_Version() << endl << endl;
00108 
00109   if (verbose1) cout << comm << endl;
00110 
00111 
00112   // Uncomment the next three lines to debug in mpi mode
00113   //int tmp;
00114   //if (MyPID==0) cin >> tmp;
00115   //comm.Barrier();
00116 
00117   Epetra_Map * map;
00118   Epetra_CrsMatrix * A; 
00119   Epetra_Vector * x; 
00120   Epetra_Vector * b;
00121   Epetra_Vector * xexact;
00122 
00123   int nx = 20*comm.NumProc();
00124   int ny = 30;
00125   int npoints = 7;
00126   int xoff[] = {-1,  0,  1, -1,  0,  1,  0};
00127   int yoff[] = {-1, -1, -1,  0,  0,  0,  1};
00128 
00129    
00130   int ierr = 0;
00131   // Call routine to read in HB problem 0-base
00132   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
00133 
00134   ierr += runTests(*map, *A, *x, *b, *xexact, verbose);
00135 
00136   delete A;
00137   delete x;
00138   delete b;
00139   delete xexact;
00140   delete map;
00141 
00142   // Call routine to read in HB problem 1-base
00143   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, 1);
00144 
00145   ierr += runTests(*map, *A, *x, *b, *xexact, verbose);
00146 
00147   delete A;
00148   delete x;
00149   delete b;
00150   delete xexact;
00151   delete map;
00152 
00153   // Call routine to read in HB problem -1-base
00154   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, -1);
00155 
00156   ierr += runTests(*map, *A, *x, *b, *xexact, verbose);
00157 
00158   delete A;
00159   delete x;
00160   delete b;
00161   delete xexact;
00162   delete map;
00163 
00164   int nx1 = 5;
00165   int ny1 = 4;
00166   Poisson2dOperator Op(nx1, ny1, comm);
00167   ierr += runOperatorTests(Op, verbose);
00168 
00169   generateHyprePrintOut("MyMatrixFile", comm);
00170 
00171   EPETRA_CHK_ERR(EpetraExt::HypreFileToCrsMatrix("MyMatrixFile", comm, A));
00172   
00173   runHypreTest(*A);
00174   delete A;
00175 
00176   #ifdef EPETRA_MPI
00177   MPI_Finalize() ;
00178 #endif
00179 
00180   return(ierr);
00181 }
00182 
00183 int runHypreTest(Epetra_CrsMatrix &A){
00184   
00185   Epetra_Vector X(A.RowMap());
00186   EPETRA_CHK_ERR(X.Random());
00187   Epetra_Vector Y(A.RowMap());
00188   EPETRA_CHK_ERR(A.Multiply(false, X, Y));
00189   
00190   return 0;
00191 }
00192 
00193 int runTests(Epetra_Map & map, Epetra_CrsMatrix & A, Epetra_Vector & x, Epetra_Vector & b, Epetra_Vector & xexact, bool verbose) {
00194 
00195   int ierr = 0;
00196 
00197   // Create MultiVectors and put x, b, xexact in both columns of X, B, and Xexact, respectively.
00198   Epetra_MultiVector X( map, 2, false );
00199   Epetra_MultiVector B( map, 2, false );
00200   Epetra_MultiVector Xexact( map, 2, false );
00201 
00202   for (int i=0; i<X.NumVectors(); ++i) {
00203     *X(i) = x;
00204     *B(i) = b;
00205     *Xexact(i) = xexact;
00206   }
00207   double residual;
00208   std::vector<double> residualmv(2);
00209   residual = A.NormInf(); double rAInf = residual;
00210   if (verbose) cout << "Inf Norm of A                                                     = " << residual << endl;
00211   residual = A.NormOne(); double rAOne = residual;
00212   if (verbose) cout << "One Norm of A                                                     = " << residual << endl;
00213   xexact.Norm2(&residual); double rxx = residual; 
00214   Xexact.Norm2(&residualmv[0]); std::vector<double> rXX( residualmv );  
00215   if (verbose) cout << "Norm of xexact                                                    = " << residual << endl;
00216   if (verbose) cout << "Norm of Xexact                                                    = (" << residualmv[0] << ", " <<residualmv[1] <<")"<< endl;
00217   Epetra_Vector tmp1(map);
00218   Epetra_MultiVector tmp1mv(map,2,false);
00219   A.Multiply(false, xexact, tmp1);
00220   A.Multiply(false, Xexact, tmp1mv);
00221   tmp1.Norm2(&residual); double rAx = residual;
00222   tmp1mv.Norm2(&residualmv[0]); std::vector<double> rAX( residualmv );
00223   if (verbose) cout << "Norm of Ax                                                        = " << residual << endl;
00224   if (verbose) cout << "Norm of AX                                                        = (" << residualmv[0] << ", " << residualmv[1] <<")"<< endl;
00225   b.Norm2(&residual); double rb = residual;
00226   B.Norm2(&residualmv[0]); std::vector<double> rB( residualmv );
00227   if (verbose) cout << "Norm of b (should equal norm of Ax)                               = " << residual << endl;
00228   if (verbose) cout << "Norm of B (should equal norm of AX)                               = (" << residualmv[0] << ", " << residualmv[1] <<")"<< endl;
00229   tmp1.Update(1.0, b, -1.0);
00230   tmp1mv.Update(1.0, B, -1.0);
00231   tmp1.Norm2(&residual);
00232   tmp1mv.Norm2(&residualmv[0]);
00233   if (verbose) cout << "Norm of difference between compute Ax and Ax from file            = " << residual << endl;
00234   if (verbose) cout << "Norm of difference between compute AX and AX from file            = (" << residualmv[0] << ", " << residualmv[1] <<")"<< endl;
00235   map.Comm().Barrier();
00236 
00237   EPETRA_CHK_ERR(EpetraExt::BlockMapToMatrixMarketFile("Test_map.mm", map, "Official EpetraExt test map", 
00238                    "This is the official EpetraExt test map generated by the EpetraExt regression tests"));
00239 
00240   EPETRA_CHK_ERR(EpetraExt::RowMatrixToMatrixMarketFile("Test_A.mm", A, "Official EpetraExt test matrix", 
00241               "This is the official EpetraExt test matrix generated by the EpetraExt regression tests"));
00242 
00243   EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_x.mm", x, "Official EpetraExt test initial guess", 
00244                  "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));
00245 
00246   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvX.mm", X, "Official EpetraExt test initial guess", 
00247                             "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));
00248                
00249   EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_xexact.mm", xexact, "Official EpetraExt test exact solution", 
00250                  "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));
00251 
00252   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvXexact.mm", Xexact, "Official EpetraExt test exact solution", 
00253                       "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));
00254                
00255   EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_b.mm", b, "Official EpetraExt test right hand side", 
00256                  "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));
00257 
00258   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvB.mm", B, "Official EpetraExt test right hand side", 
00259                       "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));
00260                
00261   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatlabFile("Test_mvB.mat", B));
00262  
00263   EPETRA_CHK_ERR(EpetraExt::RowMatrixToMatlabFile("Test_A.dat", A));
00264 
00265   Epetra_Map * map1;
00266   Epetra_CrsMatrix * A1; 
00267   Epetra_CrsMatrix * A2; 
00268   Epetra_CrsMatrix * A3; 
00269   Epetra_Vector * x1; 
00270   Epetra_Vector * b1;
00271   Epetra_Vector * xexact1;
00272   Epetra_MultiVector * X1; 
00273   Epetra_MultiVector * B1;
00274   Epetra_MultiVector * Xexact1;
00275 
00276   EpetraExt::MatrixMarketFileToMap("Test_map.mm", map.Comm(), map1);
00277 
00278   if (map.SameAs(*map1)) {
00279     if (verbose) cout << "Maps are equal.  In/Out works." << endl;
00280   }
00281   else {
00282     if (verbose) cout << "Maps are not equal.  In/Out fails." << endl;
00283     ierr += 1;
00284   }
00285   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", *map1, A1));
00286   // If map is zero-based, then we can compare to the convenient reading versions
00287   if (map1->IndexBase()==0) EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", map1->Comm(), A2));
00288   if (map1->IndexBase()==0) EPETRA_CHK_ERR(EpetraExt::MatlabFileToCrsMatrix("Test_A.dat", map1->Comm(), A3));
00289   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_x.mm", *map1, x1));
00290   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_xexact.mm", *map1, xexact1));
00291   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_b.mm", *map1, b1));
00292   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToMultiVector("Test_mvX.mm", *map1, X1));
00293   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToMultiVector("Test_mvXexact.mm", *map1, Xexact1));
00294   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToMultiVector("Test_mvB.mm", *map1, B1));
00295 
00296   residual = A1->NormInf(); double rA1Inf = residual;
00297   if (verbose) cout << "Inf Norm of A1                                                    = " << residual << endl;
00298   ierr += checkValues(rA1Inf,rAInf,"Inf Norm of A", verbose);
00299 
00300   residual = A1->NormOne(); double rA1One = residual;
00301   if (verbose) cout << "One Norm of A1                                                    = " << residual << endl;
00302   ierr += checkValues(rA1One,rAOne,"One Norm of A", verbose);
00303 
00304   xexact1->Norm2(&residual); double rxx1 = residual;
00305   if (verbose) cout << "Norm of xexact1                                                   = " << residual << endl;
00306   ierr += checkValues(rxx1,rxx,"Norm of xexact", verbose);
00307 
00308   Xexact1->Norm2(&residualmv[0]); std::vector<double> rXX1(residualmv);
00309   if (verbose) cout << "Norm of Xexact1                                                   = (" << residualmv[0] <<", " <<residualmv[1]<<")"<< endl;
00310   ierr += checkValues(rXX1[0],rXX[0],"Norm of Xexact", verbose);
00311   ierr += checkValues(rXX1[1],rXX[1],"Norm of Xexact", verbose);
00312 
00313   Epetra_Vector tmp11(*map1);
00314   A1->Multiply(false, *xexact1, tmp11);
00315 
00316   Epetra_MultiVector tmp11mv(*map1,2,false);
00317   A1->Multiply(false, *Xexact1, tmp11mv);
00318 
00319   tmp11.Norm2(&residual); double rAx1 = residual;
00320   if (verbose) cout << "Norm of A1*x1                                                     = " << residual << endl;
00321   ierr += checkValues(rAx1,rAx,"Norm of A1*x", verbose);
00322 
00323   tmp11mv.Norm2(&residualmv[0]); std::vector<double> rAX1(residualmv);
00324   if (verbose) cout << "Norm of A1*X1                                                     = (" << residualmv[0] <<", "<<residualmv[1]<<")"<< endl;
00325   ierr += checkValues(rAX1[0],rAX[0],"Norm of A1*X", verbose);
00326   ierr += checkValues(rAX1[1],rAX[1],"Norm of A1*X", verbose);
00327 
00328   if (map1->IndexBase()==0) {
00329     Epetra_Vector tmp12(*map1);
00330     A2->Multiply(false, *xexact1, tmp12);
00331     
00332     tmp12.Norm2(&residual); double rAx2 = residual;
00333     if (verbose) cout << "Norm of A2*x1                                                     = " << residual << endl;
00334     ierr += checkValues(rAx2,rAx,"Norm of A2*x", verbose);
00335 
00336     Epetra_Vector tmp13(*map1);
00337     A3->Multiply(false, *xexact1, tmp13);
00338     
00339     tmp13.Norm2(&residual); double rAx3 = residual;
00340     if (verbose) cout << "Norm of A3*x1                                                     = " << residual << endl;
00341     ierr += checkValues(rAx3,rAx,"Norm of A3*x", verbose);
00342   }
00343   b1->Norm2(&residual); double rb1 = residual;
00344   if (verbose) cout << "Norm of b1 (should equal norm of Ax)                              = " << residual << endl;
00345   ierr += checkValues(rb1,rb,"Norm of b", verbose);
00346 
00347   B1->Norm2(&residualmv[0]); std::vector<double> rB1(residualmv);
00348   if (verbose) cout << "Norm of B1 (should equal norm of AX)                              = (" << residualmv[0] <<", "<<residualmv[1]<<")"<< endl;
00349   ierr += checkValues(rB1[0],rB[0],"Norm of B", verbose);
00350   ierr += checkValues(rB1[1],rB[1],"Norm of B", verbose);
00351 
00352   tmp11.Update(1.0, *b1, -1.0);
00353   tmp11.Norm2(&residual);
00354   if (verbose) cout << "Norm of difference between computed A1x1 and A1x1 from file        = " << residual << endl;
00355   ierr += checkValues(residual,0.0,"Norm of difference between computed A1x1 and A1x1 from file", verbose);
00356 
00357   tmp11mv.Update(1.0, *B1, -1.0);
00358   tmp11mv.Norm2(&residualmv[0]);
00359   if (verbose) cout << "Norm of difference between computed A1X1 and A1X1 from file        = (" << residualmv[0] << ", "<<residualmv[1]<<")"<< endl;
00360   ierr += checkValues(residualmv[0],0.0,"Norm of difference between computed A1X1 and A1X1 from file", verbose);
00361   ierr += checkValues(residualmv[1],0.0,"Norm of difference between computed A1X1 and A1X1 from file", verbose);
00362 
00363   if (map1->IndexBase()==0) {delete A2; delete A3;}
00364   delete A1;
00365   delete x1;
00366   delete b1;
00367   delete xexact1;
00368   delete X1;
00369   delete B1;
00370   delete Xexact1;
00371   delete map1;
00372 
00373   return(ierr);
00374 }
00375 
00376 int runOperatorTests(Epetra_Operator & A, bool verbose) {
00377 
00378   int ierr = 0;
00379 
00380 
00381   double residual;
00382   EPETRA_CHK_ERR(EpetraExt::OperatorToMatrixMarketFile("Test_A1.mm", A, "Official EpetraExt test operator", 
00383               "This is the official EpetraExt test operator generated by the EpetraExt regression tests"));
00384   EPETRA_CHK_ERR(EpetraExt::OperatorToMatlabFile("Test_A1.dat", A));
00385 
00386   A.OperatorRangeMap().Comm().Barrier();
00387   A.OperatorRangeMap().Comm().Barrier();
00388   Epetra_CrsMatrix * A1; 
00389   Epetra_CrsMatrix * A2; 
00390   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A1.mm", A.OperatorRangeMap(), A1));
00391   EPETRA_CHK_ERR(EpetraExt::MatlabFileToCrsMatrix("Test_A1.dat", A.OperatorRangeMap().Comm(), A2));
00392 
00393 
00394   residual = A.NormInf(); double rAInf = residual;
00395   if (verbose) cout << "Inf Norm of Operator A                                            = " << residual << endl;
00396   residual = A1->NormInf(); double rA1Inf = residual;
00397   if (verbose) cout << "Inf Norm of Matrix A1                                             = " << residual << endl;
00398   ierr += checkValues(rA1Inf,rAInf,"Inf Norm of A", verbose);
00399 
00400 
00401   Epetra_Vector x(A.OperatorDomainMap()); x.Random();
00402   Epetra_Vector y1(A.OperatorRangeMap());
00403   Epetra_Vector y2(A.OperatorRangeMap());
00404   Epetra_Vector y3(A.OperatorRangeMap());
00405   A.Apply(x,y1);
00406   A1->Multiply(false, x, y2);
00407   A2->Multiply(false, x, y3);
00408 
00409   y1.Norm2(&residual); double rAx1 = residual;
00410   if (verbose) cout << "Norm of A*x                                                       = " << residual << endl;
00411 
00412   y2.Norm2(&residual); double rAx2 = residual;
00413   if (verbose) cout << "Norm of A1*x                                                      = " << residual << endl;
00414   ierr += checkValues(rAx1,rAx2,"Norm of A1*x", verbose);
00415 
00416   y3.Norm2(&residual); double rAx3 = residual;
00417   if (verbose) cout << "Norm of A2*x                                                      = " << residual << endl;
00418   ierr += checkValues(rAx1,rAx3,"Norm of A2*x", verbose);
00419 
00420   delete A1;
00421   delete A2;
00422 
00423   return(ierr);
00424 }
00425 
00426 int generateHyprePrintOut(const char *filename, const Epetra_Comm &comm){
00427   int MyPID = comm.MyPID();
00428   int NumProc = comm.NumProc();
00429 
00430   int N = 100;
00431   int ilower = MyPID * N;
00432   int iupper = (MyPID+1)*N-1;
00433 
00434   double filePID = (double)MyPID/(double)100000;
00435   std::ostringstream stream;
00436   // Using setprecision() puts it in the string
00437   stream << std::setiosflags(std::ios::fixed) << std::setprecision(5) << filePID;
00438   // Then just ignore the first character
00439   std::string fileName(filename);
00440   fileName += stream.str().substr(1,7);
00441 
00442   std::ofstream myfile(fileName.c_str());
00443 
00444   if(myfile.is_open()){
00445     myfile << ilower << " " << iupper << " " << ilower << " " << iupper << endl;
00446     for(int i = ilower; i <= iupper; i++){
00447       for(int j=i-5; j <= i+5; j++){
00448         if(j >= 0 && j < N*NumProc)
00449           myfile << i << " " << j << " " << (double)rand()/(double)RAND_MAX << endl;
00450       }
00451     }
00452     myfile.close();
00453     return 0;
00454   } else {
00455     cout << "\nERROR:\nCouldn't open file.\n";
00456     return -1;
00457   }
00458 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines