[petsc-users] PCJacobi applied to incompatible vector

Michael Souza souza.michael at gmail.com
Sat Jan 17 15:42:28 CST 2015

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?


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);
    ierr = MatSetUp(A); CHKERRQ(ierr);
    for (i = 0; i < na; i++) {
        if(i > 0)
            ierr = MatSetValue(A, i, (i - 1), 1.0, INSERT_VALUES);
        if(i < (na-1))
            ierr = MatSetValue(A, i, i+1, 1.0, INSERT_VALUES);
        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);
        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/3783e72d/attachment-0001.html>

More information about the petsc-users mailing list