Hello all,
First of all, thank you Matthew and Barry for your answers in July.
Your solution using PCKSP is fine.
Now it works perfectly with the examples I have tried for ksp and snes.
I made many tests with many processors.
I have implemented my own CGLS or LSQR method. In Petsc, I have seen
that CGNE and LSQR are already implemented. However, I cannot use them
because I have a problem of non square matrix.
In previous posts, I saw that the problem can come from the preconditioner.
In the attached code, I have my new solver (without preconditioner if I
made no mistake) which is based on
ierr = PCKSPGetKSP(pc,&sub_ksp);CHKERRQ(ierr);
and I need to use another ksp to minimize the residual (instead of my
own CGLS). The matrix AS is of size (n,10) with n very large.
I have tried with that (line 262 in the code)
{
KSP ksp3;
ierr = KSPCreate(PETSC_COMM_WORLD, &ksp3); CHKERRQ(ierr);
ierr = KSPSetType(ksp3, KSPLSQR); CHKERRQ(ierr);
ierr = KSPSetOperators(ksp3, AS, AS); CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp3, PETSC_DEFAULT, PETSC_DEFAULT,
PETSC_DEFAULT, iter_max_minimization);
ierr = PCSetType(ksp,PCNONE);
ierr = KSPSolve(ksp3, b, Alpha); CHKERRQ(ierr);
}
The problem on all the processors is:
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Only square matrices supported.
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015
[0]PETSC ERROR: ./ex15 on a arch-linux2-c-debug named extinction by
couturie Thu Sep 3 13:22:03 2015
[0]PETSC ERROR: Configure options --download-openmpi --download-hypre
--download-f2cblaslapack --with-fc=0
[0]PETSC ERROR: #2 MatGetDiagonalBlock_MPIDense() line 78 in
/home/couturie/work/petsc-3.6.1/src/mat/impls/dense/mpi/mpidense.c
[0]PETSC ERROR: #3 MatGetDiagonalBlock() line 159 in
/home/couturie/work/petsc-3.6.1/src/mat/interface/matrix.c
[0]PETSC ERROR: #4 PCSetUp_BJacobi() line 126 in
/home/couturie/work/petsc-3.6.1/src/ksp/pc/impls/bjacobi/bjacobi.c
[0]PETSC ERROR: #5 PCSetUp() line 982 in
/home/couturie/work/petsc-3.6.1/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #6 KSPSetUp() line 332 in
/home/couturie/work/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #7 KSPSolve() line 546 in
/home/couturie/work/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #8 KSPSolve_TSIRM() line 269 in
/home/couturie/work/petsc-3.6.1/src/ksp/ksp/impls/gmres/tsirm.c
[0]PETSC ERROR: #9 KSPSolve() line 604 in
/home/couturie/work/petsc-3.6.1/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #10 main() line 200 in
/home/couturie/work/petsc-3.6.1/src/ksp/ksp/examples/tutorials/ex15.c
Do you have an idea of the problem?
Is it a good idea to use the LSQR or CGLS method of Petsc?
Moreover I have another question. Do you think that the Petsc community
would be interested to use this solver. If yes, I would be glad to
improve it in order to integrate it in the source.
Thanks
Raphaël
#include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
#include <petscksp.h>
#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 ColS=12, col=0, iter_minimization=0, iter_max_minimization=15,outer_iter=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;
PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Tsirm","");
{
ierr = PetscOptionsInt("-ksp_tsirm_outer_iter","Number of outer iterations before minimization","",ColS,&ColS,NULL);CHKERRQ(ierr);
ierr = PetscOptionsInt("-ksp_tsirm_iter_min","Maximum number of iterations in the minimization process","",iter_max_minimization,&iter_max_minimization,NULL);CHKERRQ(ierr);
}
PetscOptionsEnd();
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);
KSP sub_ksp;
PC pc;
PCType type;
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCKSPGetKSP(pc,&sub_ksp);CHKERRQ(ierr);
ierr= KSPGetType(sub_ksp,&type);CHKERRQ(ierr);
ierr= PetscPrintf(PETSC_COMM_WORLD, "SUB KSP TYPE %s \n", type);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSetFromOptions(sub_ksp);CHKERRQ(ierr);
ierr = KSPSetTolerances(sub_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, restart); CHKERRQ(ierr);
//GMRES WITH MINIMIZATION
T1 = MPI_Wtime();
// ierr = KSPSetUp(sub_ksp); CHKERRQ(ierr);
//previously it seemed good but with SNES it seems not good....
ierr = KSP_MatMult(sub_ksp,A,x,r);CHKERRQ(ierr);
VecAXPY(r,-1,b);
VecNorm(r, NORM_2, &res);
ksp->its=0;
KSPConvergedDefault(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
//ierr = (*ksp->converged)(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "IIII converged %d residual norm %f its %d\n",ksp->reason,norm,ksp->its); CHKERRQ(ierr);
{
PetscReal rtol,abstol,dtol;
PetscInt maxits;
ierr = KSPGetTolerances(sub_ksp,&rtol,&abstol,&dtol,&maxits);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "sub ksp rtol %e abstol %e dtol %e maxits %d\n", rtol,abstol,dtol,maxits); CHKERRQ(ierr);
}
{
PetscReal rtol,abstol,dtol;
PetscInt maxits;
ierr = KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "ksp rtol %e abstol %e dtol %e maxits %d\n", rtol,abstol,dtol,maxits); CHKERRQ(ierr);
}
ierr = KSPSetInitialGuessNonzero(sub_ksp, PETSC_TRUE); CHKERRQ(ierr);
// PetscFinalize();
//exit(0);
do{
for(col=0; col<ColS && ksp->reason==0 ; col++){
//Solve
ierr = KSPSolve(sub_ksp, b, x); CHKERRQ(ierr);
ierr = KSPGetIterationNumber(sub_ksp, &its); CHKERRQ(ierr);
total += its;
//Build S'
ierr = VecGetArray(x, &array);
ierr = MatSetValues(S, Iend-Istart, ind_row, 1, &col, array, INSERT_VALUES);
VecRestoreArray(x, &array);
//ierr = KSPLogResidualHistory(sub_ksp,res);CHKERRQ(ierr);
//ierr = KSPMonitor(sub_ksp,sub_ksp->its,res);CHKERRQ(ierr);
KSPGetResidualNorm(sub_ksp,&norm);
ksp->rnorm = norm;
res=norm;
ksp->its++;
KSPConvergedDefault(ksp,ksp->its,res,&ksp->reason,ksp->cnvP);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, minimization iteration %D %g reason %D \n", res, outer_iter,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
iter_minimization = 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_minimization<iter_max_minimization){
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_minimization ++;
}
*/
{
KSP ksp3;
ierr = KSPCreate(PETSC_COMM_WORLD, &ksp3); CHKERRQ(ierr);
ierr = KSPSetType(ksp3, KSPLSQR); CHKERRQ(ierr);
ierr = KSPSetOperators(ksp3, AS, AS); CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp3, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, iter_max_minimization);
ierr = PCSetType(ksp,PCNONE);
ierr = KSPSolve(ksp3, b, Alpha); CHKERRQ(ierr);
}
outer_iter ++;
//Minimizer
MatMult(S, Alpha, x); //x = S*Alpha;
}
} while (iter_minimization<ksp->max_it && ksp->reason==0 );
ksp->its=total;
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);
}
