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