[petsc-users] MatMult with one sequential and one parallel vector

[Claudio Kozický]

Hello,

I am using PETSc in a performance comparison that evaluates the
performance of parallel sparse matrix-vector multiplication (SpMV).  For
this purpose I have implemented a simple SpMV operation using PETSc,
which multiplies parallel matrix A (type MatAIJ) with parallel vector x1
and stores the result in parallel vector y1.  Thus I perform SpMV using
PETSc as MatMult(A, x1, y1).  This part works without any problems.

I would also like to implement SpMV operation y2 = A * x2, where x2 is a
sequential vector (i.e. created using VecCreateSeq) but where A and y2
are still parallel.  The resulting implementation would be something
like:

MatCreateAIJ(..., &A);     // a parallel matrix
VecCreateSeq(..., &x2);    // a per-process _sequential_ vector
VecCreateMPI(..., &y2);    // a parallel vector
MatMult(A, x2, y2);

The motivation of storing all of x2 in each process to is remove the
need of broadcasting any elements of x2 (this approach makes sense in
the context of what I am benchmarking).  However I cannot seem to get
this approach to work in PETSc.  For example when I try this approach
with a 4-by-4 matrix, with two order-4 vectors and using two MPI
processes, then PETSc prints:

Nonconforming object sizes
Mat mat,Vec x: local dim 2 4

I have attached a minimal working example that demonstrates what I am
attempting to perform.  Could it be that PETSc does not support
combining a parallel and sequential vector in a single MatMult call?

I have found functions for scattering and gathering vectors in the
documentation of PETSc, but these do not seem to be a good match for
what I am trying to benchmark.  My intention is for each process to keep
an identical copy of vector x2 and therefore the necessity to scatter or
gather values in x2 should never arise.

I would appreciate if somebody could help point me in the right
direction regarding my failing MatMult call.  Thanks!

-- 
Claudio Kozický
-------------- next part --------------
#include <petscksp.h>

int main(int argc, char **argv)
{
    PetscErrorCode ierr = PetscInitialize(&argc, &argv, NULL, NULL);
    if (ierr != 0) {
        return 1;
    }

    Mat A;
    ierr = MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 4, 4, 1, NULL, 0, NULL, &A);
    CHKERRQ(ierr);
    PetscInt rows[] = {0, 1, 2, 3};
    PetscInt cols[] = {0, 1, 2, 3};
    PetscScalar vals[] = {10, 20, 30, 40};
    for (int i = 0; i < 4; ++i) {
        ierr = MatSetValues(A, 1, &rows[i], 1, &cols[i], &vals[i], INSERT_VALUES);
        CHKERRQ(ierr);
    }
    ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    CHKERRQ(ierr);

    Vec x1;
    ierr = VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 4, &x1);
    CHKERRQ(ierr);
    ierr = VecSet(x1, 10);
    CHKERRQ(ierr);
    Vec y1;
    ierr = VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 4, &y1);
    CHKERRQ(ierr);

    ierr = MatMult(A, x1, y1);
    CHKERRQ(ierr);
    ierr = VecView(y1, PETSC_VIEWER_STDOUT_WORLD);
    CHKERRQ(ierr);

    Vec x2;
    ierr = VecCreateSeq(PETSC_COMM_SELF, 4, &x2);
    CHKERRQ(ierr);
    ierr = VecSet(x2, 10);
    CHKERRQ(ierr);
    Vec y2;
    ierr = VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 4, &y2);
    CHKERRQ(ierr);

    ierr = MatMult(A, x2, y2);
    CHKERRQ(ierr);
    ierr = VecView(y2, PETSC_VIEWER_STDOUT_WORLD);
    CHKERRQ(ierr);

    return PetscFinalize();
}