Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Hypre/hypre_UnitTest.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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
43#include "Ifpack.h"
44#include "Ifpack_Hypre.h"
45#include "AztecOO.h"
46#include "Galeri_Maps.h"
47#include "Galeri_CrsMatrices.h"
48#include "Galeri_Utils.h"
49#include "Epetra_MultiVector.h"
50
51#include "Teuchos_UnitTestHarness.hpp"
52#include "Teuchos_RCP.hpp"
53#include "Epetra_ConfigDefs.h"
54#include "Epetra_Vector.h"
55#include "Epetra_RowMatrix.h"
56#include "Epetra_MultiVector.h"
57#include "Epetra_CrsMatrix.h"
58#include "Epetra_Map.h"
59#ifdef HAVE_MPI
60#include "mpi.h"
61#include "Epetra_MpiComm.h"
62#else
63#include "Epetra_SerialComm.h"
64#endif
65
66#include "hypre_Helpers.hpp"
67
68#include "Teuchos_Array.hpp"
69#include <string>
70#include <stdio.h>
71#include <map>
72
73using Teuchos::RCP;
74using Teuchos::rcp;
75
76
77// Tests hypre interface's ability to initialize correctly
78TEUCHOS_UNIT_TEST( Ifpack_Hypre, Construct ) {
79 const int N = 10;
80
81 Epetra_MpiComm comm(MPI_COMM_WORLD);
82 Epetra_Map RowMap(N, 0, comm);
83 Epetra_CrsMatrix Matrix(Copy, RowMap, 1);
84 for(int i = 0; i < N; i++){
85 int indices[1];
86 double values[1];
87 indices[0] = i;
88 values[0] = 1.0;
89 Matrix.InsertGlobalValues(i, 1, values, indices);
90 }
91 Matrix.FillComplete();
92 Ifpack_Hypre preconditioner(&Matrix);
93 int rv=preconditioner.Initialize();
94 TEST_EQUALITY(rv,0);
95}
96
97
98// Tests hypre's ability to work when A has a funky row map, but the vectors have
99// a contiguous row map
100TEUCHOS_UNIT_TEST( Ifpack_Hypre, ParameterList ){
101 const double tol = 1e-7;
102
103 Epetra_MpiComm Comm(MPI_COMM_WORLD);
104
105 Epetra_Map* Map;
106 // pointer to the matrix to be created
107 Epetra_CrsMatrix* Matrix;
108 // container for parameters
109 Teuchos::ParameterList GaleriList;
110 // here we specify the global dimension of the problem
111 int nx = 10 * Comm.NumProc();
112 int ny = 10 * Comm.NumProc();
113 GaleriList.set("nx", nx);
114 GaleriList.set("ny", ny);
115
116 try
117 {
118 // Create the matrix
119 Map = Galeri::CreateMap("Cartesian2D", Comm, GaleriList);
120 Matrix = Galeri::CreateCrsMatrix("Biharmonic2D", Map, GaleriList);
121
122 // Create the parameter list
123 Teuchos::ParameterList list("Preconditioner List");
124 Teuchos::ParameterList precList = list.sublist("hypre: Preconditioner functions");
125 precList.set("HYPRE_BoomerAMGSetPrintLevel", 1);// print amg solution info
126 precList.set("HYPRE_BoomerAMGSetCoarsenType", 6);
127 precList.set("HYPRE_BoomerAMGSetRelaxType", 6); //Sym G.S./Jacobi hybrid
128 precList.set("HYPRE_BoomerAMGSetNumSweeps", 1);
129 precList.set("HYPRE_BoomerAMGSetTol", 0.0); // conv. tolerance zero
130 precList.set("HYPRE_BoomerAMGSetMaxIter", 1); //do only one iteration!
131 list.set("hypre: Solver", "PCG");
132 list.set("hypre: Preconditioner", "BoomerAMG");
133 list.set("hypre: SolveOrPrecondition", "Solver");
134 list.set("hypre: SetPreconditioner", true);
135
136 // Create the preconditioner
137 Ifpack_Hypre preconditioner(Matrix);
138 TEST_EQUALITY(preconditioner.SetParameters(list),0);
139 TEST_EQUALITY(preconditioner.Compute(),0);
140
141 // Create the RHS and solution vector
142 int numVec = 2;
143 Epetra_MultiVector X(preconditioner.OperatorDomainMap(), numVec);
144 Epetra_MultiVector KnownX(preconditioner.OperatorDomainMap(), numVec);
145 TEST_EQUALITY(KnownX.Random(),0);
146 Epetra_MultiVector B(preconditioner.OperatorRangeMap(), numVec);
147 TEST_EQUALITY(preconditioner.Apply(KnownX, B),0);
148
149 TEST_EQUALITY(preconditioner.ApplyInverse(B,X),0);
150 TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,Comm.NumProc())), true);
151
152 //delete preconditioner;
153 delete Map;
154 delete Matrix;
155 }
156 catch (Galeri::Exception& rhs)
157 {
158 if (Comm.MyPID() == 0)
159 {
160 cerr << "Caught exception: ";
161 rhs.Print();
162 }
163 }
164}
165
166
167// Tests hypre's ability to work when A has a funky row map, but the vectors have
168// a contiguous row map
169TEUCHOS_UNIT_TEST( Ifpack_Hypre, ParameterList2 ){
170 const double tol = 1e-7;
171
172 Epetra_MpiComm Comm(MPI_COMM_WORLD);
173
174 Epetra_Map* Map;
175 // pointer to the matrix to be created
176 Epetra_CrsMatrix* Matrix;
177 // container for parameters
178 Teuchos::ParameterList GaleriList;
179 // here we specify the global dimension of the problem
180 int nx = 10 * Comm.NumProc();
181 int ny = 10 * Comm.NumProc();
182 GaleriList.set("nx", nx);
183 GaleriList.set("ny", ny);
184
185 try
186 {
187 // Create the matrix
188 Map = Galeri::CreateMap("Cartesian2D", Comm, GaleriList);
189 Matrix = Galeri::CreateCrsMatrix("Biharmonic2D", Map, GaleriList);
190
191 // Create the parameter list
192 Teuchos::ParameterList list("Preconditioner List");
193
194 list.set("hypre: Solver", "PCG");
195 Teuchos::ParameterList solverList = list.sublist("hypre: Solver functions");
196 solverList.set("HYPRE_PCGSetMaxIter", 1000); // max iterations
197 solverList.set("HYPRE_PCGSetTol", tol); // conv. tolerance
198 solverList.set("HYPRE_PCGSetTwoNorm", 1); // use the two norm as the stopping criteria
199 solverList.set("HYPRE_PCGSetPrintLevel", 0); // print solve info
200 solverList.set("HYPRE_PCGSetLogging", 1); // needed to get run info later
201 list.set("hypre: Solver functions",solverList);
202
203 list.set("hypre: Preconditioner", "BoomerAMG");
204 Teuchos::ParameterList precList = list.sublist("hypre: Preconditioner functions");
205 precList.set("HYPRE_BoomerAMGSetPrintLevel", 1);// print amg solution info
206 precList.set("HYPRE_BoomerAMGSetCoarsenType", 6);
207 precList.set("HYPRE_BoomerAMGSetRelaxType", 6); //Sym G.S./Jacobi hybrid
208 precList.set("HYPRE_BoomerAMGSetNumSweeps", 1);
209 precList.set("HYPRE_BoomerAMGSetTol", 0.0); // conv. tolerance zero
210 precList.set("HYPRE_BoomerAMGSetMaxIter", 1); //do only one iteration!
211 list.set("hypre: Preconditioner functions",precList);
212
213 list.set("hypre: SolveOrPrecondition", "Solver");
214 list.set("hypre: SetPreconditioner", true);
215
216 // Create the preconditioner
217 Ifpack_Hypre preconditioner(Matrix);
218 TEST_EQUALITY(preconditioner.SetParameters(list),0);
219 TEST_EQUALITY(preconditioner.Compute(),0);
220
221 // Create the RHS and solution vector
222 int numVec = 2;
223 Epetra_MultiVector X(preconditioner.OperatorDomainMap(), numVec);
224 Epetra_MultiVector KnownX(preconditioner.OperatorDomainMap(), numVec);
225 TEST_EQUALITY(KnownX.Random(),0);
226 Epetra_MultiVector B(preconditioner.OperatorRangeMap(), numVec);
227 TEST_EQUALITY(preconditioner.Apply(KnownX, B),0);
228
229 TEST_EQUALITY(preconditioner.ApplyInverse(B,X),0);
230 TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,Comm.NumProc())), true);
231
232 //delete preconditioner;
233 delete Map;
234 delete Matrix;
235 }
236 catch (Galeri::Exception& rhs)
237 {
238 if (Comm.MyPID() == 0)
239 {
240 cerr << "Caught exception: ";
241 rhs.Print();
242 }
243 }
244}
245
246
247// Tests the hypre interface's ability to work with both a preconditioner and linear
248// solver via ApplyInverse
249TEUCHOS_UNIT_TEST( Ifpack_Hypre, Ifpack ){
250 const double tol = 1E-7;
251
252 //
253 // Create Laplace2D with a contiguous row distribution
254 //
255 Epetra_CrsMatrix* Crs_Matrix = newCrsMatrix(4);
256
257 //
258 // Create the hypre preconditioner
259 //
260 Ifpack Factory;
261 RCP<Ifpack_Preconditioner> preconditioner= rcp(Factory.Create("Hypre", Crs_Matrix));
262 int rv = preconditioner->Initialize();
263 TEST_EQUALITY(rv,0);
264 int NumProc = Crs_Matrix->Comm().NumProc();
265 int MyPID = Crs_Matrix->Comm().MyPID();
266
267 //
268 // Create the solution vector and RHS
269 //
270 int numVec = 2;
271 Epetra_MultiVector X(preconditioner->OperatorRangeMap(), numVec);
272 Epetra_MultiVector KnownX(preconditioner->OperatorRangeMap(), numVec);
273 Epetra_MultiVector B(preconditioner->OperatorDomainMap(), numVec);
274 TEST_EQUALITY(KnownX.Random(), 0);
275 TEST_EQUALITY(preconditioner->Apply(KnownX, B), 0);
276
277 Teuchos::ParameterList list("New List");
278 Teuchos::ParameterList & solverList = list.sublist("hypre: Solver functions");
279 solverList.set("HYPRE_PCGSetMaxIter", 1000); // max iterations
280 solverList.set("HYPRE_PCGSetTol", 1e-9); // conv. tolerance
281 solverList.set("HYPRE_PCGSetPrintLevel", 0); // print solve info
282 solverList.set("HYPRE_PCGSetLogging", 1); // needed to get run info later
283 list.set("hypre: Solver functions",solverList);
284 list.set("hypre: SolveOrPrecondition", "Solver");
285 list.set("hypre: Solver", "PCG");
286 list.set("hypre: Preconditioner", "ParaSails");
287 list.set("hypre: SetPreconditioner", true);
288 TEST_EQUALITY(preconditioner->SetParameters(list), 0);
289 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, PCG); // Use a PCG Solver
290 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Preconditioner, ParaSails); // Use a ParaSails Preconditioner
291 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGSetMaxIter, 1000); // Set maximum iterations
292 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGSetTol, 1e-9); // Set a tolerance
293 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver); // Solve the problem
294 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(true); // Use the preconditioner
295 TEST_EQUALITY(preconditioner->Compute(),0);
296 preconditioner->Condest(Ifpack_Cheap, 1000, 1e-9, Crs_Matrix);
297 TEST_EQUALITY(preconditioner->ApplyInverse(B, X),0);
298 TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,NumProc)), true);
299 int numIters;
300 double residual;
301 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGGetNumIterations, &numIters);
302 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->SetParameter(Solver, &HYPRE_ParCSRPCGGetFinalRelativeResidualNorm, &residual);
303 (dynamic_cast<Ifpack_Hypre*> (preconditioner.get()))->CallFunctions();
304 if(MyPID == 0) printf("It took %d iterations, and achieved %e residual.\n", numIters, residual);
305 delete Crs_Matrix;
306}
307
308
309// This example uses contiguous maps, so hypre should not have problems
310TEUCHOS_UNIT_TEST( Ifpack_Hypre, DiagonalMatrixInOrder ) {
311 Epetra_MpiComm comm(MPI_COMM_WORLD);
312 int numProcs = comm.NumProc();
313 int myRank = comm.MyPID();
314
315 //
316 // Construct the contiguous map
317 //
318 Epetra_Map accMap(2*numProcs,0,comm);
319
320 //
321 // Construct the diagonal matrix
322 //
323 Epetra_CrsMatrix accMat(Copy,accMap,2,true);
324 int col;
325 double val;
326 col = 2*myRank;
327 val = 2*myRank+1;
328 accMat.InsertGlobalValues(col,1,&val,&col);
329 col = 2*myRank+1;
330 val = 2*myRank+2;
331 accMat.InsertGlobalValues(col,1,&val,&col);
332 accMat.FillComplete();
333
334 //
335 // Create the initial guess and RHS
336 //
337 int numVec = 2;
338 Epetra_MultiVector X(accMap, numVec);
339 Epetra_MultiVector KnownX(accMap, numVec);
340 KnownX.Random();
341 Epetra_MultiVector B(accMap, numVec);
342 accMat.Apply(KnownX, B);
343
344 //
345 // Create the parameter list
346 //
347 const double tol = 1e-7;
348 Teuchos::ParameterList list("Preconditioner List");
349 Teuchos::ParameterList & solverList = list.sublist("hypre: Solver functions");
350 solverList.set("HYPRE_PCGSetMaxIter", 1000); // max iterations
351 solverList.set("HYPRE_PCGSetTol", tol); // conv. tolerance
352 solverList.set("HYPRE_PCGSetPrintLevel", 0); // print solve info
353 solverList.set("HYPRE_PCGSetLogging", 1); // needed to get run info later
354 list.set("hypre: Solver", "PCG");
355 list.set("hypre: SolveOrPrecondition", "Solver");
356 list.set("hypre: SetPreconditioner", false);
357
358 //
359 // Create the preconditioner (which is actually a PCG solver)
360 //
361 Ifpack_Hypre prec(&accMat);
362 TEST_EQUALITY(prec.SetParameters(list),0);
363 TEST_EQUALITY(prec.Compute(),0);
364
365 //
366 // Solve the linear system
367 //
368 TEST_EQUALITY(prec.ApplyInverse(B,X),0);
369
370 TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,numProcs)), true);
371}
372
373
374// hypre does not like the distribution of the vectors in this example.
375// Our interface should detect that and remap as necessary.
376TEUCHOS_UNIT_TEST( Ifpack_Hypre, DiagonalMatrixOutOfOrder ) {
377 Epetra_MpiComm comm(MPI_COMM_WORLD);
378 int numProcs = comm.NumProc();
379 int myRank = comm.MyPID();
380
381 //
382 // Construct the map [0 0 1 1]
383 // Note that the rows will be out of order
384 //
385 int elementList[2];
386 elementList[0] = 2*myRank+1;
387 elementList[1] = 2*myRank;
388 Epetra_Map ncMap(2*numProcs,2,elementList,0,comm);
389
390 //
391 // Construct the diagonal matrix
392 //
393 Epetra_CrsMatrix accMat(Copy,ncMap,2,true);
394 int col;
395 double val;
396 col = 2*myRank+1;
397 val = 2*myRank+2;
398 accMat.InsertGlobalValues(col,1,&val,&col);
399 col = 2*myRank;
400 val = 2*myRank+1;
401 accMat.InsertGlobalValues(col,1,&val,&col);
402 accMat.FillComplete();
403
404 //
405 // Create the parameter list
406 //
407 const double tol = 1e-7;
408 Teuchos::ParameterList list("Preconditioner List");
409 Teuchos::ParameterList & solverList = list.sublist("hypre: Solver functions");
410 solverList.set("HYPRE_PCGSetMaxIter", 1000); // max iterations
411 solverList.set("HYPRE_PCGSetTol", tol); // conv. tolerance
412 solverList.set("HYPRE_PCGSetTwoNorm", 1); // conv. tolerance
413 solverList.set("HYPRE_PCGSetPrintLevel", 2); // print solve info
414 solverList.set("HYPRE_PCGSetLogging", 1); // needed to get run info later
415 list.set("hypre: Solver", "PCG");
416 list.set("hypre: SolveOrPrecondition", "Solver");
417 list.set("hypre: SetPreconditioner", false);
418
419
420 //
421 // Create the preconditioner (which is actually a PCG solver)
422 //
423 Ifpack_Hypre prec(&accMat);
424 TEST_EQUALITY(prec.SetParameters(list),0);
425 TEST_EQUALITY(prec.Compute(),0);
426
427 //
428 // Create the initial guess and RHS
429 // Note that these use non-contiguous maps
430 //
431 int numVec = 2;
432 Epetra_MultiVector X(ncMap, numVec);
433 Epetra_MultiVector KnownX(ncMap, numVec);
434 KnownX.Random();
435 Epetra_MultiVector B(ncMap, numVec);
436 accMat.Apply(KnownX, B);
437
438 //
439 // Solve the linear system
440 //
441 TEST_EQUALITY(prec.ApplyInverse(B,X),0);
442
443 TEST_EQUALITY(EquivalentVectors(X, KnownX, tol*10*pow(10.0,numProcs)), true);
444}
445
446
447
448// Creates the identity matrix with a non-contiguous row map
449// Even though the Epetra identity matrix has a map that hypre should not be happy with,
450// hypre should be able to redistribute it. It should also be able to accept the
451// vectors we give it, since they're using the same distribution as the hypre matrix.
452// This tests hypre's ability to perform as a linear solver via ApplyInverse.
453TEUCHOS_UNIT_TEST( Ifpack_Hypre, NonContiguousRowMap ) {
454
455 Epetra_MpiComm comm(MPI_COMM_WORLD);
456 int numProcs = comm.NumProc();
457 int myRank = comm.MyPID();
458
459 //
460 // Construct the map [0 1 0 1]
461 //
462 int elementList[2];
463 elementList[0] = myRank;
464 elementList[1] = myRank+numProcs;
465 Epetra_Map badMap(2*numProcs,2,elementList,0,comm);
466 Epetra_Map goodMap(2*numProcs,0,comm);
467
468 //
469 // Construct the identity matrix
470 //
471 Epetra_CrsMatrix badMat(Copy,badMap,2,true);
472 int col;
473 double val = 1.0;
474 col = myRank;
475 badMat.InsertGlobalValues(col,1,&val,&col);
476 col = myRank+numProcs;
477 badMat.InsertGlobalValues(col,1,&val,&col);
478 badMat.FillComplete();
479
480 //
481 // Create the parameter list
482 //
483 const double tol = 1e-7;
484 Teuchos::ParameterList list("Preconditioner List");
485 Teuchos::ParameterList & solverList = list.sublist("hypre: Solver functions");
486 solverList.set("HYPRE_PCGSetMaxIter", 1000); // max iterations
487 solverList.set("HYPRE_PCGSetTol", tol); // conv. tolerance
488 solverList.set("HYPRE_PCGSetTwoNorm", 1); // conv. tolerance
489 solverList.set("HYPRE_PCGSetPrintLevel", 2); // print solve info
490 solverList.set("HYPRE_PCGSetLogging", 1); // needed to get run info later
491 list.set("hypre: Solver", "PCG");
492 list.set("hypre: SolveOrPrecondition", "Solver");
493 list.set("hypre: SetPreconditioner", false);
494
495 //
496 // Create the preconditioner (which is actually a PCG solver)
497 //
498 Ifpack_Hypre prec(&badMat);
499 TEST_EQUALITY(prec.SetParameters(list),0);
500 TEST_EQUALITY(prec.Compute(),0);
501
502 //
503 // Create the initial guess and RHS
504 //
505 int numVec = 2;
506 Epetra_MultiVector X(goodMap, numVec);
507 Epetra_MultiVector RHS(goodMap, numVec);
508 RHS.Random();
509
510 //
511 // Solve the linear system
512 //
513 TEST_EQUALITY(prec.ApplyInverse(RHS,X),0);
514
515 TEST_EQUALITY(EquivalentVectors(X, RHS, tol*10*pow(10.0,numProcs)), true);
516}
Preconditioner
Solver
PCG
ParaSails
Copy
TEUCHOS_UNIT_TEST(Ifpack_Hypre, Construct)
@ Ifpack_Cheap
cheap estimate
std::pair< int, int > E
#define RHS(a)
Definition: MatGenFD.c:60
virtual int NumProc() const=0
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
virtual void Print(std::ostream &os) const
const Epetra_Comm & Comm() const
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int NumProc() const
int MyPID() const
static Ifpack_Preconditioner * Create(EPrecType PrecType, Epetra_RowMatrix *Matrix, const int overlap=0, bool overrideSerialDefault=false)
Creates an instance of Ifpack_Preconditioner given the enum value of the preconditioner type (can not...
Definition: Ifpack.cpp:253
bool EquivalentVectors(Epetra_MultiVector &Y1, Epetra_MultiVector &Y2, const double tol)
Epetra_CrsMatrix * newCrsMatrix(EpetraExt_HypreIJMatrix &Matrix)
const double tol
const int N