Page MenuHomec4science
Files F122202172

mexinprod.c
No OneTemporary
Actions

File Metadata

Created
Wed, Jul 16, 14:15

mexinprod.c
View Options

/***********************************************************************
* mexinprod.c : C mex file to compute inner product of 2 vectors,
*
* Ax = mexinprod(blk,Avec,x,k,p);
* where Ax = Avec{p}(:,1:k)'*x
*
* 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>
/**********************************************************
* realdot1: x dense, y dense
**********************************************************/
double realdot1(const double *x, const int col,
const double *y, const int n)
{ int i, idx;
double r;
idx=col*n;
r=0.0;
for (i=0; i<n-3; i++) { /* LEVEL 4 */
r += x[i+idx] * y[i]; i++;
r += x[i+idx] * y[i]; i++;
r += x[i+idx] * y[i]; i++;
r += x[i+idx] * y[i]; }
if (i<n-1) { /* LEVEL 2 */
r += x[i+idx] * y[i]; i++;
r += x[i+idx] * y[i]; i++; }
if (i<n) { /* LEVEL 1 */
r += x[i+idx] * y[i]; }
return r;
}
/**********************************************************
* realdot2: x sparse, y dense
**********************************************************/
double realdot2(const double *x, mwIndex *irx, mwIndex *jcx,
const int col, const double *y)
{ int i, row, idxstart, idxend;
double r;
idxstart=jcx[col]; idxend=jcx[col+1];
r=0.0;
for (i=idxstart; i<idxend-3; i++) { /* LEVEL 4 */
row = irx[i]; r += x[i]*y[row]; i++;
row = irx[i]; r += x[i]*y[row]; i++;
row = irx[i]; r += x[i]*y[row]; i++;
row = irx[i]; r += x[i]*y[row]; }
if (i<idxend-1) { /* LEVEL 2 */
row = irx[i]; r += x[i]*y[row]; i++;
row = irx[i]; r += x[i]*y[row]; i++; }
if (i<idxend) { /* LEVEL 1 */
row = irx[i]; r += x[i]*y[row]; }
return r;
}
/**********************************************************/
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{ mxArray *A_cell_pr;
double *A, *B;
mwIndex *irA, *jcA, *irB, *jcB;
double *Btmp, *tr;
int isspA, isspB, iscellA, iscellB;
mwIndex subs[2];
mwSize nsubs=2;
int mA, nA, m1, n1, m2, n2, j, index;
int rowidx, colidx, r, k, kstart, kend;
/* CHECK THE DIMENSIONS */
iscellA = mxIsCell(prhs[1]);
mA = mxGetM(prhs[1]);
if (!iscellA) { mA = 1; }
if (nrhs < 2) {
mexErrMsgTxt(" mexinprod: must have at least 3 inputs"); }
if (nlhs>2) {
mexErrMsgTxt("mexinprod: requires 1 output argument"); }
if (nrhs > 4) { rowidx = (int)*mxGetPr(prhs[4]); } else { rowidx = 1; }
if (rowidx > mA) {
mexErrMsgTxt("mexinprod: rowidx exceeds size(Avec,1)");
}
/***** assign pointers *****/
iscellB = mxIsCell(prhs[2]);
isspB = mxIsSparse(prhs[2]);
m2 = mxGetM(prhs[2]);
n2 = mxGetN(prhs[2]);
if ((n2 > 1) || (iscellB)) {
mexErrMsgTxt("mexinprod: 3RD input must be a column vector"); }
if (isspB) {
irB = mxGetIr(prhs[2]);
jcB = mxGetJc(prhs[2]);
Btmp = mxGetPr(prhs[2]);
/***** copy Btmp to B *****/
B = (double*)mxCalloc(m2,sizeof(double));
kstart = jcB[0]; kend = jcB[1];
for (k=kstart; k<kend; k++) {
r = irB[k]; B[r] = Btmp[k]; }
} else {
B = mxGetPr(prhs[2]);
}
if (iscellA) {
subs[0] = rowidx-1; /* subtract 1 to adjust for Matlab index */
subs[1] = 0;
index = mxCalcSingleSubscript(prhs[1],nsubs,subs);
A_cell_pr = mxGetCell(prhs[1],index);
A = mxGetPr(A_cell_pr);
m1 = mxGetM(A_cell_pr);
n1 = mxGetN(A_cell_pr);
isspA = mxIsSparse(A_cell_pr);
if (isspA) { irA = mxGetIr(A_cell_pr);
jcA = mxGetJc(A_cell_pr); }
} else {
A = mxGetPr(prhs[1]);
m1 = mxGetM(prhs[1]);
n1 = mxGetN(prhs[1]);
isspA = mxIsSparse(prhs[1]);
if (isspA) { irA = mxGetIr(prhs[1]);
jcA = mxGetJc(prhs[1]); }
}
if (nrhs > 3) { colidx = (int)*mxGetPr(prhs[3]); } else { colidx = 1; }
if (colidx > n1) {
mexErrMsgTxt("mexinprod: colidx exceeds size(Avec,2)"); }
if (m1 != m2) {
mexErrMsgTxt("mexinprod: 2ND and 3RD input not compatible.");
}
/***** create return argument *****/
plhs[0] = mxCreateDoubleMatrix(colidx,1,mxREAL);
tr = mxGetPr(plhs[0]);
/***** compute <Aj,B> *****/
if (isspA) {
for (j=0; j<colidx; j++){
tr[j] = realdot2(A,irA,jcA,j,B); }
} else {
for (j=0; j<colidx; j++){
tr[j] = realdot1(A,j,B,m1); }
}
if (isspB) { mxFree(B); }
return;
}
/**********************************************************/

Event Timeline

c4science · Help