On Wed, Jul 15, 2015 at 2:50 PM, Raphaël Couturier < [email protected]> 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 < > [email protected]> 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
