Search code examples
c++sparse-matrixlapackblas

Is sparse BLAS not included in BLAS?


I have a working LAPACK implementation and that, as far as I read, contains BLAS.

I want to use SPARSE BLAS and as far as I understand this website, SPARSE BLAS is part of BLAS.

But when I tried to run the code below from the sparse blas manual using

g++ -o sparse.x sparse_blas_example.c -L/usr/local/lib -lblas && ./sparse_ex.x

the compiler (or linker?) asked for blas_sparse.h. When I put that file in the working directory I got:

ludi@ludi-M17xR4:~/Desktop/tests$ g++  -o sparse.x sparse_blas_example.c -L/usr/local/lib -lblas && ./sparse_ex.x
In file included from sparse_blas_example.c:3:0:
blas_sparse.h:4:23: fatal error: blas_enum.h: No such file or directory
 #include "blas_enum.h"

What must I do to use SPARSE BLAS with LAPACK? I could start moving a lot of header files to the working directory, but I gathered I already have them with lapack!

/* C example: sparse matrix/vector multiplication */

#include "blas_sparse.h"
int main()
{
const int n = 4;
const int nz = 6;
double val[] = { 1.1, 2.2, 2.4, 3.3, 4.1, 4.4 };
int indx[] = { 0, 1, 1, 2, 3, 3};
int jndx[] = { 0, 1, 4, 2, 0, 3};
double x[] = { 1.0, 1.0, 1.0, 1.0 };
double y[] = { 0.0, 0.0, 0.0, 0.0 };
blas_sparse_matrix A;
double alpha = 1.0;
int i;

/*------------------------------------*/
/* Step 1: Create Sparse BLAS Handle */
/*------------------------------------*/

A = BLAS_duscr_begin( n, n );

/*------------------------------------*/
/* Step 2: insert entries one-by-one */
/*------------------------------------*/

for (i=0; i< nz; i++)
{
BLAS_duscr_insert_entry(A, val[i], indx[i], jndx[i]);
}

/*-------------------------------------------------*/
/* Step 3: Complete construction of sparse matrix */
/*-------------------------------------------------*/
BLAS_uscr_end(A);

/*------------------------------------------------*/
/* Step 4: Compute Matrix vector product y = A*x */
/*------------------------------------------------*/

BLAS_dusmv( blas_no_trans, alpha, A, x, 1, y, 1 );

/*---------------------------------*/
/* Step 5: Release Matrix Handle */
/*---------------------------------*/

BLAS_usds(A);

/*---------------------------*/
/* Step 6: Output Solution */
/*---------------------------*/

for (i=0; i<n; i++) printf("%12.4g ",y[i]);
printf("\n");
return 0;
}

Solution

  • You quote the Blas Technical Standard, not the LAPACK reference. LAPACK does not contain routines for sparse matrices, other than handling some banded matrices. There are other implementations such as spblas and sparse which follow the techincal standard and implement sparse BLAS. Typically, sparse operations are not considered part of BLAS, but an extension.

    I would recommend using a higher level library, such as eigen because it will save you a significant amount of development time, with usually small performance costs. There is also ublas which is part of boost, so if you are using boost as part of your project, you could give it a try, though it's not really optimised for performance. You can find a comprehensive list here (again, note that LAPACK is not listed as having support for sparse operations).