/* S = GenerateTriangleExtractorMatrixMex(numRows,triangleFlag,diagFlag)
*
* S = double sparse matrix
* numRows = integer > 0
* triangleFlag = 1 , extract lower triangle
* 2 , extract upper triangle
* diagFlag = 1 , include diagonal
* 2 , exclude diagonal
* where
*
* M = an numRows X numRows matrix of non-zero terms
* assert(isequal(S * M(:), mX(logical(tril(M, -1)))));
* assert(isequal(S * M(:), mX(logical(triu(M, 1)))));
*
* Programmer: James Tursa
* Date: 2020-April-22
*/
#include "mex.h"
void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize numRows, triangleFlag, diagFlag, numElements;
mwIndex *Ir, *Jc;
mwIndex i, j, k, m;
double *pr;
if( nrhs != 3 || !mxIsNumeric(prhs[0]) || !mxIsNumeric(prhs[1]) || !mxIsNumeric(prhs[2]) ||
mxGetNumberOfElements(prhs[0]) != 1 || mxGetNumberOfElements(prhs[1]) != 1 ||
mxGetNumberOfElements(prhs[2]) != 1 ) {
mexErrMsgTxt("Need three numeric scalar inputs");
}
if( nlhs > 1 ) {
mexErrMsgTxt("Too many outputs");
}
numRows = mxGetScalar(prhs[0]);
triangleFlag = mxGetScalar(prhs[1]);
diagFlag = mxGetScalar(prhs[2]);
if( numRows < 1 ) {
mexErrMsgTxt("Invalid numRows, should be > 0");
}
if( triangleFlag != 1 && triangleFlag != 2 ) {
mexErrMsgTxt("Invalid triangleFlag, should be 1 or 2");
}
if( diagFlag != 1 && diagFlag != 2 ) {
mexErrMsgTxt("Invalid diagFlag, should be 1 or 2");
}
if( diagFlag == 1 ) {
numElements = numRows * (numRows + 1) / 2; /* include diagonal */
} else {
numElements = (numRows - 1) * numRows / 2; /* exclude diagonal */
}
plhs[0] = mxCreateSparse(numElements, numRows*numRows, numElements, mxREAL);
pr = (double *) mxGetData(plhs[0]);
Ir = mxGetIr(plhs[0]);
Jc = mxGetJc(plhs[0]);
Jc[0] = 0;
diagFlag--;
k = 0;
m = 1;
if( triangleFlag == 1 ) { /* Lower */
for( j=0; j<numRows; j++ ) {
for( i=0; i<numRows; i++ ) {
if( i >= j+diagFlag ) {
*pr++ = 1.0;
*Ir++ = k++;
Jc[m] = Jc[m-1] + 1;
} else {
Jc[m] = Jc[m-1];
}
m++;
}
}
} else { /* Upper */
for( j=0; j<numRows; j++ ) {
for( i=0; i<numRows; i++ ) {
if( i+diagFlag <= j ) {
*pr++ = 1.0;
*Ir++ = k++;
Jc[m] = Jc[m-1] + 1;
} else {
Jc[m] = Jc[m-1];
}
m++;
}
}
}
}
n = 300;
disp('m-code timing')
tic
GenerateTriangleExtractorMatrix(10000,1,1);
toc
disp('mex code timing')
tic
GenerateTriangleExtractorMatrixMex(10000,1,1);
toc
for k=1:n
numRows = ceil(rand*5000+100);
numCols = numRows;
triangleFlag = (rand<0.5) + 1;
diagFlag = (rand<0.5) + 1;
Mm = GenerateTriangleExtractorMatrix(numRows,triangleFlag,diagFlag);
Mx = GenerateTriangleExtractorMatrixMex(numRows,triangleFlag,diagFlag);
if( ~isequal(Mm,Mx) )
error('Not equal');
end
end
disp('Random tests passed')