[petsc-users] PCJacobi applied to incompatible vector

[Dave May]

Many of the Mat-Vec operations will check that the local and global sizes
of the matrix and vector are compatible, for example MatMult. If you
changed your outer KSP type to anything other than PREONLY, an error about
the mismatched sizes would get thrown.

It seems that none of the methods being called from within PCSetUp_BJACOBI
perform such a check (which is a bit of a surprise). I cannot tell you
exactly what the result of calling PCApply() with BJACOBI will be in the
case when the local/global sizes of x,y are larger than the those for the
operator. (I would guess that the "extra" entries might be ignored by
PCApply_BJACOBI, but without going through the function in detail this
really is a guess).

The call to PCApply() does not check that the local/global sizes of the
input and output vectors are compatible with the sizes associated with the
preconditioned operator. It seems neither does PCSetUp().

Shouldn't PCSetUp() perform such a check to detect this error early on and
to not have the in-built assumption that a method called within
PCSetUp_XXX() will in fact catch the error??


Cheers,
  Dave



On 17 January 2015 at 22:42, Michael Souza <souza.michael at gmail.com> wrote:

> I'm not sure if this is a bug, but the PCJACOBI can be applied in a vector
> with size different from its operators.
> In the code below, I define a PCJACOBI for a matrix A with size equal to
> 20 and I apply it in a 30 length vector x.
>
> Is this just an unforeseen mistake (unexpected application)? If it isn't,
> what exactly does PCJACOBI do in this situation?
>
> Cheers,
> Michael
>
>
> -------------------------------------------------------------------------------------
> int main(int argc, char **args){
>     PetscErrorCode ierr;
>     PC pc;
>     KSP ksp;
>     MPI_Comm comm;
>     Mat A;
>     PetscInt i, na=20, nx=30;
>     Vec x, y;
>
>     ierr = PetscInitialize(&argc, &args, (char *) 0, help); CHKERRQ(ierr);
>
>     // matrix creation and setup
>     ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRQ(ierr);
>     ierr = MatSetType(A, MATMPIAIJ); CHKERRQ(ierr);
>     ierr = MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, na, na);
> CHKERRQ(ierr);
>     ierr = MatSetUp(A); CHKERRQ(ierr);
>     for (i = 0; i < na; i++) {
>         if(i > 0)
>             ierr = MatSetValue(A, i, (i - 1), 1.0, INSERT_VALUES);
> CHKERRQ(ierr);
>         if(i < (na-1))
>             ierr = MatSetValue(A, i, i+1, 1.0, INSERT_VALUES);
> CHKERRQ(ierr);
>         ierr = MatSetValue(A, i, i, -2.0, INSERT_VALUES); CHKERRQ(ierr);
>     }
>     ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>     ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
>     // vector creation and setup
>     ierr = VecCreate(PETSC_COMM_WORLD, &x); CHKERRQ(ierr);
>     ierr = VecSetSizes(x,PETSC_DECIDE,nx);CHKERRQ(ierr);
>     ierr = VecSetFromOptions(x);CHKERRQ(ierr);
>     ierr = VecSet(x,1.0);CHKERRQ(ierr);
>     ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
>
>     // set KSP
>     ierr = PetscObjectGetComm((PetscObject) A, &comm); CHKERRQ(ierr);
>     ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr);
>     ierr = KSPSetType(ksp, KSPPREONLY); CHKERRQ(ierr);
>     ierr = KSPSetOperators(ksp, A, A); CHKERRQ(ierr);
>     ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
>     ierr = PCSetType(pc, PCBJACOBI); CHKERRQ(ierr);
>     ierr = PCSetUp(pc); CHKERRQ(ierr);
>     { // setup the SubKSP
>         PetscInt nlocal, dummy;
>         KSP *subksp;
>         PC subpc;
>         ierr = PCBJacobiGetSubKSP(pc, &nlocal, &dummy, &subksp);
> CHKERRQ(ierr);
>         ierr = KSPSetType(*subksp, KSPPREONLY); CHKERRQ(ierr);
>         ierr = KSPGetPC(*subksp, &subpc); CHKERRQ(ierr);
>         ierr = PCSetType(subpc, PCILU); CHKERRQ(ierr);
>         ierr = PCFactorSetLevels(subpc, 1); CHKERRQ(ierr);
>         ierr = PCSetFromOptions(subpc); CHKERRQ(ierr);
>         ierr = PCSetUp(subpc); CHKERRQ(ierr);
>     }
>
>     // solve
>     ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);
>     ierr = KSPSolve(ksp,x,y); CHKERRQ(ierr);
>     ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
>
>     // destroy objects
>     ierr = MatDestroy(&A); CHKERRQ(ierr);
>     ierr = VecDestroy(&x); CHKERRQ(ierr);
>     ierr = VecDestroy(&y); CHKERRQ(ierr);
>     ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
>
>     return EXIT_SUCCESS;
> }
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150117/45949c32/attachment.html>