[petsc-dev] problem when creating a new ksp that calls another one

Matthew Knepley knepley at gmail.com
Wed Jul 15 15:01:23 CDT 2015 

On Wed, Jul 15, 2015 at 2:50 PM, Raphaël Couturier <
raphael.couturier at univ-fcomte.fr> wrote:

>  Hello,
>
> Sorry for my bad example...
>
> So maybe it is simpler I give you my complete code called TSIRM for A
> Two-Stage Iteration with least-squares Residual Minimization algorithm to
> solve large sparse linear systems (
> https://www.researchgate.net/publication/277571149_TSIRM_A_Two-Stage_Iteration_with_least-squares_Residual_Minimization_algorithm_to_solve_large_sparse_linear_systems
> )
> It is still under in development as a solver for PETSc (with few comments
> in the code...)
>
> In this code, I use FGMRES as an inner solver only for 30 iterations (I
> want to put this as a parameter after)
> For the outer iteration I apply a minimization method on the residual each
> 12 outer iterations (this is the parameter colS in the code, this will also
> be a parameter after). The minimization is CGLS conjugate gradient with
> least square.
> There are useless messages when running the code (but they are interesting
> to understand the behavior of the code)
>
> Some examples:
>
> With SNES ex35
>
> With FGMRES
> time  $PETSC_DIR/arch-linux2-c-debug/bin/mpirun -np 3 ./ex35 -snes_rtol
> 1.e-12 -snes_monitor  -ksp_type fgmres  -da_grid_x 660 -da_grid_y 660
>   0 SNES Function norm 3.808595655784e+02
>   1 SNES Function norm 3.785299355288e-03
>   2 SNES Function norm 3.766792406275e-08
>   3 SNES Function norm 1.428490284661e-09
>
> real    2m35.576s
> user    7m44.008s
> sys    0m2.504s
>
>
> With TSIRM
> time  $PETSC_DIR/arch-linux2-c-debug/bin/mpirun -np 3 ./ex35 -snes_rtol
> 1.e-12 -snes_monitor  -ksp_type tsirm  -da_grid_x 660 -da_grid_y 660
>   0 SNES Function norm 3.808595655784e+02
> Enter tsirm
> Size of vector 435600
> converged 0
> Norm of error 270.484, outer iteration 0 270.484 reason 0
> Norm of error 236.805, outer iteration 0 236.805 reason 0
> Norm of error 213.261, outer iteration 0 213.261 reason 0
> ....
> Norm of error 5.84204e-13, outer iteration 3 5.84204e-13 reason 0
> Norm of error 4.61092e-13, outer iteration 3 4.61092e-13 reason 0
> Norm of error 3.74294e-13, outer iteration 3 3.74294e-13 reason 2
>              -- Execution time : 22.2333 (s)
>              -- Total number of iterations : 1199
>   3 SNES Function norm 1.432052860297e-09
>
> real    0m53.639s
> user    2m39.120s
> sys    0m1.652s
>
>
>
>
>
> With KSP 15
>
> with FGMRES
>
>  time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun  -np 3    ./ex15 -m 800 -n
> 800 -ksp_type fgmres  -ksp_rtol 1e-6 -pc_type mg
> Norm of error 1.79776 iterations 1102
>
> real    0m51.147s
> user    2m32.788s
> sys    0m0.508s
>
>
> With TSIRM
>
>  time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun  -np 3    ./ex15 -m 800 -n
> 800 -ksp_type tsirm  -ksp_rtol 1e-6 -pc_type mg
> Enter tsirm
> Size of vector 640000
> converged 0
> Norm of error 0.0970815, outer iteration 0 0.0970815 reason 0
> Norm of error 0.0368809, outer iteration 0 0.0368809 reason 0
> Norm of error 0.0223725, outer iteration 0 0.0223725 reason 0
> Norm of error 0.0163448, outer iteration 0 0.0163448 reason 0
> Norm of error 0.013248, outer iteration 0 0.013248 reason 0
> Norm of error 0.0112125, outer iteration 0 0.0112125 reason 0
> Norm of error 0.00961675, outer iteration 0 0.00961675 reason 0
> Norm of error 0.00825048, outer iteration 0 0.00825048 reason 0
> Norm of error 0.00707865, outer iteration 0 0.00707865 reason 0
> Norm of error 0.00606273, outer iteration 0 0.00606273 reason 0
> Norm of error 0.00516968, outer iteration 0 0.00516968 reason 0
> Norm of error 0.00439978, outer iteration 0 0.00439978 reason 0
> Norm of error 5.22567e-05, outer iteration 1 5.22567e-05 reason 2
>              -- Execution time : 18.8489 (s)
>              -- Total number of iterations : 386
> Norm of error 0.901491 iterations 26
>
> real    0m19.263s
> user    0m57.248s
> sys    0m0.400s
>
>
>
>
> With KSP ex45 (with two KSPs not as in this code, we obtain faster results
> TSIRM with the ASM preconditioner than FGMRES with HYPRE or with ASM)
> With one KSP (this code)
> It crashes...
>

1) The traces for proc 0 and 1 are different, which is not good.

