Page MenuHomec4science

mexqops.c
No OneTemporary

File Metadata

Created
Thu, Jul 10, 01:40

mexqops.c

/***********************************************************************
* mexqops.c : C mex file
*
* z = mexqops(blk,x,y,options);
*
* Input: blk = [n1, n2, ... nk]
* x = n-vector or k-vector
* y = n-vector, where n = n1+...+nk
*
* options = 1, z = k-vector, z(i) = <xi,yi>
* = 2, z = k-vector, z(i) = 2*xi(1)yi(1) - <xi,yi>
* = 3, z = n-vector, zi = x(i)*yi
* = 4, z = n-vector, zi = x(i)*yi, zi(1) = -zi(1).
*
* SDPT3: version 3.0
* Copyright (c) 1997 by
* K.C. Toh, M.J. Todd, R.H. Tutuncu
* Last Modified: 2 Feb 01
***********************************************************************/
#include <mex.h>
#include <math.h>
/**********************************************************
* ops1
**********************************************************/
void ops1(double *x, double *y, double *z,
int numblk, int *cumblk, int options)
{ int j, l, jstart, jend;
double tmp;
if (options == 1) {
for (l=0; l<numblk; ++l) {
jstart = cumblk[l]; jend = cumblk[l+1];
tmp = 0;
for (j=jstart; j<jend; ++j) { tmp += x[j]*y[j]; }
z[l] = tmp;
}
}
else if (options == 2) {
for (l=0; l<numblk; ++l) {
jstart = cumblk[l]; jend = cumblk[l+1];
tmp = x[jstart]*y[jstart];
for (j=jstart+1; j<jend; ++j) { tmp -= x[j]*y[j]; }
z[l] = tmp;
}
}
return;
}
/**********************************************************
* ops3
**********************************************************/
void ops3(double *x, double *y, double *z,
int numblk, int *cumblk, int options)
{ int j, l, jstart, jend;
double tmp;
if (options == 3) {
for (l=0; l<numblk; ++l) {
jstart = cumblk[l]; jend = cumblk[l+1];
tmp = x[l];
for (j=jstart; j<jend; ++j) {
z[j] = tmp*y[j]; }
}
}
else if (options == 4) {
for (l=0; l<numblk; ++l) {
jstart = cumblk[l]; jend = cumblk[l+1];
tmp = x[l];
for (j=jstart; j<jend; ++j) {
z[j] = tmp*y[j]; }
z[jstart] = -z[jstart];
}
}
return;
}
/**********************************************************
*
***********************************************************/
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{ double *blksize, *x, *y, *z, *xtmp, *ytmp;
mwIndex *irx, *jcx, *iry, *jcy;
int *cumblk;
int mblk, options;
int l, n, k, r, numblk, cols;
/* CHECK FOR PROPER NUMBER OF ARGUMENTS */
if (nrhs !=4){
mexErrMsgTxt("mexqops: requires 4 input arguments."); }
if (nlhs > 1){
mexErrMsgTxt("mexqops: requires 1 output argument."); }
/* CHECK THE DIMENSIONS */
numblk = mxGetN(prhs[0]);
blksize = mxGetPr(prhs[0]);
cumblk = (int*)mxCalloc(numblk+1,sizeof(int));
for (l=0; l<numblk; l++) {
cols = (int)blksize[l];
cumblk[l+1] = cumblk[l] + cols;
}
n = cumblk[numblk];
/***** assign pointers *****/
options = (int)*mxGetPr(prhs[3]);
if (mxGetM(prhs[2]) != n) {
mexErrMsgTxt("mexqops: dim not compatible."); }
if (options < 3 && mxGetM(prhs[1]) != n) {
mexErrMsgTxt("mexqops: dim not compatible.."); }
if (options >= 3 && mxGetM(prhs[1]) != numblk) {
mexErrMsgTxt("mexqops: dim not compatible..."); }
if (mxIsSparse(prhs[1])) {
irx = mxGetIr(prhs[1]); jcx = mxGetJc(prhs[1]); xtmp = mxGetPr(prhs[1]);
x = (double*)mxCalloc(n,sizeof(double));
for (k=0; k<jcx[1]; ++k) { r=irx[k]; x[r]=xtmp[k]; } }
else {
x = mxGetPr(prhs[1]);
}
if (mxIsSparse(prhs[2])) {
iry = mxGetIr(prhs[2]); jcy = mxGetJc(prhs[2]); ytmp = mxGetPr(prhs[2]);
y = (double*)mxCalloc(n,sizeof(double));
for (k=0; k<jcy[1]; ++k) { r=iry[k]; y[r]=ytmp[k]; } }
else {
y = mxGetPr(prhs[2]);
}
/***** Do the computations in a subroutine *****/
if (options < 3) {
plhs[0] = mxCreateDoubleMatrix(numblk,1,mxREAL);
z = mxGetPr(plhs[0]);
ops1(x,y,z,numblk,cumblk,options);
}
else if (options >= 3) {
plhs[0] = mxCreateDoubleMatrix(n,1,mxREAL);
z = mxGetPr(plhs[0]);
ops3(x,y,z,numblk,cumblk,options); }
return;
}
/**********************************************************/

Event Timeline