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
<[email protected]
<mailto:[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
#include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
#undef __FUNCT__
#define __FUNCT__ "KSPSetUp_TSIRM"
static PetscErrorCode KSPSetUp_TSIRM(KSP ksp)
{
PetscErrorCode ierr;
PetscBool diagonalscale;
PetscFunctionBegin;
//Maybe this need to be improved, I copied this from other solvers
ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
ierr = KSPSetWorkVecs(ksp,9);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "KSPSolve_TSIRM"
PetscErrorCode KSPSolve_TSIRM(KSP ksp) {
//Variables
PetscScalar gamma, alpha, oldgamma, beta;
PetscReal norm=20, cgprec=1e-40;
PetscInt giter=0, ColS=12, col=0, iter=0, iterations=15, Iiter=0;
PetscErrorCode ierr;
PetscScalar T1, T2;
PetscInt total=0;
PetscInt size;
PetscInt Istart,Iend;
PetscInt i,its;
Vec residu;
Mat S, AS;
PetscScalar *array;
PetscInt *ind_row;
Vec Alpha, p, ss, vect, r, q, Ax,x, b;
Mat A;
PetscInt first_iteration=1;
PetscReal res=99999999;
ierr = PetscPrintf(PETSC_COMM_WORLD, "Enter tsirm\n"); CHKERRQ(ierr);
KSPGetRhs(ksp,&b);
KSPGetOperators(ksp,&A,NULL);
KSPGetSolution(ksp,&x);
VecGetSize(b,&size);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Size of vector %D\n", size); CHKERRQ(ierr);
PetscInt aa,bb;
MatGetOwnershipRange(A,&aa,&bb);
ierr = MatCreate(PETSC_COMM_WORLD, &S);CHKERRQ(ierr);
ierr = MatSetSizes(S, bb-aa, PETSC_DECIDE, size, ColS); CHKERRQ(ierr);
ierr = MatSetType(S, MATMPIDENSE); CHKERRQ(ierr);
ierr = MatSetUp(S); CHKERRQ(ierr);
ierr = MatGetOwnershipRange(S, &Istart, &Iend); CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD, &Alpha); CHKERRQ(ierr);
ierr = VecSetSizes(Alpha, PETSC_DECIDE, ColS); CHKERRQ(ierr);
ierr = VecSetFromOptions(Alpha); CHKERRQ(ierr);
ierr = VecDuplicate(Alpha, &vect); CHKERRQ(ierr);
ierr = VecDuplicate(Alpha, &p); CHKERRQ(ierr);
ierr = VecDuplicate(Alpha, &ss); CHKERRQ(ierr);
ierr = VecDuplicate(b, &r); CHKERRQ(ierr);
ierr = VecDuplicate(b, &q); CHKERRQ(ierr);
ierr = VecDuplicate(b, &Ax); CHKERRQ(ierr);
ierr = VecDuplicate(b,&residu);CHKERRQ(ierr);
//indexes of row (these indexes are global)
ind_row = (PetscInt*)malloc(sizeof(PetscInt)*(Iend-Istart));
for(i=0; i<Iend-Istart; i++) ind_row[i] = i + Istart;
PetscInt restart=30;
PetscReal aa1,aa2,aa3;
PetscInt old_maxits;
ierr = KSPGetTolerances(ksp, &aa1, &aa2, &aa3, &old_maxits); CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, restart); CHKERRQ(ierr);
KSPSetType(ksp,KSPFGMRES);
ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr);
//GMRES WITH MINIMIZATION
T1 = MPI_Wtime();
ierr = KSPSetUp(ksp); CHKERRQ(ierr);
//previously it seemed good but with SNES it seems not good....
/* ierr = KSP_MatMult(ksp,A,x,r);CHKERRQ(ierr);
VecAXPY(r,-1,b);
VecNorm(r, NORM_2, &res);
ksp->reason=0;
ierr = PetscPrintf(PETSC_COMM_WORLD, "converged %d\n",ksp->reason); CHKERRQ(ierr);
*/
ierr = PetscPrintf(PETSC_COMM_WORLD, "converged %d\n",ksp->reason); CHKERRQ(ierr);
do{
for(col=0; col<ColS && ksp->reason==0 ; col++){
//Solve
ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
total += its; Iiter ++;
//Build S'
ierr = VecGetArray(x, &array);
ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
VecRestoreArray(x, &array);
KSPGetResidualNorm(ksp,&norm);
res=norm;
KSPConvergedDefault(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, outer iteration %D %g reason %D \n", norm, giter,ksp->rnorm, ksp->reason); CHKERRQ(ierr);
}
//minimization step
if(ksp->reason==0)
{
MatAssemblyBegin(S, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(S, MAT_FINAL_ASSEMBLY);
//Build the AS matrix
if(first_iteration) {
MatMatMult(A,S, MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AS);
first_iteration=0;
}
else
MatMatMult(A,S, MAT_REUSE_MATRIX,PETSC_DEFAULT,&AS);
//Minimization with the CGLS conjugate gradient least square method
MatMult(AS, Alpha, r); VecAYPX(r, -1, b); //r_0 = b-AS*x_0
MatMultTranspose(AS, r, p); //p_0 = AS'*r_0
ierr = VecCopy(p, ss); CHKERRQ(ierr); //p_0 = s_0
VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
while(gamma>cgprec && iter<iterations){
MatMult(AS, p, q); //q = AS*p
VecNorm(q, NORM_2, &alpha); alpha = alpha * alpha; //alpha = norm2(q)^2
alpha = gamma / alpha;
VecAXPY(Alpha, alpha, p); //Alpha += alpha*p
VecAXPY(r, -alpha, q); //r -= alpha*q
MatMultTranspose(AS, r, ss); // ss = AS'*r
oldgamma = gamma;
VecNorm(ss, NORM_2, &gamma); gamma = gamma * gamma; //gamma = norm2(s)^2
beta = gamma / oldgamma;
VecAYPX(p, beta, ss); //p = s+beta*p;
iter ++;
}
iter = 0; giter ++;
//Minimizer
MatMult(S, Alpha, x); //x = S*Alpha;
}
} while (giter<ksp->max_it && ksp->reason==0 );
//we need to put TSIRM otherwise for the next step, ksp will use FGMRES !!
KSPSetType(ksp,KSPTSIRM);
//we put the old max number of iterations (if an application need to change the ksp)
ierr = KSPSetTolerances(ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, old_maxits); CHKERRQ(ierr);
T2 = MPI_Wtime();
ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Execution time : %g (s)\n", T2-T1); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "\t\t\t -- Total number of iterations : %D\n", total); CHKERRQ(ierr);
//free allocated variables
ierr = VecDestroy(&Alpha);CHKERRQ(ierr);
ierr = MatDestroy(&S);CHKERRQ(ierr);
ierr = VecDestroy(&vect);CHKERRQ(ierr);
ierr = VecDestroy(&p);CHKERRQ(ierr);
ierr = VecDestroy(&ss);CHKERRQ(ierr);
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = VecDestroy(&q);CHKERRQ(ierr);
ierr = VecDestroy(&Ax);CHKERRQ(ierr);
ierr = VecDestroy(&residu);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "KSPCreate_TSIRM"
PETSC_EXTERN PetscErrorCode KSPCreate_TSIRM(KSP ksp)
{
PetscErrorCode ierr;
PetscFunctionBegin;
ksp->data = (void*)0;
ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);CHKERRQ(ierr);
ksp->ops->setup = KSPSetUp_TSIRM;
ksp->ops->solve = KSPSolve_TSIRM;
ksp->ops->destroy = KSPDestroyDefault;
ksp->ops->buildsolution = KSPBuildSolutionDefault;
ksp->ops->buildresidual = KSPBuildResidualDefault;
ksp->ops->setfromoptions = 0;
ksp->ops->view = 0;
#if defined(PETSC_USE_COMPLEX)
SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"This is not supported for complex numbers");
#endif
PetscFunctionReturn(0);
}
