idamnjanovic@70: /**************************************************************************
idamnjanovic@70:  *
idamnjanovic@70:  * File name: rowlincomb.c
idamnjanovic@70:  *
idamnjanovic@70:  * Ron Rubinstein
idamnjanovic@70:  * Computer Science Department
idamnjanovic@70:  * Technion, Haifa 32000 Israel
idamnjanovic@70:  * ronrubin@cs
idamnjanovic@70:  *
idamnjanovic@70:  * Last Updated: 21.5.2009
idamnjanovic@70:  *
idamnjanovic@70:  *************************************************************************/
idamnjanovic@70: 
idamnjanovic@70: #include "mex.h"
idamnjanovic@70: 
idamnjanovic@70: 
idamnjanovic@70: /* Input Arguments */
idamnjanovic@70: 
idamnjanovic@70: #define	X_IN	   prhs[0]
idamnjanovic@70: #define A_IN     prhs[1]
idamnjanovic@70: #define ROWS_IN  prhs[2]
idamnjanovic@70: #define COLS_IN  prhs[3]
idamnjanovic@70: 
idamnjanovic@70: 
idamnjanovic@70: /* Output Arguments */
idamnjanovic@70: 
idamnjanovic@70: #define	Y_OUT	plhs[0]
idamnjanovic@70: 
idamnjanovic@70: 
idamnjanovic@70: void mexFunction(int nlhs, mxArray *plhs[],
idamnjanovic@70: int nrhs, const mxArray*prhs[])
idamnjanovic@70: 
idamnjanovic@70: {
idamnjanovic@70:   double *A, *x, *y, *rows, *cols;
idamnjanovic@70:   mwSize m,n,m1,n1,m2,n2,rownum,colnum;
idamnjanovic@70:   mwIndex *row_ids,*col_ids,i,j;
idamnjanovic@70:   int colnumspecified=0;
idamnjanovic@70:   
idamnjanovic@70:   
idamnjanovic@70:   /* Check for proper number of arguments */
idamnjanovic@70:   
idamnjanovic@70:   if (nrhs!=3 && nrhs!=4) {
idamnjanovic@70:     mexErrMsgTxt("Invalid number of input arguments.");
idamnjanovic@70:   } else if (nlhs > 1) {
idamnjanovic@70:     mexErrMsgTxt("Too many output arguments.");
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   
idamnjanovic@70:   /* Check the input dimensions */
idamnjanovic@70:   
idamnjanovic@70:   m = mxGetM(A_IN);
idamnjanovic@70:   n = mxGetN(A_IN);
idamnjanovic@70:   if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) {
idamnjanovic@70:     mexErrMsgTxt("ROWLINCOMB requires that A be a double matrix.");
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   m1 = mxGetM(ROWS_IN);
idamnjanovic@70:   n1 = mxGetN(ROWS_IN);
idamnjanovic@70:   if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) {
idamnjanovic@70:     mexErrMsgTxt("ROWLINCOMB requires that ROWS be an index vector of type double.");
idamnjanovic@70:   }
idamnjanovic@70:   rownum = (m1 > n1) ? m1 : n1;   /* the number of rows in the linear combination */
idamnjanovic@70:   
idamnjanovic@70:   m2 = mxGetM(X_IN);
idamnjanovic@70:   n2 = mxGetN(X_IN);
idamnjanovic@70:   if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ((m2!=1) && (n2!=1))) {
idamnjanovic@70:     mexErrMsgTxt("ROWLINCOMB requires that X be a double vector.");
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   if (m2 != rownum && n2 != rownum) {
idamnjanovic@70:     mexErrMsgTxt("The length of X does not match the number of rows in ROWS.");
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   if (nrhs==4) {
idamnjanovic@70:     m1 = mxGetM(COLS_IN);
idamnjanovic@70:     n1 = mxGetN(COLS_IN);
idamnjanovic@70:     if (!mxIsDouble(COLS_IN) || mxIsComplex(COLS_IN) || (m1!=1 && n1!=1)) {
idamnjanovic@70:       mexErrMsgTxt("ROWLINCOMB requires that COLS be an index vector of type double.");
idamnjanovic@70:     }
idamnjanovic@70:     colnum = (m1 > n1) ? m1 : n1;   /* the number of columns */
idamnjanovic@70:     colnumspecified = 1;
idamnjanovic@70:     cols = mxGetPr(COLS_IN);
idamnjanovic@70:     
idamnjanovic@70:     Y_OUT = mxCreateDoubleMatrix(1, colnum, mxREAL);
idamnjanovic@70:   }
idamnjanovic@70:   else {
idamnjanovic@70:     cols = 0;
idamnjanovic@70:     Y_OUT = mxCreateDoubleMatrix(1, n, mxREAL);
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   
idamnjanovic@70:   /* Assign pointers to the various parameters */
idamnjanovic@70:   A = mxGetPr(A_IN);
idamnjanovic@70:   rows = mxGetPr(ROWS_IN);
idamnjanovic@70:   x = mxGetPr(X_IN);
idamnjanovic@70:   y = mxGetPr(Y_OUT);
idamnjanovic@70:   
idamnjanovic@70:   
idamnjanovic@70:   /* check row indices */
idamnjanovic@70:   
idamnjanovic@70:   row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex));
idamnjanovic@70:   
idamnjanovic@70:   for (i=0; i<rownum; ++i) {
idamnjanovic@70:     row_ids[i] = (mwIndex)(rows[i]+0.1)-1;
idamnjanovic@70:     if (row_ids[i]<0 || row_ids[i]>=m) {
idamnjanovic@70:       mexErrMsgTxt("Row index in ROWS is out of range.");
idamnjanovic@70:     }
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   
idamnjanovic@70:   
idamnjanovic@70:   if (colnumspecified) {
idamnjanovic@70:     
idamnjanovic@70:     /* check column indices */
idamnjanovic@70:     col_ids = (mwIndex*)mxMalloc(colnum*sizeof(mwIndex));
idamnjanovic@70:     
idamnjanovic@70:     for (i=0; i<colnum; ++i) {
idamnjanovic@70:       col_ids[i] = (mwIndex)(cols[i]+0.1)-1;
idamnjanovic@70:       if (col_ids[i]<0 || col_ids[i]>=n) {
idamnjanovic@70:         mexErrMsgTxt("Column index in COLS is out of range.");
idamnjanovic@70:       }
idamnjanovic@70:     }
idamnjanovic@70:     
idamnjanovic@70:     /* Do the actual computation */
idamnjanovic@70:     for (j=0; j<colnum; ++j) {
idamnjanovic@70:       for (i=0; i<rownum; ++i) {
idamnjanovic@70:         y[j] += A[m*col_ids[j]+row_ids[i]]*x[i];
idamnjanovic@70:       }
idamnjanovic@70:     }
idamnjanovic@70:     
idamnjanovic@70:     mxFree(col_ids);
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   else {
idamnjanovic@70:     
idamnjanovic@70:     /* Do the actual computation */
idamnjanovic@70:     for (j=0; j<n; ++j) {
idamnjanovic@70:       for (i=0; i<rownum; ++i) {
idamnjanovic@70:         y[j] += A[m*j+row_ids[i]]*x[i];
idamnjanovic@70:       }
idamnjanovic@70:     }
idamnjanovic@70:   }
idamnjanovic@70:   
idamnjanovic@70:   
idamnjanovic@70:   mxFree(row_ids);
idamnjanovic@70:   
idamnjanovic@70:   return;
idamnjanovic@70: }