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 (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "Teuchos_UnitTestHarness.hpp"
00043 #include "HYPRE_IJ_mv.h"
00044 #include "EpetraExt_HypreIJMatrix.h"
00045 #include "EpetraExt_MatrixMatrix.h"
00046 #include "Epetra_ConfigDefs.h"
00047 #include "Epetra_Vector.h"
00048 #include "Epetra_RowMatrix.h"
00049 #include "Epetra_MultiVector.h"
00050 #include "Epetra_CrsMatrix.h"
00051 #include "Epetra_Map.h"
00052 #ifdef HAVE_MPI
00053 #include "mpi.h"
00054 #include "Epetra_MpiComm.h"
00055 #else
00056 #include "Epetra_SerialComm.h"
00057 #endif
00058 
00059 #include "hypre_Helpers.hpp"
00060 #include "Teuchos_ParameterList.hpp"
00061 #include "Teuchos_ParameterEntry.hpp"
00062 #include "Teuchos_ParameterListExceptions.hpp"
00063 #include "Teuchos_Array.hpp"
00064 #include <string>
00065 #include <stdio.h>
00066 #include <map>
00067 
00068 namespace EpetraExt {
00069 
00070 const int N = 100;
00071 const int MatType = 4; //0 -> Unit diagonal, 1 -> Random diagonal, 2 -> Dense, val=col, 3 -> Random Dense, 4 -> Random Sparse
00072 const double tol = 1E-6;
00073 
00074 TEUCHOS_UNIT_TEST( EpetraExt_hypre, Construct ) {
00075 
00076   int ierr = 0;
00077   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 0);
00078   
00079   TEST_EQUALITY(Matrix->Filled(), true);
00080   
00081   for(int i = 0; i < Matrix->NumMyRows(); i++){
00082     int entries;
00083     ierr += Matrix->NumMyRowEntries(i, entries);
00084     TEST_EQUALITY(entries, 1);
00085     int numentries;
00086     Teuchos::Array<double> Values; Values.resize(entries);
00087     Teuchos::Array<int> Indices; Indices.resize(entries);
00088     ierr += Matrix->ExtractMyRowCopy(i, entries, numentries, &Values[0], &Indices[0]);
00089     TEST_EQUALITY(ierr, 0);
00090     TEST_EQUALITY(numentries,1);
00091     for(int j = 0; j < numentries; j++){
00092       TEST_FLOATING_EQUALITY(Values[j],1.0,tol);
00093       TEST_EQUALITY(Indices[j],i);
00094     }
00095   }
00096   delete Matrix;
00097 }
00098 
00099 TEUCHOS_UNIT_TEST( EpetraExt_hypre, MatVec ) {
00100   int ierr = 0;
00101   int num_vectors = 5;
00102   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 0);
00103   
00104   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00105   ierr += X.Random();
00106   Epetra_MultiVector Y(Matrix->RowMatrixRowMap(), num_vectors, true);
00107   
00108   ierr += Matrix->Multiply(false, X, Y);
00109   
00110   TEST_EQUALITY(EquivalentVectors(X,Y,tol),true);
00111   TEST_EQUALITY(ierr, 0);
00112   delete Matrix;
00113 }
00114 
00115 TEUCHOS_UNIT_TEST( EpetraExt_hypre, BetterMatVec ) {
00116   int ierr = 0;
00117   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00118   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00119   //TestMat->Print(std::cout);
00120   int num_vectors = 5;
00121   
00122   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00123   ierr += X.Random();
00124   Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
00125   Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
00126   
00127   ierr += Matrix->Multiply(false, X, Y1);
00128   ierr += TestMat->Multiply(false, X, Y2);
00129   
00130   TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol),true);
00131 
00132   ierr += Matrix->Multiply(false, Y1, X);
00133   ierr += TestMat->Multiply(false, Y1, Y2);
00134 
00135   TEST_EQUALITY(EquivalentVectors(X,Y2,tol),true);
00136 
00137   ierr += Matrix->Multiply(false, Y2, X);
00138   ierr += TestMat->Multiply(false, Y2, Y1);
00139 
00140   TEST_EQUALITY(EquivalentVectors(X,Y1,tol),true);
00141   TEST_EQUALITY_CONST(ierr, 0);
00142   delete Matrix;
00143   delete TestMat;
00144 }
00145 
00146 TEUCHOS_UNIT_TEST( EpetraExt_hypre, TransposeMatVec ) {
00147   int ierr = 0;
00148   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00149   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00150   
00151   int num_vectors = 5;
00152   
00153   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00154   ierr += X.Random();
00155   Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
00156   Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
00157   
00158   ierr += Matrix->Multiply(true, X, Y1);
00159   ierr += TestMat->Multiply(true, X, Y2);
00160   
00161   TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol),true);
00162   TEST_EQUALITY(ierr, 0);
00163   
00164   delete Matrix;
00165   delete TestMat;
00166 }
00167 
00168 TEUCHOS_UNIT_TEST( EpetraExt_hypre, LeftScale ) {
00169   int ierr = 0;
00170   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00171   //EpetraExt_HypreIJMatrix* BackUp = EpetraExt_HypreIJMatrix(Matrix);
00172   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00173   
00174   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00175   ierr += X.Random();
00176   Matrix->NumMyNonzeros();
00177   ierr += Matrix->LeftScale(X);
00178   ierr += TestMat->LeftScale(X);
00179   TEST_EQUALITY(EquivalentMatrices(*Matrix, *TestMat,tol), true);
00180   TEST_EQUALITY(ierr, 0);
00181   delete Matrix;
00182   delete TestMat;
00183 }
00184 
00185 TEUCHOS_UNIT_TEST( EpetraExt_hypre, RightScale ) {
00186   int ierr = 0;
00187   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00188   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00189   
00190   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00191   ierr += X.Random();
00192   Matrix->NumMyNonzeros();
00193   ierr += Matrix->RightScale(X);
00194   ierr += TestMat->RightScale(X);
00195   
00196   TEST_EQUALITY(EquivalentMatrices(*Matrix,*TestMat,tol), true);
00197   TEST_EQUALITY(ierr, 0);
00198   
00199   delete Matrix;
00200   delete TestMat;
00201 }
00202 
00203 TEUCHOS_UNIT_TEST( EpetraExt_hypre, ExtractDiagonalCopy ) {
00204   int ierr = 0;
00205   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00206   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00207   
00208   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00209   Epetra_Vector Y(TestMat->RowMatrixRowMap(),true);
00210   
00211   ierr += Matrix->ExtractDiagonalCopy(X);
00212   ierr += TestMat->ExtractDiagonalCopy(Y);
00213   
00214   TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
00215   TEST_EQUALITY(ierr, 0);
00216   
00217   delete Matrix;
00218   delete TestMat;
00219 }
00220 
00221 TEUCHOS_UNIT_TEST( EpetraExt_hypre, InvRowSums ) {
00222   int ierr = 0;
00223   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00224   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00225   
00226   Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
00227   Epetra_Vector Y(TestMat->RowMatrixRowMap(),true);
00228   
00229   ierr += Matrix->InvRowSums(X);
00230   ierr += TestMat->InvRowSums(Y);
00231   
00232   TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
00233   TEST_EQUALITY(ierr, 0);
00234   
00235   delete Matrix;
00236   delete TestMat;
00237 }
00238 
00239 TEUCHOS_UNIT_TEST( EpetraExt_hypre, InvColSums ) {
00240   int ierr = 0;
00241   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00242   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00243   
00244   Epetra_Vector X(Matrix->RowMatrixColMap(), true);
00245   Epetra_Vector Y(TestMat->RowMatrixColMap(),true);
00246   
00247   ierr += Matrix->InvColSums(X);
00248   ierr += TestMat->InvColSums(Y);
00249   
00250   TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
00251   TEST_EQUALITY(ierr, 0);
00252 
00253   delete Matrix;
00254   delete TestMat;
00255 }
00256 
00257 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NormInf ) {
00258   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00259   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00260   
00261   double norm1 = Matrix->NormInf();
00262   double norm2 = TestMat->NormInf();
00263   
00264   TEST_FLOATING_EQUALITY(norm1, norm2, tol);
00265   delete Matrix;
00266   delete TestMat;
00267 }
00268 
00269 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NormOne ) {
00270   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00271   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00272   
00273   double norm1 = Matrix->NormOne();
00274   double norm2 = TestMat->NormOne();
00275   
00276   TEST_FLOATING_EQUALITY(norm1, norm2, tol);
00277 
00278   delete Matrix;
00279   delete TestMat;
00280 }
00281 
00282 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalNonzeros ) {
00283   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00284   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00285   
00286   int nnz1 = Matrix->NumGlobalNonzeros();
00287   int nnz2 = TestMat->NumGlobalNonzeros();
00288   
00289   TEST_EQUALITY(nnz1, nnz2);
00290 
00291   delete Matrix;
00292   delete TestMat;
00293 }
00294 
00295 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalRows ) {
00296   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00297   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00298   
00299   int rows1 = Matrix->NumGlobalRows();
00300   int rows2 = TestMat->NumGlobalRows();
00301   
00302   TEST_EQUALITY(rows1, rows2);
00303 
00304   delete Matrix;
00305   delete TestMat;
00306 }
00307 
00308 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalCols ) {
00309   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00310   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00311   
00312   int cols1 = Matrix->NumGlobalCols();
00313   int cols2 = TestMat->NumGlobalCols();
00314   
00315   TEST_EQUALITY(cols1, cols2);
00316 
00317   delete Matrix;
00318   delete TestMat;
00319 }
00320 
00321 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalDiagonals ) {
00322   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00323   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00324   
00325   int hdiag1 = Matrix->NumGlobalDiagonals();
00326   int Ediag2 = TestMat->NumGlobalDiagonals();
00327   
00328   TEST_EQUALITY(hdiag1, Ediag2);
00329 
00330   delete Matrix;
00331   delete TestMat;
00332 }
00333 
00334 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyNonzeros ) {
00335   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00336   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00337   
00338   int nnz1 = Matrix->NumMyNonzeros();
00339   int nnz2 = TestMat->NumMyNonzeros();
00340   
00341   TEST_EQUALITY(nnz1, nnz2);
00342 
00343   delete Matrix;
00344   delete TestMat;
00345 }
00346 
00347 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyRows ) {
00348   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00349   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00350   
00351   int rows1 = Matrix->NumMyRows();
00352   int rows2 = TestMat->NumMyRows();
00353   
00354   TEST_EQUALITY(rows1, rows2);
00355 
00356   delete Matrix;
00357   delete TestMat;
00358 }
00359 
00360 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyCols ) {
00361   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00362   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00363   
00364   int cols1 = Matrix->NumMyCols();
00365   int cols2 = TestMat->NumMyCols();
00366   
00367   TEST_EQUALITY(cols1, cols2);
00368 
00369   delete Matrix;
00370   delete TestMat;
00371 }
00372 
00373 TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyDiagonals ) {
00374   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00375   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00376   
00377   int diag1 = Matrix->NumMyDiagonals();
00378   int diag2 = TestMat->NumMyDiagonals();
00379   
00380   TEST_EQUALITY(diag1, diag2);
00381 
00382   delete Matrix;
00383   delete TestMat;
00384 }
00385 
00386 TEUCHOS_UNIT_TEST( EpetraExt_hypre, MaxNumEntries ) {
00387   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00388   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00389   
00390   int ent1 = Matrix->MaxNumEntries();
00391   int ent2 = TestMat->MaxNumEntries();
00392   
00393   TEST_EQUALITY(ent1, ent2);
00394 
00395   delete Matrix;
00396   delete TestMat;
00397 }
00398 
00399 TEUCHOS_UNIT_TEST( EpetraExt_hypre, ApplyInverse ) {
00400   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 1);
00401   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00402 
00403   int num_vectors = 1;
00404   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00405   X.Random();
00406   Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
00407   Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
00408   
00409   TestMat->ApplyInverse(X,Y2);
00410   Matrix->ApplyInverse(X,Y1);
00411   TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol), true);
00412 
00413   delete Matrix;
00414   delete TestMat;
00415 }
00416 
00417 TEUCHOS_UNIT_TEST( EpetraExt_hypre, SameMatVec ) {
00418   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, MatType);
00419   Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
00420   
00421   int num_vectors = 3;
00422   Epetra_MultiVector X1(Matrix->RowMatrixRowMap(), num_vectors, true);
00423   X1.Random();
00424   Epetra_MultiVector X2(X1);
00425   
00426   Matrix->Multiply(false,X1,X1);
00427   TestMat->Multiply(false,X2,X2);
00428 
00429   TEST_EQUALITY(EquivalentVectors(X1,X2,tol), true); 
00430 
00431   delete Matrix;
00432   delete TestMat;
00433 }
00434 
00435 TEUCHOS_UNIT_TEST( EpetraExt_hypre, Solve ) {
00436   int num_vectors = 2;
00437   EpetraExt_HypreIJMatrix* Matrix = newHypreMatrix(N, 1);
00438   Epetra_MultiVector TrueX(Matrix->RowMatrixRowMap(), num_vectors, true);
00439   TrueX.Random();
00440   Epetra_MultiVector RHS(Matrix->RowMatrixRowMap(), num_vectors, true);
00441   Matrix->Multiply(false,TrueX, RHS);
00442   Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
00443   
00444   Matrix->SetParameter(Solver, PCG);
00445   Matrix->SetParameter(Preconditioner, BoomerAMG);
00446 
00447   /* Set some parameters (See Reference Manual for more parameters) */
00448   Matrix->SetParameter(Solver, &HYPRE_PCGSetMaxIter, 1000); /* max iterations */
00449   Matrix->SetParameter(Solver, &HYPRE_PCGSetTol, 1e-7); /* conv. tolerance */
00450   Matrix->SetParameter(Solver, &HYPRE_PCGSetTwoNorm, 1); /* use the two norm as the stopping criteria */
00451   Matrix->SetParameter(Solver, &HYPRE_PCGSetPrintLevel, 0); /* print solve info */
00452   Matrix->SetParameter(Solver, &HYPRE_PCGSetLogging, 0); /* needed to get run info later */
00453 
00454   /* Now set up the AMG preconditioner and specify any parameters */
00455   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetPrintLevel, 0); /* print amg solution info */
00456   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetCoarsenType, 6);
00457   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetRelaxType, 6); /* Sym G.S./Jacobi hybrid */ 
00458   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetNumSweeps, 1);
00459   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetTol, 0.0); /* conv. tolerance zero */
00460   Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetMaxIter, 10); /* do only one iteration! */
00461 
00462   Matrix->SetParameter(true);
00463   /* Now setup and solve! */
00464   Matrix->Solve(false, false, false, RHS, X);
00465   TEST_EQUALITY(EquivalentVectors(X,TrueX,tol), true);
00466   
00467   delete Matrix;
00468 }
00469 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines