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

MLX mlx82 at fastwebnet.it
Mon Mar 3 11:45:26 MST 2008


Re: [Trilinos-Users] Problem with AztecOO (with AZ_reuse in nested loops)Thanks Mike. I have tried that way last week, but nothing change.

By the way, working on it another day, I've finally solved the problem. I think that this could be useful for someone in the future, so here is the correct order of the comands to use:

1) Before the first loop (computation of RHS for Schur complement) this parameters must be set for the inner loop:

List_Solver_int->set("pre_calc", "calc");
List_Solver_int->set("keep_info",1);

2) The first call to iterative method must be to the Aztec00 function member Iterate(maxIT, tol). So it compute the preconditioner and retain it.

3) Then the external loop begins. And here is the difference! In the inner loop is necessary to call the member function: recursiveIterate(). This will use the previous computed preconditioner, without giving any error.

4) At the end of the nested loop there is a last step to do: compute velocity. This step is similar to the first one:

- First is necessary to set: List_Solver_int->set("keep_info",0);
- Then is necessary to call: Iterate() and not recursiveIterate() (infact calling recursiveIterate() bring to the memory leak).

This way will avoid any error or memory leak.

I hope to be useful. Maybe this should be written somewere in the manual.

Thanks,

   Cristiano Malossi

  ----- Original Message ----- 
  From: Heroux, Michael A 
  To: mlx82 at fastwebnet.it ; trilinos-users at software.sandia.gov 
  Sent: Monday, March 03, 2008 5:56 PM
  Subject: Re: [Trilinos-Users] Problem with AztecOO (with AZ_reuse in nested loops)


  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-
    linux/vgpreload_memcheck.so)
    ==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
    error about AZ_reuse.

    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/a85c93f6/attachment-0001.html


More information about the Trilinos-Users mailing list