[petsc-users] PCJacobi applied to incompatible vector

[Barry Smith]

  As Dave points out the error is not detected because of a confluence of options you happen to be using; normally a MatMult() in the Krylov methods triggers an error but in your choices there is no MatMult(). The reason PCSetUp() and PCSetUp_BJacobi() cannot detect the problem is because they may not yet have available the vector that you pass into KSPSolve() hence cannot detect the size difference.

    I have added an error check in PCApply() to make sure the sizes of the objects passed in are correct; it is in the branch barry/fix-bjacobi-vector-lengths and next (I learned my lesson and did not stick it immediately into next and master :-(). I've also attached the patch file.

  Thanks

   Barry

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Preconditioner number of local rows 20 does not equal resulting vector number of rows 30
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.5.2, unknown 
[0]PETSC ERROR: ./ex1 on a arch-maint named Barrys-MacBook-Pro-3.local by barrysmith Sat Jan 17 17:15:12 2015
[0]PETSC ERROR: Configure options --download-mpich --with-fc=0 --with-cxx=0 --download-sowing=0
[0]PETSC ERROR: #1 PCApply() line 439 in /Users/barrysmith/Src/petsc/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #2 KSP_PCApply() line 230 in /Users/barrysmith/Src/petsc/include/petsc-private/kspimpl.h
[0]PETSC ERROR: #3 KSPSolve_PREONLY() line 26 in /Users/barrysmith/Src/petsc/src/ksp/ksp/impls/preonly/preonly.c
[0]PETSC ERROR: #4 KSPSolve() line 460 in /Users/barrysmith/Src/petsc/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #5 main() line 59 in /Users/barrysmith/Src/petsc/test-dir/ex1.c
[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint at mcs.anl.gov----------
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fix-bjacobi-vector-lengths.patch
Type: application/octet-stream
Size: 1269 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150117/541fe47a/attachment-0001.obj>
-------------- next part --------------

> On Jan 17, 2015, at 4:13 PM, Dave May <dave.mayhem23 at gmail.com> wrote:
> 
> 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;
> }
>