2) We lock vectors during a solve, so maybe you are reusing the rhs vector
in the inner solve? This is not
    a problem unless the example is using the ComputeRHS() interface, which
assembles the rhs again for
    each liner solve.

   Matt


>  time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun  -np 3    ./ex45
> -da_grid_x 160 -da_grid_y 160 -da_grid_z 160 -ksp_type tsirm  -ksp_rtol
> 1e-10 -pc_type asm -ksp_monitor
> Enter tsirm
> Size of vector 4096000
> [1]PETSC ERROR: --------------------- Error Message
> --------------------------------------------------------------
> [1]PETSC ERROR: Object is in wrong state
> [1]PETSC ERROR:  Vec is locked read only, argument # 1
> [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
> [1]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
> [1]PETSC ERROR: ./ex45 on a arch-linux2-c-debug named extinction by
> couturie Wed Jul 15 21:43:51 2015
> [1]PETSC ERROR: Configure options --download-openmpi --download-hypre
> --download-f2cblaslapack --with-fc=0
> [1]PETSC ERROR: #1 VecGetArray() line 1646 in
> /home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c
> [1]PETSC ERROR: #2 VecGetArray3d() line 2411 in
> /home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c
> [1]PETSC ERROR: #3 DMDAVecGetArray() line 77 in
> /home/couturie/work/petsc-3.6.0/src/dm/impls/da/dagetarray.c
> [1]PETSC ERROR: #4 ComputeRHS() line 84 in
> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c
> [1]PETSC ERROR: #5 KSPSetUp() line 277 in
> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c
> [1]PETSC ERROR: #6 KSPSolve_TSIRM() line 135 in
> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/impls/gmres/tsirm.c
>
>
>
> Thank you in advance,
>
> Raphaël
>
>
>   On Wed, Jul 15, 2015 at 12:04 PM, Raphaël Couturier <
> raphael.couturier at univ-fcomte.fr> wrote:
>
>> Hello,
>>
>> I want to make a new solver inside petsc.
>> This solver is based on a 2 stage method. So I have one outer iteration
>> and one inner iteration. For the inner iteration I can use gmres or one of
>> its variants.
>> I succeeded to make my algorithm by creating a new ksp (that copies the
>> linear system). So I have tried to do that without creating a new ksp. It
>> works well with ksp examples that don't use dmda and with snes examples.
>>
>
>  The code you sent does not make any sense to me. How can KSPSolve(ksp,
> NULL, NULL) be meaningful?
>
>  Also, can't you get exactly this effect using PCKSP?
>
>   Thanks,
>
>      Matt
>
>
>> So I have tried to make a really simple ksp (only to test) in which I
>> simply call another ksp.
>> Here is the part of the code (and attached is the complete code):
>>    PetscErrorCode ierr;
>>    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
>>    KSPSetType(ksp,KSPFGMRES);
>>    ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);
>>    PetscFunctionReturn(0);
>>
>>
>> For example with the ksp example 15:
>>  time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun  -np 3    ./ex15 -m 200
>> -n 200 -ksp_type test  -ksp_rtol 1e-10 -pc_type mg -ksp_monitor
>> there is no problem, I obtain the good result.
>>
>> But with the ksp example 45 (based on dmda):
>> time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun  -np 1    ./ex45 -ksp_type
>> test  -ksp_rtol 1e-10 -pc_type mg -ksp_monitor
>>
>> It crashes with that (I didn't change anything in ex45)
>> [0]PETSC ERROR: Object is in wrong state
>> [0]PETSC ERROR:  Vec is locked read only, argument # 1
>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
>> for trouble shooting.
>> [0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015
>> [0]PETSC ERROR: ./ex45 on a arch-linux2-c-debug named extinction by
>> couturie Wed Jul 15 18:45:39 2015
>> [0]PETSC ERROR: Configure options --download-openmpi --download-hypre
>> --download-f2cblaslapack --with-fc=0
>> [0]PETSC ERROR: #1 VecGetArray() line 1646 in
>> /home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c
>> [0]PETSC ERROR: #2 VecGetArray3d() line 2411 in
>> /home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c
>> [0]PETSC ERROR: #3 DMDAVecGetArray() line 77 in
>> /home/couturie/work/petsc-3.6.0/src/dm/impls/da/dagetarray.c
>> [0]PETSC ERROR: #4 ComputeRHS() line 84 in
>> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c
>> [0]PETSC ERROR: #5 KSPSetUp() line 277 in
>> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: #6 KSPSolve() line 546 in
>> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: #7 KSPSolve_TEST() line 43 in
>> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/impls/gmres/test.c
>> [0]PETSC ERROR: #8 KSPSolve() line 604 in
>> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: #9 main() line 51 in
>> /home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c
>>
>>
>> I don't understand since many examples in snes are good with my code and
>> they use dmda.
>> Do any of you have an idea of the problem?
>>
>> Of course, I have a bad solution that consists in creating a new ksp but
>> I would like to do without, if possible.
>>
>> Thank you in advance,
>> Regards,
>>
>> Raphaël Couturier
>>
>
>
>
>  --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20150715/a82bd303/attachment.html>

More information about the petsc-dev mailing list