EpetraExt Package Browser (Single Doxygen Collection) Development
hypre_UnitTest.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2009) 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 #include "Teuchos_UnitTestHarness.hpp"
00030 #include "HYPRE_IJ_mv.h"
00031 #include "EpetraExt_HypreIJMatrix.h"
00032 #include "EpetraExt_MatrixMatrix.h"
00033 #include "Epetra_ConfigDefs.h"
00034 #include "Epetra_Vector.h"
00035 #include "Epetra_RowMatrix.h"
00036 #include "Epetra_MultiVector.h"
00037 #include "Epetra_CrsMatrix.h"
00038 #include "Epetra_Map.h"
00039 #ifdef HAVE_MPI
00040 #include "mpi.h"
00041 #include "Epetra_MpiComm.h"
00042 #else
00043 #include "Epetra_SerialComm.h"
00044 #endif
00045 
00046 #include "hypre_Helpers.hpp"
00047 #include "Teuchos_ParameterList.hpp"
00048 #include "Teuchos_ParameterEntry.hpp"
00049 #include "Teuchos_ParameterListExceptions.hpp"
00050 #include "Teuchos_Array.hpp"
00051 #include <string>
00052 #include <stdio.h>
00053 #include <map>
00054 
00055 namespace EpetraExt {
00056 
00057 const int N = 100;
00058 const int MatType = 4; //0 -> Unit diagonal, 1 -> Random diagonal, 2 -> Dense, val=col, 3 -> Random Dense, 4 -> Random Sparse
00059 const double tol = 1E-6;
00060 
00061 TEUCHOS_UNIT_TEST( EpetraExt_hypre, Construct ) {
00062 
00063   int ierr = 0;
00064   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 0);
00065   
00066   TEST_EQUALITY(Matrix->Filled(), true);
00067   
00068   for(int i = 0; i < Matrix->NumMyRows(); i++){
00069     int entries;
00070     ierr += Matrix->NumMyRowEntries(i, entries);
00071     TEST_EQUALITY(entries, 1);
00072     int numentries;
00073     Teuchos::Array<double> Values; Values.resize(entries);
00074     Teuchos::Array<int> Indices; Indices.resize(entries);
00075     ierr += Matrix->ExtractMyRowCopy(i, entries, numentries, &Values[0], &Indices[0]);
00076     TEST_EQUALITY(ierr, 0);
00077     TEST_EQUALITY(numentries,1);
00078     for(int j = 0; j < numentries; j++){
00079       TEST_FLOATING_EQUALITY(Values[j],1.0,tol);
00080       TEST_EQUALITY(Indices[j],i);
00081     }
00082   }
00083   delete Matrix;
00084 }
00085 
00086 TEUCHOS_UNIT_TEST( EpetraExt_hypre, MatVec ) {
00087   int ierr = 0;
00088   int num_vectors = 5;
00089   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 0);
00090   
00091   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00092   ierr += X.Random();
00093   Epetra_MultiVector Y(Matrix->RowMatrixRowMap(), num_vectors, true);
00094   
00095   ierr += Matrix->Multiply(false, X, Y);
00096   
00097   TEST_EQUALITY(EquivalentVectors(X,Y,tol),true);
00098   TEST_EQUALITY(ierr, 0);
00099   delete Matrix;
00100 }
00101 
00102 TEUCHOS_UNIT_TEST( EpetraExt_hypre, BetterMatVec ) {
00103   int ierr = 0;
00104   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00105   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00106   //TestMat->Print(std::cout);
00107   int num_vectors = 5;
00108   
00109   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00110   ierr += X.Random();
00111   Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
00112   Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
00113   
00114   ierr += Matrix->Multiply(false, X, Y1);
00115   ierr += TestMat->Multiply(false, X, Y2);
00116   
00117   TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol),true);
00118 
00119   ierr += Matrix->Multiply(false, Y1, X);
00120   ierr += TestMat->Multiply(false, Y1, Y2);
00121 
00122   TEST_EQUALITY(EquivalentVectors(X,Y2,tol),true);
00123 
00124   ierr += Matrix->Multiply(false, Y2, X);
00125   ierr += TestMat->Multiply(false, Y2, Y1);
00126 
00127   TEST_EQUALITY(EquivalentVectors(X,Y1,tol),true);
00128   TEST_EQUALITY_CONST(ierr, 0);
00129   delete Matrix;
00130   delete TestMat;
00131 }
00132 
00133 TEUCHOS_UNIT_TEST( EpetraExt_hypre, TransposeMatVec ) {
00134   int ierr = 0;
00135   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00136   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00137   
00138   int num_vectors = 5;
00139   
00140   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00141   ierr += X.Random();
00142   Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
00143   Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
00144   
00145   ierr += Matrix->Multiply(true, X, Y1);
00146   ierr += TestMat->Multiply(true, X, Y2);
00147   
00148   TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol),true);
00149   TEST_EQUALITY(ierr, 0);
00150   
00151   delete Matrix;
00152   delete TestMat;
00153 }
00154 
00155 TEUCHOS_UNIT_TEST( EpetraExt_hypre, LeftScale ) {
00156   int ierr = 0;
00157   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00158   //EpetraExt_HypreIJMatrix* BackUp = EpetraExt_HypreIJMatrix(Matrix);
00159   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00160   
00161   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00162   ierr += X.Random();
00163   Matrix->NumMyNonzeros();
00164   ierr += Matrix->LeftScale(X);
00165   ierr += TestMat->LeftScale(X);
00166   TEST_EQUALITY(EquivalentMatrices(*Matrix, *TestMat,tol), true);
00167   TEST_EQUALITY(ierr, 0);
00168   delete Matrix;
00169   delete TestMat;
00170 }
00171 
00172 TEUCHOS_UNIT_TEST( EpetraExt_hypre, RightScale ) {
00173   int ierr = 0;
00174   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00175   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00176   
00177   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00178   ierr += X.Random();
00179   Matrix->NumMyNonzeros();
00180   ierr += Matrix->RightScale(X);
00181   ierr += TestMat->RightScale(X);
00182   
00183   TEST_EQUALITY(EquivalentMatrices(*Matrix,*TestMat,tol), true);
00184   TEST_EQUALITY(ierr, 0);
00185   
00186   delete Matrix;
00187   delete TestMat;
00188 }
00189 
00190 TEUCHOS_UNIT_TEST( EpetraExt_hypre, ExtractDiagonalCopy ) {
00191   int ierr = 0;
00192   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00193   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00194   
00195   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00196   Epetra_Vector Y(TestMat->RowMatrixRowMap(),true);
00197   
00198   ierr += Matrix->ExtractDiagonalCopy(X);
00199   ierr += TestMat->ExtractDiagonalCopy(Y);
00200   
00201   TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
00202   TEST_EQUALITY(ierr, 0);
00203   
00204   delete Matrix;
00205   delete TestMat;
00206 }
00207 
00208 TEUCHOS_UNIT_TEST( EpetraExt_hypre, InvRowSums ) {
00209   int ierr = 0;
00210   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00211   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00212   
00213   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00214   Epetra_Vector Y(TestMat->RowMatrixRowMap(),true);
00215   
00216   ierr += Matrix->InvRowSums(X);
00217   ierr += TestMat->InvRowSums(Y);
00218   
00219   TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
00220   TEST_EQUALITY(ierr, 0);
00221   
00222   delete Matrix;
00223   delete TestMat;
00224 }
00225 
00226 TEUCHOS_UNIT_TEST( EpetraExt_hypre, InvColSums ) {
00227   int ierr = 0;
00228   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00229   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00230   
00231   Epetra_Vector X(Matrix->RowMatrixColMap(), true);
00232   Epetra_Vector Y(TestMat->RowMatrixColMap(),true);
00233   
00234   ierr += Matrix->InvColSums(X);
00235   ierr += TestMat->InvColSums(Y);
00236   
00237   TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
00238   TEST_EQUALITY(ierr, 0);
00239 
00240   delete Matrix;
00241   delete TestMat;
00242 }
00243 
00244 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NormInf ) {
00245   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00246   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00247   
00248   double norm1 = Matrix->NormInf();
00249   double norm2 = TestMat->NormInf();
00250   
00251   TEST_FLOATING_EQUALITY(norm1, norm2, tol);
00252   delete Matrix;
00253   delete TestMat;
00254 }
00255 
00256 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NormOne ) {
00257   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00258   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00259   
00260   double norm1 = Matrix->NormOne();
00261   double norm2 = TestMat->NormOne();
00262   
00263   TEST_FLOATING_EQUALITY(norm1, norm2, tol);
00264 
00265   delete Matrix;
00266   delete TestMat;
00267 }
00268 
00269 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalNonzeros ) {
00270   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00271   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00272   
00273   int nnz1 = Matrix->NumGlobalNonzeros();
00274   int nnz2 = TestMat->NumGlobalNonzeros();
00275   
00276   TEST_EQUALITY(nnz1, nnz2);
00277 
00278   delete Matrix;
00279   delete TestMat;
00280 }
00281 
00282 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalRows ) {
00283   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00284   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00285   
00286   int rows1 = Matrix->NumGlobalRows();
00287   int rows2 = TestMat->NumGlobalRows();
00288   
00289   TEST_EQUALITY(rows1, rows2);
00290 
00291   delete Matrix;
00292   delete TestMat;
00293 }
00294 
00295 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalCols ) {
00296   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00297   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00298   
00299   int cols1 = Matrix->NumGlobalCols();
00300   int cols2 = TestMat->NumGlobalCols();
00301   
00302   TEST_EQUALITY(cols1, cols2);
00303 
00304   delete Matrix;
00305   delete TestMat;
00306 }
00307 
00308 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalDiagonals ) {
00309   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00310   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00311   
00312   int hdiag1 = Matrix->NumGlobalDiagonals();
00313   int Ediag2 = TestMat->NumGlobalDiagonals();
00314   
00315   TEST_EQUALITY(hdiag1, Ediag2);
00316 
00317   delete Matrix;
00318   delete TestMat;
00319 }
00320 
00321 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyNonzeros ) {
00322   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00323   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00324   
00325   int nnz1 = Matrix->NumMyNonzeros();
00326   int nnz2 = TestMat->NumMyNonzeros();
00327   
00328   TEST_EQUALITY(nnz1, nnz2);
00329 
00330   delete Matrix;
00331   delete TestMat;
00332 }
00333 
00334 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyRows ) {
00335   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00336   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00337   
00338   int rows1 = Matrix->NumMyRows();
00339   int rows2 = TestMat->NumMyRows();
00340   
00341   TEST_EQUALITY(rows1, rows2);
00342 
00343   delete Matrix;
00344   delete TestMat;
00345 }
00346 
00347 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyCols ) {
00348   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00349   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00350   
00351   int cols1 = Matrix->NumMyCols();
00352   int cols2 = TestMat->NumMyCols();
00353   
00354   TEST_EQUALITY(cols1, cols2);
00355 
00356   delete Matrix;
00357   delete TestMat;
00358 }
00359 
00360 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyDiagonals ) {
00361   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00362   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00363   
00364   int diag1 = Matrix->NumMyDiagonals();
00365   int diag2 = TestMat->NumMyDiagonals();
00366   
00367   TEST_EQUALITY(diag1, diag2);
00368 
00369   delete Matrix;
00370   delete TestMat;
00371 }
00372 
00373 TEUCHOS_UNIT_TEST( EpetraExt_hypre, MaxNumEntries ) {
00374   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00375   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00376   
00377   int ent1 = Matrix->MaxNumEntries();
00378   int ent2 = TestMat->MaxNumEntries();
00379   
00380   TEST_EQUALITY(ent1, ent2);
00381 
00382   delete Matrix;
00383   delete TestMat;
00384 }
00385 
00386 TEUCHOS_UNIT_TEST( EpetraExt_hypre, ApplyInverse ) {
00387   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 1);
00388   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00389 
00390   int num_vectors = 1;
00391   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00392   X.Random();
00393   Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
00394   Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
00395   
00396   TestMat->ApplyInverse(X,Y2);
00397   Matrix->ApplyInverse(X,Y1);
00398   TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol), true);
00399 
00400   delete Matrix;
00401   delete TestMat;
00402 }
00403 
00404 TEUCHOS_UNIT_TEST( EpetraExt_hypre, SameMatVec ) {
00405   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00406   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00407   
00408   int num_vectors = 3;
00409   Epetra_MultiVector X1(Matrix->RowMatrixRowMap(), num_vectors, true);
00410   X1.Random();
00411   Epetra_MultiVector X2(X1);
00412   
00413   Matrix->Multiply(false,X1,X1);
00414   TestMat->Multiply(false,X2,X2);
00415 
00416   TEST_EQUALITY(EquivalentVectors(X1,X2,tol), true); 
00417 
00418   delete Matrix;
00419   delete TestMat;
00420 }
00421 
00422 TEUCHOS_UNIT_TEST( EpetraExt_hypre, Solve ) {
00423   int num_vectors = 2;
00424   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 1);
00425   Epetra_MultiVector TrueX(Matrix->RowMatrixRowMap(), num_vectors, true);
00426   TrueX.Random();
00427   Epetra_MultiVector RHS(Matrix->RowMatrixRowMap(), num_vectors, true);
00428   Matrix->Multiply(false,TrueX, RHS);
00429   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00430   
00431   Matrix->SetParameter(Solver, PCG);
00432   Matrix->SetParameter(Preconditioner, BoomerAMG);
00433 
00434   /* Set some parameters (See Reference Manual for more parameters) */
00435   Matrix->SetParameter(Solver, &HYPRE_PCGSetMaxIter, 1000); /* max iterations */
00436   Matrix->SetParameter(Solver, &HYPRE_PCGSetTol, 1e-7); /* conv. tolerance */
00437   Matrix->SetParameter(Solver, &HYPRE_PCGSetTwoNorm, 1); /* use the two norm as the stopping criteria */
00438   Matrix->SetParameter(Solver, &HYPRE_PCGSetPrintLevel, 0); /* print solve info */
00439   Matrix->SetParameter(Solver, &HYPRE_PCGSetLogging, 0); /* needed to get run info later */
00440 
00441   /* Now set up the AMG preconditioner and specify any parameters */
00442   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetPrintLevel, 0); /* print amg solution info */
00443   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetCoarsenType, 6);
00444   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetRelaxType, 6); /* Sym G.S./Jacobi hybrid */ 
00445   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetNumSweeps, 1);
00446   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetTol, 0.0); /* conv. tolerance zero */
00447   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetMaxIter, 10); /* do only one iteration! */
00448 
00449   Matrix->SetParameter(true);
00450   /* Now setup and solve! */
00451   Matrix->Solve(false, false, false, RHS, X);
00452   TEST_EQUALITY(EquivalentVectors(X,TrueX,tol), true);
00453   
00454   delete Matrix;
00455 }
00456 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines