EpetraExt Package Browser (Single Doxygen Collection): cxx

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 
00084 int main(int argc, char *argv[]) {
00085 
00086 #ifdef EPETRA_MPI
00087   MPI_Init(&argc,&argv);
00088   Epetra_MpiComm comm (MPI_COMM_WORLD);
00089 #else
00090   Epetra_SerialComm comm;
00091 #endif
00092 
00093   int MyPID = comm.MyPID();
00094 
00095   bool verbose = false;
00096   bool verbose1 = false; 
00097   // Check if we should print results to standard out
00098   if (argc > 1) {
00099     if ((argv[1][0] == '-') && (argv[1][1] == 'v')) {
00100       verbose1 = true;
00101       if (MyPID==0) verbose = true;
00102     }
00103   }
00104   if (verbose)
00105     cout << EpetraExt::EpetraExt_Version() << endl << endl;
00106 
00107   if (verbose1) cout << comm << endl;
00108 
00109 
00110   // Uncomment the next three lines to debug in mpi mode
00111   //int tmp;
00112   //if (MyPID==0) cin >> tmp;
00113   //comm.Barrier();
00114 
00115   Epetra_Map * map;
00116   Epetra_CrsMatrix * A; 
00117   Epetra_Vector * x; 
00118   Epetra_Vector * b;
00119   Epetra_Vector * xexact;
00120 
00121   int nx = 20*comm.NumProc();
00122   int ny = 30;
00123   int npoints = 7;
00124   int xoff[] = {-1,  0,  1, -1,  0,  1,  0};
00125   int yoff[] = {-1, -1, -1,  0,  0,  0,  1};
00126 
00127    
00128   int ierr = 0;
00129   // Call routine to read in HB problem 0-base
00130   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
00131 
00132   ierr += runTests(*map, *A, *x, *b, *xexact, verbose);
00133 
00134   delete A;
00135   delete x;
00136   delete b;
00137   delete xexact;
00138   delete map;
00139 
00140   // Call routine to read in HB problem 1-base
00141   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, 1);
00142 
00143   ierr += runTests(*map, *A, *x, *b, *xexact, verbose);
00144 
00145   delete A;
00146   delete x;
00147   delete b;
00148   delete xexact;
00149   delete map;
00150 
00151   // Call routine to read in HB problem -1-base
00152   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, -1);
00153 
00154   ierr += runTests(*map, *A, *x, *b, *xexact, verbose);
00155 
00156   delete A;
00157   delete x;
00158   delete b;
00159   delete xexact;
00160   delete map;
00161 
00162   int nx1 = 5;
00163   int ny1 = 4;
00164   Poisson2dOperator Op(nx1, ny1, comm);
00165   ierr += runOperatorTests(Op, verbose);
00166 
00167   #ifdef EPETRA_MPI
00168   MPI_Finalize() ;
00169 #endif
00170 
00171   return(ierr);
00172 }
00173 
00174 int runTests(Epetra_Map & map, Epetra_CrsMatrix & A, Epetra_Vector & x, Epetra_Vector & b, Epetra_Vector & xexact, bool verbose) {
00175 
00176   int ierr = 0;
00177 
00178   // Create MultiVectors and put x, b, xexact in both columns of X, B, and Xexact, respectively.
00179   Epetra_MultiVector X( map, 2, false );
00180   Epetra_MultiVector B( map, 2, false );
00181   Epetra_MultiVector Xexact( map, 2, false );
00182 
00183   for (int i=0; i<X.NumVectors(); ++i) {
00184     *X(i) = x;
00185     *B(i) = b;
00186     *Xexact(i) = xexact;
00187   }
00188   double residual;
00189   std::vector<double> residualmv(2);
00190   residual = A.NormInf(); double rAInf = residual;
00191   if (verbose) cout << "Inf Norm of A                                                     = " << residual << endl;
00192   residual = A.NormOne(); double rAOne = residual;
00193   if (verbose) cout << "One Norm of A                                                     = " << residual << endl;
00194   xexact.Norm2(&residual); double rxx = residual; 
00195   Xexact.Norm2(&residualmv[0]); std::vector<double> rXX( residualmv );  
00196   if (verbose) cout << "Norm of xexact                                                    = " << residual << endl;
00197   if (verbose) cout << "Norm of Xexact                                                    = (" << residualmv[0] << ", " <<residualmv[1] <<")"<< endl;
00198   Epetra_Vector tmp1(map);
00199   Epetra_MultiVector tmp1mv(map,2,false);
00200   A.Multiply(false, xexact, tmp1);
00201   A.Multiply(false, Xexact, tmp1mv);
00202   tmp1.Norm2(&residual); double rAx = residual;
00203   tmp1mv.Norm2(&residualmv[0]); std::vector<double> rAX( residualmv );
00204   if (verbose) cout << "Norm of Ax                                                        = " << residual << endl;
00205   if (verbose) cout << "Norm of AX                                                        = (" << residualmv[0] << ", " << residualmv[1] <<")"<< endl;
00206   b.Norm2(&residual); double rb = residual;
00207   B.Norm2(&residualmv[0]); std::vector<double> rB( residualmv );
00208   if (verbose) cout << "Norm of b (should equal norm of Ax)                               = " << residual << endl;
00209   if (verbose) cout << "Norm of B (should equal norm of AX)                               = (" << residualmv[0] << ", " << residualmv[1] <<")"<< endl;
00210   tmp1.Update(1.0, b, -1.0);
00211   tmp1mv.Update(1.0, B, -1.0);
00212   tmp1.Norm2(&residual);
00213   tmp1mv.Norm2(&residualmv[0]);
00214   if (verbose) cout << "Norm of difference between compute Ax and Ax from file            = " << residual << endl;
00215   if (verbose) cout << "Norm of difference between compute AX and AX from file            = (" << residualmv[0] << ", " << residualmv[1] <<")"<< endl;
00216   map.Comm().Barrier();
00217 
00218   EPETRA_CHK_ERR(EpetraExt::BlockMapToMatrixMarketFile("Test_map.mm", map, "Official EpetraExt test map", 
00219                    "This is the official EpetraExt test map generated by the EpetraExt regression tests"));
00220 
00221   EPETRA_CHK_ERR(EpetraExt::RowMatrixToMatrixMarketFile("Test_A.mm", A, "Official EpetraExt test matrix", 
00222               "This is the official EpetraExt test matrix generated by the EpetraExt regression tests"));
00223 
00224   EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_x.mm", x, "Official EpetraExt test initial guess", 
00225                  "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));
00226 
00227   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvX.mm", X, "Official EpetraExt test initial guess", 
00228                             "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));
00229                
00230   EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_xexact.mm", xexact, "Official EpetraExt test exact solution", 
00231                  "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));
00232 
00233   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvXexact.mm", Xexact, "Official EpetraExt test exact solution", 
00234                       "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));
00235                
00236   EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_b.mm", b, "Official EpetraExt test right hand side", 
00237                  "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));
00238 
00239   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvB.mm", B, "Official EpetraExt test right hand side", 
00240                       "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));
00241                
00242   EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatlabFile("Test_mvB.mat", B));
00243  
00244   EPETRA_CHK_ERR(EpetraExt::RowMatrixToMatlabFile("Test_A.dat", A));
00245 
00246   Epetra_Map * map1;
00247   Epetra_CrsMatrix * A1; 
00248   Epetra_CrsMatrix * A2; 
00249   Epetra_CrsMatrix * A3; 
00250   Epetra_Vector * x1; 
00251   Epetra_Vector * b1;
00252   Epetra_Vector * xexact1;
00253   Epetra_MultiVector * X1; 
00254   Epetra_MultiVector * B1;
00255   Epetra_MultiVector * Xexact1;
00256 
00257   EpetraExt::MatrixMarketFileToMap("Test_map.mm", map.Comm(), map1);
00258 
00259   if (map.SameAs(*map1)) {
00260     if (verbose) cout << "Maps are equal.  In/Out works." << endl;
00261   }
00262   else {
00263     if (verbose) cout << "Maps are not equal.  In/Out fails." << endl;
00264     ierr += 1;
00265   }
00266   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", *map1, A1));
00267   // If map is zero-based, then we can compare to the convenient reading versions
00268   if (map1->IndexBase()==0) EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", map1->Comm(), A2));
00269   if (map1->IndexBase()==0) EPETRA_CHK_ERR(EpetraExt::MatlabFileToCrsMatrix("Test_A.dat", map1->Comm(), A3));
00270   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_x.mm", *map1, x1));
00271   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_xexact.mm", *map1, xexact1));
00272   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_b.mm", *map1, b1));
00273   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToMultiVector("Test_mvX.mm", *map1, X1));
00274   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToMultiVector("Test_mvXexact.mm", *map1, Xexact1));
00275   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToMultiVector("Test_mvB.mm", *map1, B1));
00276 
00277 
00278   residual = A1->NormInf(); double rA1Inf = residual;
00279   if (verbose) cout << "Inf Norm of A1                                                    = " << residual << endl;
00280   ierr += checkValues(rA1Inf,rAInf,"Inf Norm of A", verbose);
00281 
00282   residual = A1->NormOne(); double rA1One = residual;
00283   if (verbose) cout << "One Norm of A1                                                    = " << residual << endl;
00284   ierr += checkValues(rA1One,rAOne,"One Norm of A", verbose);
00285 
00286   xexact1->Norm2(&residual); double rxx1 = residual;
00287   if (verbose) cout << "Norm of xexact1                                                   = " << residual << endl;
00288   ierr += checkValues(rxx1,rxx,"Norm of xexact", verbose);
00289 
00290   Xexact1->Norm2(&residualmv[0]); std::vector<double> rXX1(residualmv);
00291   if (verbose) cout << "Norm of Xexact1                                                   = (" << residualmv[0] <<", " <<residualmv[1]<<")"<< endl;
00292   ierr += checkValues(rXX1[0],rXX[0],"Norm of Xexact", verbose);
00293   ierr += checkValues(rXX1[1],rXX[1],"Norm of Xexact", verbose);
00294 
00295   Epetra_Vector tmp11(*map1);
00296   A1->Multiply(false, *xexact1, tmp11);
00297 
00298   Epetra_MultiVector tmp11mv(*map1,2,false);
00299   A1->Multiply(false, *Xexact1, tmp11mv);
00300 
00301   tmp11.Norm2(&residual); double rAx1 = residual;
00302   if (verbose) cout << "Norm of A1*x1                                                     = " << residual << endl;
00303   ierr += checkValues(rAx1,rAx,"Norm of A1*x", verbose);
00304 
00305   tmp11mv.Norm2(&residualmv[0]); std::vector<double> rAX1(residualmv);
00306   if (verbose) cout << "Norm of A1*X1                                                     = (" << residualmv[0] <<", "<<residualmv[1]<<")"<< endl;
00307   ierr += checkValues(rAX1[0],rAX[0],"Norm of A1*X", verbose);
00308   ierr += checkValues(rAX1[1],rAX[1],"Norm of A1*X", verbose);
00309 
00310   if (map1->IndexBase()==0) {
00311     Epetra_Vector tmp12(*map1);
00312     A2->Multiply(false, *xexact1, tmp12);
00313     
00314     tmp12.Norm2(&residual); double rAx2 = residual;
00315     if (verbose) cout << "Norm of A2*x1                                                     = " << residual << endl;
00316     ierr += checkValues(rAx2,rAx,"Norm of A2*x", verbose);
00317 
00318     Epetra_Vector tmp13(*map1);
00319     A3->Multiply(false, *xexact1, tmp13);
00320     
00321     tmp13.Norm2(&residual); double rAx3 = residual;
00322     if (verbose) cout << "Norm of A3*x1                                                     = " << residual << endl;
00323     ierr += checkValues(rAx3,rAx,"Norm of A3*x", verbose);
00324   }
00325   b1->Norm2(&residual); double rb1 = residual;
00326   if (verbose) cout << "Norm of b1 (should equal norm of Ax)                              = " << residual << endl;
00327   ierr += checkValues(rb1,rb,"Norm of b", verbose);
00328 
00329   B1->Norm2(&residualmv[0]); std::vector<double> rB1(residualmv);
00330   if (verbose) cout << "Norm of B1 (should equal norm of AX)                              = (" << residualmv[0] <<", "<<residualmv[1]<<")"<< endl;
00331   ierr += checkValues(rB1[0],rB[0],"Norm of B", verbose);
00332   ierr += checkValues(rB1[1],rB[1],"Norm of B", verbose);
00333 
00334   tmp11.Update(1.0, *b1, -1.0);
00335   tmp11.Norm2(&residual);
00336   if (verbose) cout << "Norm of difference between computed A1x1 and A1x1 from file        = " << residual << endl;
00337   ierr += checkValues(residual,0.0,"Norm of difference between computed A1x1 and A1x1 from file", verbose);
00338 
00339   tmp11mv.Update(1.0, *B1, -1.0);
00340   tmp11mv.Norm2(&residualmv[0]);
00341   if (verbose) cout << "Norm of difference between computed A1X1 and A1X1 from file        = (" << residualmv[0] << ", "<<residualmv[1]<<")"<< endl;
00342   ierr += checkValues(residualmv[0],0.0,"Norm of difference between computed A1X1 and A1X1 from file", verbose);
00343   ierr += checkValues(residualmv[1],0.0,"Norm of difference between computed A1X1 and A1X1 from file", verbose);
00344 
00345   if (map1->IndexBase()==0) {delete A2; delete A3;}
00346   delete A1;
00347   delete x1;
00348   delete b1;
00349   delete xexact1;
00350   delete X1;
00351   delete B1;
00352   delete Xexact1;
00353   delete map1;
00354 
00355   return(ierr);
00356 }
00357 
00358 int runOperatorTests(Epetra_Operator & A, bool verbose) {
00359 
00360   int ierr = 0;
00361 
00362 
00363   double residual;
00364   EPETRA_CHK_ERR(EpetraExt::OperatorToMatrixMarketFile("Test_A1.mm", A, "Official EpetraExt test operator", 
00365               "This is the official EpetraExt test operator generated by the EpetraExt regression tests"));
00366   EPETRA_CHK_ERR(EpetraExt::OperatorToMatlabFile("Test_A1.dat", A));
00367 
00368   A.OperatorRangeMap().Comm().Barrier();
00369   A.OperatorRangeMap().Comm().Barrier();
00370   Epetra_CrsMatrix * A1; 
00371   Epetra_CrsMatrix * A2; 
00372   EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A1.mm", A.OperatorRangeMap(), A1));
00373   EPETRA_CHK_ERR(EpetraExt::MatlabFileToCrsMatrix("Test_A1.dat", A.OperatorRangeMap().Comm(), A2));
00374 
00375 
00376   residual = A.NormInf(); double rAInf = residual;
00377   if (verbose) cout << "Inf Norm of Operator A                                            = " << residual << endl;
00378   residual = A1->NormInf(); double rA1Inf = residual;
00379   if (verbose) cout << "Inf Norm of Matrix A1                                             = " << residual << endl;
00380   ierr += checkValues(rA1Inf,rAInf,"Inf Norm of A", verbose);
00381 
00382 
00383   Epetra_Vector x(A.OperatorDomainMap()); x.Random();
00384   Epetra_Vector y1(A.OperatorRangeMap());
00385   Epetra_Vector y2(A.OperatorRangeMap());
00386   Epetra_Vector y3(A.OperatorRangeMap());
00387   A.Apply(x,y1);
00388   A1->Multiply(false, x, y2);
00389   A2->Multiply(false, x, y3);
00390 
00391   y1.Norm2(&residual); double rAx1 = residual;
00392   if (verbose) cout << "Norm of A*x                                                       = " << residual << endl;
00393 
00394   y2.Norm2(&residual); double rAx2 = residual;
00395   if (verbose) cout << "Norm of A1*x                                                      = " << residual << endl;
00396   ierr += checkValues(rAx1,rAx2,"Norm of A1*x", verbose);
00397 
00398   y3.Norm2(&residual); double rAx3 = residual;
00399   if (verbose) cout << "Norm of A2*x                                                      = " << residual << endl;
00400   ierr += checkValues(rAx1,rAx3,"Norm of A2*x", verbose);
00401 
00402   delete A1;
00403   delete A2;
00404 
00405   return(ierr);
00406 }
test/inout/cxx_main.cpp
