/*-- function C = [parmatV(V,p,K)] returns X'*V for partition matrix X. Input: V = M times N matrix p = M times 1 label vector, indicating a partition matrix X = sparse([1:M]',p,ones(M,1),M,K) K = scalor, ie. X's dimension (max label), default K = N Output: C = X' * V, K x N matrix Comments: This operation speeds up by avoiding multiplication of a 0-1 partition matrix with any matrix. It also assumes a compact representation for X. V, p, and C are all real matrices. --*/ # include "math.h" # include "mex.h" void mexFunction( int nargout, mxArray *out[], int nargin, const mxArray *in[] ) { /* declare variables */ int M,N,K,i,j,s,t,u; double *V, *p, *C; /* check argument */ if (nargin <2) { mexErrMsgTxt("Two input arguments required"); } if (nargout>1) { mexErrMsgTxt("Too many output arguments."); } /* get all inputs */ M = mxGetM(in[0]); N = mxGetN(in[0]); i = mxGetM(in[1]) * mxGetN(in[1]); if (i!=M) { mexErrMsgTxt("Dimension mismatch."); } if (nargin==3) { K = mxGetScalar(in[2]); } else { K = N; } V = mxGetPr(in[0]); p = mxGetPr(in[1]); /* computation starts */ out[0] = mxCreateDoubleMatrix(K,N,0); if (out[0]==NULL) { mexErrMsgTxt("Not enough space for the output matrix."); } C = mxGetPr(out[0]); /* multiply */ for (i=0; i