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

Barry Smith bsmith at mcs.anl.gov
Wed Jul 15 15:03:41 CDT 2015 

  You cannot have your KSPSolve_TSIRM() change KSP options on ITS KSP and then call KSPSetUp() on its KSP then call KSPSolve() on its KSP. If you want to have an "inner" solver you need to have a different KSP for that inner solver.

  Barry



> On 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...
> 
>  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
> 
> <tsirm.c>

More information about the petsc-dev mailing list