A sparse matrix vector multiplication in C mex programming
显示 更早的评论
Hi everyone,
I'm new in C mex programming and I'm really stuck. I would like to write a code that computes a sparse matrix vector multiplication (sparse_matrix*vector=result_vector). Given how sparse matrices are managed in C mex, I have difficulties getting the non zero elements of each row (to do the multiplication with the vector), because the pointers Ir=mxGetIr(sparse_matrix) and Jc=mxGetJc(sparse_matrix) give the information about the columns only and not the rows.
I would be really grateful if somebody helps me out because I can't move on without solving this problem.
Thanks in advance.
Sincerely,
J
3 个评论
Matt J
2017-11-28
Why must you re-implement an operation that MATLAB has already implemented for you?
In any case, why must you implement sparse_matrix(i,:)*vector? Why can't you implement sum_i sparse_matrix(:,i)*vector(i)?
James Tursa
2017-11-28
I'm with Matt J. Why are you doing this? And if for some reason you really do need to do this, what have you done so far? Where specifically are you stuck? The mxGetIr( ) gives the row indexing, btw. Is the vector full or sparse?
JM
2017-11-28
采纳的回答
James Tursa
2017-11-29
编辑：James Tursa
2017-12-6
Here is the bare bones C code for this. Most of the code consists of argument checking etc. The actual computation part of the code is only a few lines and is fairly trivial. It simply accumulates the individual element multiplies into the appropriate resulting element. The j for-loop loops over the columns of M, and the while loop does the individual element multiplies within each column.
/* Bare bones sparse matrix * full vector */
/* Produces a full vector result */
/* Assumes all input elements are finite */
/* E.g., does not account for NaN or 0*inf etc. */
/* Programmer: James Tursa */
/* Includes ----------------------------------------------------------- */
#include "mex.h"
/* Gateway ------------------------------------------------------------ */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize m, n, j, nrow;
double *Mpr, *Vpr, *Cpr;
mwIndex *Mir, *Mjc;
/* Argument checks */
if( nlhs > 1 ) {
mexErrMsgTxt("Too many outputs");
}
if( nrhs != 2 ) {
mexErrMsgTxt("Need exactly two inputs");
}
if( !mxIsDouble(prhs[0]) || !mxIsSparse(prhs[0]) || mxIsComplex(prhs[0]) ) {
mexErrMsgTxt("1st argument must be real double sparse matrix");
}
if( !mxIsDouble(prhs[1]) || mxIsSparse(prhs[1]) || mxIsComplex(prhs[1]) ||
mxGetNumberOfDimensions(prhs[1])!=2 || mxGetN(prhs[1])!=1 ) {
mexErrMsgTxt("2nd argument must be real double full column vector");
}
m = mxGetM(prhs[0]);
n = mxGetN(prhs[0]);
if( mxGetM(prhs[1])!=n ) {
mexErrMsgTxt("Inner matrix dimensions must agree.");
}
/* Create output */
plhs[0] = mxCreateDoubleMatrix( m, 1, mxREAL );
/* Get data pointers */
Mpr = mxGetPr(prhs[0]);
Vpr = mxGetPr(prhs[1]);
Cpr = mxGetPr(plhs[0]);
Mir = mxGetIr(prhs[0]);
Mjc = mxGetJc(prhs[0]);
/* Calculate result */
for( j=0; j<n; j++ ) {
nrow = Mjc[j+1] - Mjc[j]; /* Number of row elements for this column */
while( nrow-- ) {
Cpr[*Mir++] += *Mpr++ * Vpr[j]; /* Accumulate contribution of Vpr[j] */
}
}
}
2 个评论
JM
2017-12-1
Royi Avital
2021-8-26
Isn't there an internal function to do this efficiently (Multi Threaded)?
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
