<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, Jul 15, 2015 at 2:50 PM, Raphaël Couturier <span dir="ltr"><<a href="mailto:raphael.couturier@univ-fcomte.fr" target="_blank">raphael.couturier@univ-fcomte.fr</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div text="#000000" bgcolor="#FFFFFF">
<div>Hello,<br>
<br>
Sorry for my bad example...<br>
<br>
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
(<a href="https://www.researchgate.net/publication/277571149_TSIRM_A_Two-Stage_Iteration_with_least-squares_Residual_Minimization_algorithm_to_solve_large_sparse_linear_systems" target="_blank">https://www.researchgate.net/publication/277571149_TSIRM_A_Two-Stage_Iteration_with_least-squares_Residual_Minimization_algorithm_to_solve_large_sparse_linear_systems</a>)<br>
It is still under in development as a solver for PETSc (with few
comments in the code...)<br>
<br>
In this code, I use FGMRES as an inner solver only for 30
iterations (I want to put this as a parameter after)<br>
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.<br>
There are useless messages when running the code (but they are
interesting to understand the behavior of the code)<br>
<br>
Some examples:<br>
<br>
With SNES ex35<br>
<br>
With FGMRES<br>
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 <br>
0 SNES Function norm 3.808595655784e+02 <br>
1 SNES Function norm 3.785299355288e-03 <br>
2 SNES Function norm 3.766792406275e-08 <br>
3 SNES Function norm 1.428490284661e-09 <br>
<br>
real 2m35.576s<br>
user 7m44.008s<br>
sys 0m2.504s<br>
<br>
<br>
With TSIRM<br>
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 <br>
0 SNES Function norm 3.808595655784e+02 <br>
Enter tsirm<br>
Size of vector 435600<br>
converged 0<br>
Norm of error 270.484, outer iteration 0 270.484 reason 0 <br>
Norm of error 236.805, outer iteration 0 236.805 reason 0 <br>
Norm of error 213.261, outer iteration 0 213.261 reason 0 <br>
....<br>
Norm of error 5.84204e-13, outer iteration 3 5.84204e-13 reason 0
<br>
Norm of error 4.61092e-13, outer iteration 3 4.61092e-13 reason 0
<br>
Norm of error 3.74294e-13, outer iteration 3 3.74294e-13 reason 2
<br>
-- Execution time : 22.2333 (s)<br>
-- Total number of iterations : 1199<br>
3 SNES Function norm 1.432052860297e-09 <br>
<br>
real 0m53.639s<br>
user 2m39.120s<br>
sys 0m1.652s<br>
<br>
<br>
<br>
<br>
<br>
With KSP 15<br>
<br>
with FGMRES<br>
<br>
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<br>
Norm of error 1.79776 iterations 1102<br>
<br>
real 0m51.147s<br>
user 2m32.788s<br>
sys 0m0.508s<br>
<br>
<br>
With TSIRM<br>
<br>
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<br>
Enter tsirm<br>
Size of vector 640000<br>
converged 0<br>
Norm of error 0.0970815, outer iteration 0 0.0970815 reason 0 <br>
Norm of error 0.0368809, outer iteration 0 0.0368809 reason 0 <br>
Norm of error 0.0223725, outer iteration 0 0.0223725 reason 0 <br>
Norm of error 0.0163448, outer iteration 0 0.0163448 reason 0 <br>
Norm of error 0.013248, outer iteration 0 0.013248 reason 0 <br>
Norm of error 0.0112125, outer iteration 0 0.0112125 reason 0 <br>
Norm of error 0.00961675, outer iteration 0 0.00961675 reason 0 <br>
Norm of error 0.00825048, outer iteration 0 0.00825048 reason 0 <br>
Norm of error 0.00707865, outer iteration 0 0.00707865 reason 0 <br>
Norm of error 0.00606273, outer iteration 0 0.00606273 reason 0 <br>
Norm of error 0.00516968, outer iteration 0 0.00516968 reason 0 <br>
Norm of error 0.00439978, outer iteration 0 0.00439978 reason 0 <br>
Norm of error 5.22567e-05, outer iteration 1 5.22567e-05 reason 2
<br>
-- Execution time : 18.8489 (s)<br>
-- Total number of iterations : 386<br>
Norm of error 0.901491 iterations 26<br>
<br>
real 0m19.263s<br>
user 0m57.248s<br>
sys 0m0.400s<br>
<br>
<br>
<br>
<br>
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)<br>
With one KSP (this code)<br>
It crashes...<br></div></div></blockquote><div><br></div><div>1) The traces for proc 0 and 1 are different, which is not good.</div><div><br></div><div>2) We lock vectors during a solve, so maybe you are reusing the rhs vector in the inner solve? This is not</div><div> a problem unless the example is using the ComputeRHS() interface, which assembles the rhs again for</div><div> each liner solve.</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div text="#000000" bgcolor="#FFFFFF"><div>
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<br>
Enter tsirm<br>
Size of vector 4096000<br>
[1]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------<br>
[1]PETSC ERROR: Object is in wrong state<br>
[1]PETSC ERROR: Vec is locked read only, argument # 1<br>
[1]PETSC ERROR: See
<a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble
shooting.<br>
[1]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015 <br>
[1]PETSC ERROR: ./ex45 on a arch-linux2-c-debug named extinction
by couturie Wed Jul 15 21:43:51 2015<br>
[1]PETSC ERROR: Configure options --download-openmpi
--download-hypre --download-f2cblaslapack --with-fc=0<br>
[1]PETSC ERROR: #1 VecGetArray() line 1646 in
/home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c<br>
[1]PETSC ERROR: #2 VecGetArray3d() line 2411 in
/home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c<br>
[1]PETSC ERROR: #3 DMDAVecGetArray() line 77 in
/home/couturie/work/petsc-3.6.0/src/dm/impls/da/dagetarray.c<br>
[1]PETSC ERROR: #4 ComputeRHS() line 84 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c<br>
[1]PETSC ERROR: #5 KSPSetUp() line 277 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c<br>
[1]PETSC ERROR: #6 KSPSolve_TSIRM() line 135 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/impls/gmres/tsirm.c<br>
<br>
<br>
<br>
Thank you in advance,<br>
<br>
Raphaël<br>
<br>
<br>
</div>
<blockquote type="cite">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote">On Wed, Jul 15, 2015 at 12:04 PM,
Raphaël Couturier <span dir="ltr"><<a href="mailto:raphael.couturier@univ-fcomte.fr" target="_blank">raphael.couturier@univ-fcomte.fr</a>></span>
wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hello,<br>
<br>
I want to make a new solver inside petsc.<br>
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.<br>
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.<br>
</blockquote>
<div><br>
</div>
<div>The code you sent does not make any sense to me. How
can KSPSolve(ksp, NULL, NULL) be meaningful?</div>
<div><br>
</div>
<div>Also, can't you get exactly this effect using PCKSP?</div>
<div><br>
</div>
<div> Thanks,</div>
<div><br>
</div>
<div> Matt</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
So I have tried to make a really simple ksp (only to test)
in which I simply call another ksp.<br>
Here is the part of the code (and attached is the complete
code):<br>
PetscErrorCode ierr;<br>
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);<br>
KSPSetType(ksp,KSPFGMRES);<br>
ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);<br>
PetscFunctionReturn(0);<br>
<br>
<br>
For example with the ksp example 15:<br>
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<br>
there is no problem, I obtain the good result.<br>
<br>
But with the ksp example 45 (based on dmda):<br>
time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun -np 1
./ex45 -ksp_type test -ksp_rtol 1e-10 -pc_type mg
-ksp_monitor<br>
<br>
It crashes with that (I didn't change anything in ex45)<br>
[0]PETSC ERROR: Object is in wrong state<br>
[0]PETSC ERROR: Vec is locked read only, argument # 1<br>
[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a>
for trouble shooting.<br>
[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015<br>
[0]PETSC ERROR: ./ex45 on a arch-linux2-c-debug named
extinction by couturie Wed Jul 15 18:45:39 2015<br>
[0]PETSC ERROR: Configure options --download-openmpi
--download-hypre --download-f2cblaslapack --with-fc=0<br>
[0]PETSC ERROR: #1 VecGetArray() line 1646 in
/home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c<br>
[0]PETSC ERROR: #2 VecGetArray3d() line 2411 in
/home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c<br>
[0]PETSC ERROR: #3 DMDAVecGetArray() line 77 in
/home/couturie/work/petsc-3.6.0/src/dm/impls/da/dagetarray.c<br>
[0]PETSC ERROR: #4 ComputeRHS() line 84 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c<br>
[0]PETSC ERROR: #5 KSPSetUp() line 277 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: #6 KSPSolve() line 546 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: #7 KSPSolve_TEST() line 43 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/impls/gmres/test.c<br>
[0]PETSC ERROR: #8 KSPSolve() line 604 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: #9 main() line 51 in
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c<br>
<br>
<br>
I don't understand since many examples in snes are good
with my code and they use dmda.<br>
Do any of you have an idea of the problem?<br>
<br>
Of course, I have a bad solution that consists in creating
a new ksp but I would like to do without, if possible.<br>
<br>
Thank you in advance,<br>
Regards,<br>
<br>
Raphaël Couturier<br>
</blockquote>
</div>
<br>
<br clear="all"><span class="HOEnZb"><font color="#888888">
<div><br>
</div>
-- <br>
<div>What most experimenters take for
granted before they begin their experiments is infinitely
more interesting than any results to which their experiments
lead.<br>
-- Norbert Wiener</div>
</font></span></div>
</div>
</blockquote>
<br>
<div></div>
</div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>