# [Trilinos-Users] Problem with AztecOO (with AZ_reuse in nested loops)

Heroux, Michael A maherou at sandia.gov
Mon Mar 3 09:56:20 MST 2008

```Cristiano,

I cannot see anything obviously wrong with your approach, but you might consider a different one.  Once you have registered your linear system with the AztecOO object, you can call ConstructPreconditioner() to explicitly build the preconditioner.  The preconditioner should be retained until you call the DestroyPreconditioner() method.

Mike

On 2/29/08 10:03 AM, "mlx82 at fastwebnet.it" <mlx82 at fastwebnet.it> wrote:

Good evening,

I'm using Trilinos 8.0.5 (without MPI) to solve a Stokes problem. I'm
solving the system with the Schur complement and two nested loops.
As it is a "difficult" problem I need also some preconditioners.
In particoular I'm using two different ILUT preconditioners one for
the
inner loop and the other for the external one.

Everything is already implemented with AztecOO package (and it is
really
farster than IFPACK). I want to compute the preconditioner
for the inner loop only once, and reuse it for every iteration of
external loop.

To be as clear as possible I'm introducing here some notation; suppose
that the linear system is:

K * x = b

where:

K = [A B
C 0]

and:

x = [ V  P]^T;
b = [bV bP]^T;

The Schur problem is:
(-C * A^-1 * B) P = bP - C * A^-1 * bV

The resolution of the system is made by two steps:

1) Compute the RHS for the external loop: this is done solving a first
time a linear problem with A matrix and it's ILUT preconditioner (that
is computed now for the first time).

2) Compute P with the two nested loops: this time I want to reuse the
solver (and the preconditioner) used in step 1 for the inner loop. So
I
don't want to recompute the preconditioner.

Unfortunately after the first step I can't reuse the preconditioner of
matrix A in the inner loop. This is the error message that appears:

Error:  Did not find previous factorization (requested
by setting options[AZ_pre_calc] to AZ_reuse).
To find this factorization, the following
parameters must match the previous factorization:
1) Total number of unknowns.
2) options[AZ_overlap]
3) options[AZ_scaling]
4) options[AZ_precond]
5) options[AZ_reorder]
6) options[AZ_type_overlap]
7) options[AZ_subdomain_solve]
8) options[AZ_graph_fill]
9) params[AZ_ilut_fill]
10) params[AZ_drop]
11) data_org[AZ_name]
XXXP4302 '*! 1.0277XXX -7926 208
(   -7926,     208  32352016) ==> P4302 '*! 1.0277
( -914901,  269872  20123904) ==> vblock in gmres0
( -914901,    1632  19968448) ==> ptrs in gmres0
( -914901,   90760  15778480) ==> general in gmres0
(    2076,      80  32351840) ==> T332 '*! 1.0277
(    2076,     208  32351536) ==> P332 '*! 1.0277
(    2074,   34416  32097392) ==> x_reord P4302 '*! 1.0277
(    2074,  701688  14846224) ==> bindx P4302 '*! 1.0277
(    2074, 1403376  39390720) ==> val P4302 '*! 1.0277
(    2074,   17212  20028912) ==> iu P4302 '*! 1.0277
(    2074,   17212  20011600) ==> ordering P4302 '*! 1.0277
(    2074,   17212  19994288) ==> inv_ordering P4302 '*! 1.0277
(    2074,     224  32349920) ==> A_over P4302 '*! 1.0277
(    2074,      80  32349104) ==> T4302 '*! 1.0277
(    2074,     208  32348800) ==> P4302 '*! 1.0277

I have spent three days tryng to solve this problem without success.
The VERY strange thing is that if I try to solve again a problem with
matrix A before the beginning of step two (and external loop) the code
works, and reuse correctly the preconditioner. Also, if I re-compute
again
the preconditioner(after step one) during the first iteration of
external loop, and than reuse it for all the next iterations (inside
the external loop) the code works again. Unfortunately, with this
latter way I've not only to recompute the same preconditioner two
times, but also I've got a memory leakage! This is certanly introduced
by the
options keep_info = 1 in the internal loop. Valgrind gives me this
output:

==2038== searching for pointers to 9 not-freed blocks.
==2038== checked 3,851,536 bytes.
==2038==
==2038==
==2038== 2,192,376 bytes in 9 blocks are still reachable in loss
record 1 of 1
==2038==    at 0x4C21D06: malloc (in /usr/lib64/valgrind/amd64-
==2038==    by 0x503E9A: AZ_allocate (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x505C20: AZ_manage_memory (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x510626: AZ_mk_context (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x51092B: AZ_oldsolve_setup (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x5117F8: AZ_oldsolve (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x51257A: AZ_iterate (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x4FA178: AztecOO::Iterate(int, double) (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x454FF7: shurMatrix::Apply(Epetra_MultiVector const&,
Epetra_MultiVector&) const (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x4F5FF0: Epetra_Aztec_operatorvec (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x503360: AZ_compute_residual (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==    by 0x5083EF: AZ_scale_true_residual (in
/home/mlx/Steam3D/Steam3D_MLX/bin/Steam3D)
==2038==
==2038== LEAK SUMMARY:
==2038==    definitely lost: 0 bytes in 0 blocks.
==2038==      possibly lost: 0 bytes in 0 blocks.
==2038==    still reachable: 2,192,376 bytes in 9 blocks.
==2038==         suppressed: 0 bytes in 0 blocks.

To better understood this leakage I've tryed to solve the same problem
recomputing always the preconditioner and setting AZ_keep_info = 0,
AZ_pre_calc = "calc". And the leakage disapear.

Below I've reported a short description of the code that I'm using
now.

At the beginning I'm setting some parameters for the loops:

//Setting parameters for internal loop
List_Solver_int = new Teuchos::ParameterList("Solver: internal loop");
List_Solver_int->set("solver",   Sol_int_method);
List_Solver_int->set("conv",     "r0");
List_Solver_int->set("tol",      Sol_int_tol);
List_Solver_int->set("max_iter", Sol_int_ITmax);
List_Solver_int->set("kspace",   Sol_int_Kril);
List_Solver_int->set("output",   "last");
List_Solver_int->set("precond", "dom_decomp");
List_Solver_int->set("reorder", 1);
List_Solver_int->set("subdomain_solve", Prec_type);
List_Solver_int->set("overlap", 0);
List_Solver_int->set("type_overlap", "standard");
List_Solver_int->set("graph_fill", Prec_lof);
List_Solver_int->set("ilut_fill",  Prec_lof);
List_Solver_int->set("drop", Prec_dt);
List_Solver_int->set("pre_calc", "calc");
List_Solver_int->set("keep_info",1);

//Setting parameters for external loop
List_Solver_ext = new Teuchos::ParameterList("Solver: external loop");
List_Solver_ext->set("solver",   Sol_ext_method);
List_Solver_ext->set("conv",     "rhs");
List_Solver_ext->set("tol",      Sol_ext_tol);
List_Solver_ext->set("max_iter", Sol_ext_ITmax);
List_Solver_ext->set("kspace",   Sol_ext_Kril);
List_Solver_ext->set("output",   "none");
List_Solver_ext->set("precond", "dom_decomp");
List_Solver_ext->set("reorder", 1);
List_Solver_ext->set("subdomain_solve", Prec_type);
List_Solver_ext->set("overlap", 0);
List_Solver_ext->set("type_overlap", "standard");
List_Solver_ext->set("graph_fill", Prec_lof);
List_Solver_ext->set("ilut_fill",  Prec_lof);
List_Solver_ext->set("drop", Prec_dt);
List_Solver_ext->set("pre_calc", "calc");
List_Solver_ext->set("keep_info", 0);

Than I have a C++ class "invA" to solve the general problem: A * LHS =
RHS. The constructor of this class is:

solver.SetUserMatrix(A, flag);
solver.SetLHS(X);
solver.SetRHS(RHS);
solver.SetParameters(List_Solver_int);

where A, X, RHS are addresses, and solver is a private member (AztecOO
solver).
Therefore defining once the class "invA" I create only one solver for
this kind of problem. The main idea for this class is to change RHS
values outside the class without changing it's address.

To perform the first step I call a member function of "invA" that
does:

solver.Iterate(maxIT,tol);

this compute automatically the preconditioner.
Than I call a member function of "invA" that uses:

solver.SetAztecOption(AZ_pre_calc, AZ_reuse);

so, for the future call of solver.Iterate(maxIT,tol), it should reuse
the preconditioner!

To perform the second step and compute the pressure P, I define an
epetra operator "schurMatrix" in which is present the function
member "apply" that contains the operation that computes (-C * A^-1 *
B) * P,
and pointers to "invA" class and it's member.

I start the external loop (Aztec00 Sys_solver, with List_Solver_ext
options and appropiate LHS and RHS), after setting Sys_solver.
SetPrecMatrix(Mass) where Mass is a Mass matrix for compute external
ILUT preconditioner. Inside apply there is a call to solver.Iterate
(maxIT,tol)
that is inside "invA". And here the program ends with that

Thanks for any help.

Cristiano Malossi

_______________________________________________
Trilinos-Users mailing list
Trilinos-Users at software.sandia.gov
http://software.sandia.gov/mailman/listinfo/trilinos-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://software.sandia.gov/mailman/private/trilinos-users/attachments/20080303/41c4b7fe/attachment-0001.html
```