Sparse Matrix Multiplication

Hi,

I know there are many mex based routines are available for sparse matrix multiplication (such as SSMULT), however, I am writing my own as a part of mex exercise. Here is the code that I have written.

#include "mex.h"
int scatter(mxArray *A,int j,double beta,int *w,double *x,int mark,mxArray *C,int nz)
{
  mwIndex *Ap,*Ai,*Ci;
  int i,p;
  double *Ax = mxGetPr(A);
  Ap = mxGetJc(A);
  Ai = mxGetIr(A);
  Ci = mxGetIr(C);
  for(p = Ap[j]; p < Ap[j+1];p++)
  {
    i = Ap[p];
    if(w[i] < mark)
    {
      w[i] = mark;
      Ci[nz++] = i;
      if(x)
      {
         x[i] = beta * Ax[p];
      }
    }  
    else if(x)
    {
       x[i] += beta * Ax[p];  
    }
  }    
  return nz;
}
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
  
  double *A,*B;
  mxArray *temp,*C;
  
  A = mxGetPr(prhs[0]);
  B = mxGetPr(prhs[1]);

  int i,j,nz = 0;

  int m_a = mxGetM(prhs[0]);
  int n_a = mxGetN(prhs[0]);
  int m_b = mxGetM(prhs[1]);
  int n_b = mxGetN(prhs[1]);
      
  mwIndex *jc_a,*ir_a,*jc_b,*ir_b,*ir,*jc;

  jc_a = mxGetJc(prhs[0]);
  ir_a = mxGetIr(prhs[0]);
  
  jc_b = mxGetJc(prhs[1]);
  ir_b = mxGetIr(prhs[1]);
  
  int m = m_a;
  int anz = jc[n_a];
  int values = (mxGetPr(prhs[0]) != NULL) && (mxGetPr(prhs[1]) !=NULL);
  int n = n_b;
  int bnz = jc_b[n];
  
  int *w = (int*) mxMalloc(m*sizeof(int));
  double *x = (double *) mxMalloc(m*sizeof(double));

  double *Cx = mxGetPr(C);
  
  temp = mxCreateSparse(m_a,n_a,anz,mxREAL);
  mxSetPr(temp,A);

  C = mxCreateSparse(m,n,anz+bnz,mxREAL);
  jc = mxGetJc(C);
  ir = mxGetIr(C);
  
  for(j = 0;j < n;j++)
  {
    jc[j] = nz;
    for(i = ir_b[j]; i < ir_b [j+1]; j++)
    {
      nz = scatter(temp,ir_b[i],B ? B[i] : 1,w,x,j+1,C,nz);      
    }    
    if(values)
    {
       for(i = jc[j]; i < nz;i++)
       {
         Cx[i] = x[ir[i]];
       }
    }
  }
  jc[n] = nz;

  mwSize ni, nrow;
        int y,z;
        double *pr;

  if( !mxIsSparse(C) ) 
  {
    mexPrintf("Error in calculation");
    return;
        }
  else
  {
   ni = mxGetN(C);
         pr = mxGetPr(C);

        for( y=0; y<n; y++ ) {
            nrow = jc[y+1] - jc[y];
            for( z=0; z<nrow; z++ ) {
          mexPrintf("   (%d,%d)    %g\n",(*ir++)+1,y+1,*pr++);
            }
        }
   }
  mxFree(x);
  mxFree(w);
  }

However, whenever I run this code by giving two sparse matrices as input, it just crashes (and possibly gives segmentation fault).

I am using cs_mult and scatter methods of SSMULT library as my baseline algorithms.

Kindly help guys. Cheers.