Symmetric matrix-vector multiplication with Accelerate Sparse BLAS

48 views Asked by At

I'd like to multiply a vector with a symmetric sparse matrix using Apple's Accelerate framework that also includes facilities for sparse matrices. Despite setting the symmetric property for the matrix, only the upper triangle is multiplied in contrast to my expectations. Unfortunately the documentation is sparse (pun intended) on what the expected behaviour is supposed to be. For example the Intel MKL library will multiply the whole matrix, and not just half of it.

Here is the minimal example and output when using the Accelerate framework:

// spmv.c
// compiled with: clang -Wall spmv.c -framework Accelerate

#include <stdio.h>
#include <Accelerate/Accelerate.h>

int main(void) {

    sparse_dimension n = 4;
    sparse_dimension nnz = 8;

    double va[] = {1.0, 2.0, 3.0, 4.0,
                        5.0,      6.0,
                             7.0,     
                                  8.0};
    sparse_index ia[] = {0,0,0,0,1,1,2,3};
    sparse_index ja[] = {0,1,2,3,1,3,2,3};

    sparse_matrix_double A = sparse_matrix_create_double(n, n);
    sparse_set_matrix_property(A,SPARSE_UPPER_SYMMETRIC);
    sparse_insert_entries_double(A,nnz,va,ia,ja);
    sparse_commit(A);

    double x[] = {1.0,1.0,1.0,1.0};
    double y[] = {0.0,0.0,0.0,0.0};

    // y = triu(A) x (however (A x) was expected)
    sparse_matrix_vector_product_dense_double(CblasNoTrans,1.0,A,x,1,y,1);

    for(int i = 0; i < n; ++i)
        printf("y[%d] = %f\n",i,y[i]);

    sparse_matrix_destroy(A);
}

Output using Accelerate:

$ clang -Wall spmv.c -framework Accelerate
$ ./a.out
y[0] = 10.000000
y[1] = 11.000000
y[2] = 7.000000
y[3] = 8.000000

I was expecting Accelerate would behave similar to Intel MKL, and multiply with the full matrix A, and not just the upper half:

// spmv_mkl.c
// compile with: icc -std=c99 -qmkl spmv_mkl.c

#include <stdio.h>
#include <mkl_spblas.h>

int main(void) {

    MKL_INT n = 4;
    MKL_INT nnz = 8;

    // symmetric matrix
    double va[] = {1.0, 2.0, 3.0, 4.0,
                        5.0,      6.0,
                             7.0,     
                                  8.0};

    MKL_INT ia[] = {0,0,0,0,1,1,2,3};
    MKL_INT ja[] = {0,1,2,3,1,3,2,3};

    sparse_matrix_t A;
    mkl_sparse_d_create_coo(&A,SPARSE_INDEX_BASE_ZERO,n,n,nnz,ia,ja,va);

    double x[] = {1.0,1.0,1.0,1.0};
    double y[] = {0.0,0.0,0.0,0.0};

    struct matrix_descr dA;
    dA.type = SPARSE_MATRIX_TYPE_SYMMETRIC;
    dA.mode = SPARSE_FILL_MODE_UPPER;
    dA.diag = SPARSE_DIAG_NON_UNIT;

    // y = A x
    const double alpha = 1.0;
    const double beta  = 0.0;
    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE,alpha,A,dA,x,beta,y);

    for(int i = 0; i < n; ++i)
        printf("y[%d] = %f\n",i,y[i]);

    mkl_sparse_destroy(A);
}

Output using MKL:

$ icc -Wall -std=c99 -qmkl -diag-disable=10441 spmv_mkl.c
$ ./a.out
y[0] = 10.000000
y[1] = 13.000000
y[2] = 10.000000
y[3] = 18.000000
0

There are 0 answers