dimanche 26 juin 2016

sparse BLAS solver used in a shared library doesn't work (returns '-1')


I have tried to program a solver-function in C (with Ubuntu in Code::Blocks). I want to call that solver from a Python script, so I have build it as a shared library and called it with 'ctypes' in Python. At first have tried the example from there http://www.netlib.org/blas/blast-forum/chapter3.pdf (page 12). This is a simple matrix(sparse)-vector(dense) multiplication. And everything is working just fine. Here the multiplication code:

const int n =4;     //dimension of sparse matrix
const int nz=6;     //amount of non-zero elements
double val[] = {1.1,2.2,3.3,4.1,2.4,4.4,20.0};  //value
int indx[]={0,1,2,3,3,3,3};     //row index
int jndx[]={0,1,2,0,1,3,3};     //column index

double x[]={1.0,1.0,1.0,1.0};   //x-vector
double y[]={0.0,0.0,0.0,0.0};   //output vector

blas_sparse_matrix A; 
double alph=1.0;


A=BLAS_duscr_begin(n,n);    // give matrix A the dimension
for(i=0;i<nz;i++){
    BLAS_duscr_insert_entry(A,val[i],indx[i],jndx[i]);  //write entries in A
}
//BLAS_ussp(A,blas_lower_symmetric);  
//BLAS_ussp(A,blas_lower_triangular);
BLAS_duscr_end(A);

BLAS_dusmv(blas_no_trans, alph,A,x,1,y,1); //BLAS_dussv(blas_no_trans,alph,A,x,1);
BLAS_usds(A);

for(i=0;i<n;i++){
    printf("%fn",y[i]);
}

It performs y<-A*x and returns 0 when everything was successful.

Now I have just changed this line:

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

To:

BLAS_dussv(blas_no_trans,alph,A,x,1);

This is a triangular solver which performs x <- A^(-1)*x. But this function just returns -1 (something went wrong) and doesn't calculate anything.

Does someone know how to use the solver of sparse BLAS correctly?

And does the solver just solve a triangular matrix or does it interpret it as a square matrix? (takes the lower half of the matrix and assumes it as a symmetric matrix)


Aucun commentaire:

Enregistrer un commentaire