EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
hypre_UnitTest.cpp
Go to the documentation of this file.
1//@HEADER
2// ***********************************************************************
3//
4// EpetraExt: Epetra Extended - Linear Algebra Services Package
5// Copyright (2011) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41
42#include "Teuchos_UnitTestHarness.hpp"
43#include "HYPRE_IJ_mv.h"
46#include "Epetra_ConfigDefs.h"
47#include "Epetra_Vector.h"
48#include "Epetra_RowMatrix.h"
49#include "Epetra_MultiVector.h"
50#include "Epetra_CrsMatrix.h"
51#include "Epetra_Map.h"
52#ifdef HAVE_MPI
53#include "mpi.h"
54#include "Epetra_MpiComm.h"
55#else
56#include "Epetra_SerialComm.h"
57#endif
58
59#include "hypre_Helpers.hpp"
60#include "Teuchos_ParameterList.hpp"
61#include "Teuchos_ParameterEntry.hpp"
62#include "Teuchos_ParameterListExceptions.hpp"
63#include "Teuchos_Array.hpp"
64#include <string>
65#include <stdio.h>
66#include <map>
67
68namespace EpetraExt {
69
70const int N = 100;
71const int MatType = 4; //0 -> Unit diagonal, 1 -> Random diagonal, 2 -> Dense, val=col, 3 -> Random Dense, 4 -> Random Sparse
72const double tol = 1E-6;
73
74TEUCHOS_UNIT_TEST( EpetraExt_hypre, Construct ) {
75
76 int ierr = 0;
78
79 TEST_EQUALITY(Matrix->Filled(), true);
80
81 for(int i = 0; i < Matrix->NumMyRows(); i++){
82 int entries;
83 ierr += Matrix->NumMyRowEntries(i, entries);
84 TEST_EQUALITY(entries, 1);
85 int numentries;
86 Teuchos::Array<double> Values; Values.resize(entries);
87 Teuchos::Array<int> Indices; Indices.resize(entries);
88 ierr += Matrix->ExtractMyRowCopy(i, entries, numentries, &Values[0], &Indices[0]);
89 TEST_EQUALITY(ierr, 0);
90 TEST_EQUALITY(numentries,1);
91 for(int j = 0; j < numentries; j++){
92 TEST_FLOATING_EQUALITY(Values[j],1.0,tol);
93 TEST_EQUALITY(Indices[j],i);
94 }
95 }
96 delete Matrix;
97}
98
99TEUCHOS_UNIT_TEST( EpetraExt_hypre, MatVec ) {
100 int ierr = 0;
101 int num_vectors = 5;
103
104 Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
105 ierr += X.Random();
106 Epetra_MultiVector Y(Matrix->RowMatrixRowMap(), num_vectors, true);
107
108 ierr += Matrix->Multiply(false, X, Y);
109
110 TEST_EQUALITY(EquivalentVectors(X,Y,tol),true);
111 TEST_EQUALITY(ierr, 0);
112 delete Matrix;
113}
114
115TEUCHOS_UNIT_TEST( EpetraExt_hypre, BetterMatVec ) {
116 int ierr = 0;
118 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
119 //TestMat->Print(std::cout);
120 int num_vectors = 5;
121
122 Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
123 ierr += X.Random();
124 Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
125 Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
126
127 ierr += Matrix->Multiply(false, X, Y1);
128 ierr += TestMat->Multiply(false, X, Y2);
129
130 TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol),true);
131
132 ierr += Matrix->Multiply(false, Y1, X);
133 ierr += TestMat->Multiply(false, Y1, Y2);
134
135 TEST_EQUALITY(EquivalentVectors(X,Y2,tol),true);
136
137 ierr += Matrix->Multiply(false, Y2, X);
138 ierr += TestMat->Multiply(false, Y2, Y1);
139
140 TEST_EQUALITY(EquivalentVectors(X,Y1,tol),true);
141 TEST_EQUALITY_CONST(ierr, 0);
142 delete Matrix;
143 delete TestMat;
144}
145
146TEUCHOS_UNIT_TEST( EpetraExt_hypre, TransposeMatVec ) {
147 int ierr = 0;
149 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
150
151 int num_vectors = 5;
152
153 Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
154 ierr += X.Random();
155 Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
156 Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
157
158 ierr += Matrix->Multiply(true, X, Y1);
159 ierr += TestMat->Multiply(true, X, Y2);
160
161 TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol),true);
162 TEST_EQUALITY(ierr, 0);
163
164 delete Matrix;
165 delete TestMat;
166}
167
168TEUCHOS_UNIT_TEST( EpetraExt_hypre, LeftScale ) {
169 int ierr = 0;
171 //EpetraExt_HypreIJMatrix* BackUp = EpetraExt_HypreIJMatrix(Matrix);
172 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
173
174 Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
175 ierr += X.Random();
176 Matrix->NumMyNonzeros();
177 ierr += Matrix->LeftScale(X);
178 ierr += TestMat->LeftScale(X);
179 TEST_EQUALITY(EquivalentMatrices(*Matrix, *TestMat,tol), true);
180 TEST_EQUALITY(ierr, 0);
181 delete Matrix;
182 delete TestMat;
183}
184
185TEUCHOS_UNIT_TEST( EpetraExt_hypre, RightScale ) {
186 int ierr = 0;
188 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
189
190 Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
191 ierr += X.Random();
192 Matrix->NumMyNonzeros();
193 ierr += Matrix->RightScale(X);
194 ierr += TestMat->RightScale(X);
195
196 TEST_EQUALITY(EquivalentMatrices(*Matrix,*TestMat,tol), true);
197 TEST_EQUALITY(ierr, 0);
198
199 delete Matrix;
200 delete TestMat;
201}
202
203TEUCHOS_UNIT_TEST( EpetraExt_hypre, ExtractDiagonalCopy ) {
204 int ierr = 0;
206 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
207
208 Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
209 Epetra_Vector Y(TestMat->RowMatrixRowMap(),true);
210
211 ierr += Matrix->ExtractDiagonalCopy(X);
212 ierr += TestMat->ExtractDiagonalCopy(Y);
213
214 TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
215 TEST_EQUALITY(ierr, 0);
216
217 delete Matrix;
218 delete TestMat;
219}
220
221TEUCHOS_UNIT_TEST( EpetraExt_hypre, InvRowSums ) {
222 int ierr = 0;
224 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
225
226 Epetra_Vector X(Matrix->RowMatrixRowMap(), true);
227 Epetra_Vector Y(TestMat->RowMatrixRowMap(),true);
228
229 ierr += Matrix->InvRowSums(X);
230 ierr += TestMat->InvRowSums(Y);
231
232 TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
233 TEST_EQUALITY(ierr, 0);
234
235 delete Matrix;
236 delete TestMat;
237}
238
239TEUCHOS_UNIT_TEST( EpetraExt_hypre, InvColSums ) {
240 int ierr = 0;
242 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
243
244 Epetra_Vector X(Matrix->RowMatrixColMap(), true);
245 Epetra_Vector Y(TestMat->RowMatrixColMap(),true);
246
247 ierr += Matrix->InvColSums(X);
248 ierr += TestMat->InvColSums(Y);
249
250 TEST_EQUALITY(EquivalentVectors(X,Y,tol), true);
251 TEST_EQUALITY(ierr, 0);
252
253 delete Matrix;
254 delete TestMat;
255}
256
257TEUCHOS_UNIT_TEST( EpetraExt_hypre, NormInf ) {
259 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
260
261 double norm1 = Matrix->NormInf();
262 double norm2 = TestMat->NormInf();
263
264 TEST_FLOATING_EQUALITY(norm1, norm2, tol);
265 delete Matrix;
266 delete TestMat;
267}
268
269TEUCHOS_UNIT_TEST( EpetraExt_hypre, NormOne ) {
271 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
272
273 double norm1 = Matrix->NormOne();
274 double norm2 = TestMat->NormOne();
275
276 TEST_FLOATING_EQUALITY(norm1, norm2, tol);
277
278 delete Matrix;
279 delete TestMat;
280}
281
282TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalNonzeros ) {
284 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
285
286 int nnz1 = Matrix->NumGlobalNonzeros();
287 int nnz2 = TestMat->NumGlobalNonzeros();
288
289 TEST_EQUALITY(nnz1, nnz2);
290
291 delete Matrix;
292 delete TestMat;
293}
294
295TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalRows ) {
297 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
298
299 int rows1 = Matrix->NumGlobalRows();
300 int rows2 = TestMat->NumGlobalRows();
301
302 TEST_EQUALITY(rows1, rows2);
303
304 delete Matrix;
305 delete TestMat;
306}
307
308TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalCols ) {
310 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
311
312 int cols1 = Matrix->NumGlobalCols();
313 int cols2 = TestMat->NumGlobalCols();
314
315 TEST_EQUALITY(cols1, cols2);
316
317 delete Matrix;
318 delete TestMat;
319}
320
321TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumGlobalDiagonals ) {
323 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
324
325 int hdiag1 = Matrix->NumGlobalDiagonals();
326 int Ediag2 = TestMat->NumGlobalDiagonals();
327
328 TEST_EQUALITY(hdiag1, Ediag2);
329
330 delete Matrix;
331 delete TestMat;
332}
333
334TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyNonzeros ) {
336 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
337
338 int nnz1 = Matrix->NumMyNonzeros();
339 int nnz2 = TestMat->NumMyNonzeros();
340
341 TEST_EQUALITY(nnz1, nnz2);
342
343 delete Matrix;
344 delete TestMat;
345}
346
347TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyRows ) {
349 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
350
351 int rows1 = Matrix->NumMyRows();
352 int rows2 = TestMat->NumMyRows();
353
354 TEST_EQUALITY(rows1, rows2);
355
356 delete Matrix;
357 delete TestMat;
358}
359
360TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyCols ) {
362 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
363
364 int cols1 = Matrix->NumMyCols();
365 int cols2 = TestMat->NumMyCols();
366
367 TEST_EQUALITY(cols1, cols2);
368
369 delete Matrix;
370 delete TestMat;
371}
372
373TEUCHOS_UNIT_TEST( EpetraExt_hypre, NumMyDiagonals ) {
375 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
376
377 int diag1 = Matrix->NumMyDiagonals();
378 int diag2 = TestMat->NumMyDiagonals();
379
380 TEST_EQUALITY(diag1, diag2);
381
382 delete Matrix;
383 delete TestMat;
384}
385
386TEUCHOS_UNIT_TEST( EpetraExt_hypre, MaxNumEntries ) {
388 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
389
390 int ent1 = Matrix->MaxNumEntries();
391 int ent2 = TestMat->MaxNumEntries();
392
393 TEST_EQUALITY(ent1, ent2);
394
395 delete Matrix;
396 delete TestMat;
397}
398
399TEUCHOS_UNIT_TEST( EpetraExt_hypre, ApplyInverse ) {
401 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
402
403 int num_vectors = 1;
404 Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
405 X.Random();
406 Epetra_MultiVector Y1(Matrix->RowMatrixRowMap(), num_vectors, true);
407 Epetra_MultiVector Y2(Matrix->RowMatrixRowMap(), num_vectors, true);
408
409 TestMat->ApplyInverse(X,Y2);
410 Matrix->ApplyInverse(X,Y1);
411 TEST_EQUALITY(EquivalentVectors(Y1,Y2,tol), true);
412
413 delete Matrix;
414 delete TestMat;
415}
416
417TEUCHOS_UNIT_TEST( EpetraExt_hypre, SameMatVec ) {
419 Epetra_CrsMatrix* TestMat = newCrsMatrix(*Matrix);
420
421 int num_vectors = 3;
422 Epetra_MultiVector X1(Matrix->RowMatrixRowMap(), num_vectors, true);
423 X1.Random();
424 Epetra_MultiVector X2(X1);
425
426 Matrix->Multiply(false,X1,X1);
427 TestMat->Multiply(false,X2,X2);
428
429 TEST_EQUALITY(EquivalentVectors(X1,X2,tol), true);
430
431 delete Matrix;
432 delete TestMat;
433}
434
435TEUCHOS_UNIT_TEST( EpetraExt_hypre, Solve ) {
436 int num_vectors = 2;
438 Epetra_MultiVector TrueX(Matrix->RowMatrixRowMap(), num_vectors, true);
439 TrueX.Random();
440 Epetra_MultiVector RHS(Matrix->RowMatrixRowMap(), num_vectors, true);
441 Matrix->Multiply(false,TrueX, RHS);
442 Epetra_MultiVector X(Matrix->RowMatrixRowMap(), num_vectors, true);
443
444 Matrix->SetParameter(Solver, PCG);
446
447 /* Set some parameters (See Reference Manual for more parameters) */
448 Matrix->SetParameter(Solver, &HYPRE_PCGSetMaxIter, 1000); /* max iterations */
449 Matrix->SetParameter(Solver, &HYPRE_PCGSetTol, 1e-7); /* conv. tolerance */
450 Matrix->SetParameter(Solver, &HYPRE_PCGSetTwoNorm, 1); /* use the two norm as the stopping criteria */
451 Matrix->SetParameter(Solver, &HYPRE_PCGSetPrintLevel, 0); /* print solve info */
452 Matrix->SetParameter(Solver, &HYPRE_PCGSetLogging, 0); /* needed to get run info later */
453
454 /* Now set up the AMG preconditioner and specify any parameters */
455 Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetPrintLevel, 0); /* print amg solution info */
456 Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetCoarsenType, 6);
457 Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetRelaxType, 6); /* Sym G.S./Jacobi hybrid */
458 Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetNumSweeps, 1);
459 Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetTol, 0.0); /* conv. tolerance zero */
460 Matrix->SetParameter(Preconditioner, &HYPRE_BoomerAMGSetMaxIter, 10); /* do only one iteration! */
461
462 Matrix->SetParameter(true);
463 /* Now setup and solve! */
464 Matrix->Solve(false, false, false, RHS, X);
465 TEST_EQUALITY(EquivalentVectors(X,TrueX,tol), true);
466
467 delete Matrix;
468}
469} // namespace EpetraExt
int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a EpetraExt_HypreIJMatrix multiplied by a Epetra_MultiVector X in Y.
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
int LeftScale(const Epetra_Vector &X)
Scales the EpetraExt_HypreIJMatrix on the left with a Epetra_Vector x.
int SetParameter(Hypre_Chooser chooser, int(*pt2Func)(HYPRE_Solver, int), int parameter)
Set a parameter that takes a single int.
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a EpetraExt_HypreIJMatrix solving a Epetra_MultiVector X in Y.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
int RightScale(const Epetra_Vector &X)
Scales the EpetraExt_HypreIJMatrix on the right with a Epetra_Vector x.
virtual int InvColSums(Epetra_Vector &x) const
virtual int NumGlobalDiagonals() const
virtual int NumMyDiagonals() const
virtual int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
virtual int NumMyNonzeros() const
virtual const Epetra_Map & RowMatrixColMap() const
virtual int InvRowSums(Epetra_Vector &x) const
virtual int NumMyCols() const
virtual double NormInf() const
virtual int MaxNumEntries() const
virtual int NumGlobalCols() const
virtual int NumMyRows() const
virtual double NormOne() const
virtual int NumGlobalRows() const
virtual const Epetra_Map & RowMatrixRowMap() const
virtual bool Filled() const
virtual int NumGlobalNonzeros() const
bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol)
Epetra_CrsMatrix * newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
EpetraExt_HypreIJMatrix * newHypreMatrix(const int N, const int type)
bool EquivalentMatrices(Epetra_RowMatrix &HypreMatrix, Epetra_RowMatrix &CrsMatrix, const double tol)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
const int MatType
TEUCHOS_UNIT_TEST(EpetraExt_hypre, Construct)
const double tol
const int N