<div dir="ltr"><div><div>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.<br><br></div><div>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). <br></div><div><br></div>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().<br><br>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??<br><br><br></div><div>Cheers,<br></div><div> Dave<br></div><br><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On 17 January 2015 at 22:42, Michael Souza <span dir="ltr"><<a href="mailto:souza.michael@gmail.com" target="_blank">souza.michael@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_default" style="font-family:courier new,monospace">I'm not sure if this is a bug, but the PCJACOBI can be applied in a vector with size different from its operators.<br>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.<br><br>Is this just an unforeseen mistake (unexpected application)? If it isn't, what exactly does PCJACOBI do in this situation?<br><br>Cheers,<br>Michael<br><br>-------------------------------------------------------------------------------------<br>int main(int argc, char **args){<br> PetscErrorCode ierr;<br> PC pc;<br> KSP ksp;<br> MPI_Comm comm;<br> Mat A;<br> PetscInt i, na=20, nx=30;<br> Vec x, y;<br> <br> ierr = PetscInitialize(&argc, &args, (char *) 0, help); CHKERRQ(ierr);<br><br> // matrix creation and setup<br> ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRQ(ierr);<br> ierr = MatSetType(A, MATMPIAIJ); CHKERRQ(ierr);<br> ierr = MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, na, na); CHKERRQ(ierr);<br> ierr = MatSetUp(A); CHKERRQ(ierr);<br> for (i = 0; i < na; i++) {<br> if(i > 0)<br> ierr = MatSetValue(A, i, (i - 1), 1.0, INSERT_VALUES); CHKERRQ(ierr);<br> if(i < (na-1))<br> ierr = MatSetValue(A, i, i+1, 1.0, INSERT_VALUES); CHKERRQ(ierr);<br> ierr = MatSetValue(A, i, i, -2.0, INSERT_VALUES); CHKERRQ(ierr);<br> }<br> ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br> ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br> <br> // vector creation and setup<br> ierr = VecCreate(PETSC_COMM_WORLD, &x); CHKERRQ(ierr);<br> ierr = VecSetSizes(x,PETSC_DECIDE,nx);CHKERRQ(ierr);<br> ierr = VecSetFromOptions(x);CHKERRQ(ierr);<br> ierr = VecSet(x,1.0);CHKERRQ(ierr);<br> ierr = VecDuplicate(x,&y);CHKERRQ(ierr);<br> <br> // set KSP<br> ierr = PetscObjectGetComm((PetscObject) A, &comm); CHKERRQ(ierr);<br> ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr);<br> ierr = KSPSetType(ksp, KSPPREONLY); CHKERRQ(ierr);<br> ierr = KSPSetOperators(ksp, A, A); CHKERRQ(ierr);<br> ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);<br> ierr = PCSetType(pc, PCBJACOBI); CHKERRQ(ierr);<br> ierr = PCSetUp(pc); CHKERRQ(ierr);<br> { // setup the SubKSP<br> PetscInt nlocal, dummy;<br> KSP *subksp;<br> PC subpc;<br> ierr = PCBJacobiGetSubKSP(pc, &nlocal, &dummy, &subksp); CHKERRQ(ierr);<br> ierr = KSPSetType(*subksp, KSPPREONLY); CHKERRQ(ierr);<br> ierr = KSPGetPC(*subksp, &subpc); CHKERRQ(ierr);<br> ierr = PCSetType(subpc, PCILU); CHKERRQ(ierr);<br> ierr = PCFactorSetLevels(subpc, 1); CHKERRQ(ierr);<br> ierr = PCSetFromOptions(subpc); CHKERRQ(ierr);<br> ierr = PCSetUp(subpc); CHKERRQ(ierr);<br> }<br> <br> // solve<br> ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);<br> ierr = KSPSolve(ksp,x,y); CHKERRQ(ierr);<br> ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);<br> <br> // destroy objects<br> ierr = MatDestroy(&A); CHKERRQ(ierr);<br> ierr = VecDestroy(&x); CHKERRQ(ierr);<br> ierr = VecDestroy(&y); CHKERRQ(ierr);<br> ierr = KSPDestroy(&ksp); CHKERRQ(ierr);<br> <br> return EXIT_SUCCESS;<br>}<br></div></div>
</blockquote></div><br></div>