// Extract the elements in X (a dense matrix) indexed by jc and ir.
// Usage: b = proj_omega_sub(X,jc,ir);

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){

  mwSize m = mxGetM(prhs[0]);
  mwSize n = mxGetN(prhs[0]);
  
  double* X = mxGetPr(prhs[0]);

  mwIndex* jc = (mwIndex*)mxGetData(prhs[1]);
  mwIndex* ir = (mwIndex*)mxGetData(prhs[2]);
  int nb_samples = mxGetNumberOfElements(prhs[2]);

  plhs[0] = mxCreateDoubleMatrix(nb_samples, 1, mxREAL);
  double *b = mxGetPr(plhs[0]), *p;

  mwSize nnz;
  for(mwIndex j = 0; j < n; j++) 
  { 
    nnz = jc[j+1] - jc[j];
    p = X + j*m;
    for (mwIndex i = 0; i < nnz; i++)
      *b++ = *(p + *(ir++));
  }
}