Page MenuHomec4science
Files F121743539

mexMatvec.c
No OneTemporary
Actions

File Metadata

Created
Sun, Jul 13, 13:11

mexMatvec.c
View Options

/***********************************************************************
* mexMatvec: compute
*
* mexMatvec(A,y,options)
*
* options = 0, compute A*y
* = 1, compute (y'*A)'
*
* Copyright (c) 2004 by
* K.C. Toh
* Last Modified: 120404
************************************************************************/
#include <mex.h>
#include <math.h>
/********************************
* realdotde: x dense matrix,
* y dense vector
*********************************/
double realdotde(const double *x, const int idx,
const double *y, const int n)
{ int i;
double r;
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;
}
/********************************
* saxpymat: z = z + alpha*y
* y dense matrix, z dense vector
********************************/
void saxpymat(const double *y, const int idx1,
const int istart, const int iend,
const double alp, double *z, const int idx2)
{ int i;
for(i=istart; i< iend-3; i++){ /* LEVEL 4 */
z[i+idx2] += alp * y[i+idx1]; i++;
z[i+idx2] += alp * y[i+idx1]; i++;
z[i+idx2] += alp * y[i+idx1]; i++;
z[i+idx2] += alp * y[i+idx1];
}
if(i < iend-1){ /* LEVEL 2 */
z[i+idx2] += alp * y[i+idx1]; i++;
z[i+idx2] += alp * y[i+idx1]; i++;
}
if(i < iend){ /* LEVEL 1 */
z[i+idx2] += alp * y[i+idx1];
}
return;
}
/**********************************************************/
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{ double *A, *y;
mwIndex *irA, *jcA, *iry, *jcy;
double *ytmp, *Ay;
int isspA, isspy, options;
int m1, n1, m2, n2, j, jm1;
int i, r, k, istart, iend, kstart, kend;
double tmp;
/* CHECK THE DIMENSIONS */
if (mxIsCell(prhs[0]) || mxIsCell(prhs[1])) {
mexErrMsgTxt(" mexMatvec: A, x must be a double array"); }
if (nrhs <2) {
mexErrMsgTxt(" mexMatvec: must have at least 2 inputs"); }
if (nlhs > 2) {
mexErrMsgTxt("mexMatvec: requires 1 output argument"); }
if (nrhs == 2) { options = 0; }
else { options = (int) *mxGetPr(prhs[2]); }
/***** assign pointers *****/
A = mxGetPr(prhs[0]);
m1 = mxGetM(prhs[0]);
n1 = mxGetN(prhs[0]);
isspA = mxIsSparse(prhs[0]);
if (isspA) { irA = mxGetIr(prhs[0]);
jcA = mxGetJc(prhs[0]); }
isspy = mxIsSparse(prhs[1]);
m2 = mxGetM(prhs[1]);
n2 = mxGetN(prhs[1]);
if (n2 > 1) {
mexErrMsgTxt("mexMatvec: 2ND input must be a column vector"); }
if (isspy) {
iry = mxGetIr(prhs[1]);
jcy = mxGetJc(prhs[1]);
ytmp = mxGetPr(prhs[1]);
/***** copy ytmp to y *****/
y = (double*)mxCalloc(m2,sizeof(double));
kstart = jcy[0]; kend = jcy[1];
for (k=kstart; k<kend; k++) {
r = iry[k]; y[r] = ytmp[k]; }
} else {
y = mxGetPr(prhs[1]);
}
if (options == 0 && n1 != m2) {
mexErrMsgTxt("mexMatvec: 1ST and 2ND input not compatible.");
} else if (options && m1 != m2) {
mexErrMsgTxt("mexMatvec: 1ST and 2ND input not compatible.");
}
/***** create return argument *****/
if (options==0) {
plhs[0] = mxCreateDoubleMatrix(m1,1,mxREAL);
} else {
plhs[0] = mxCreateDoubleMatrix(n1,1,mxREAL);
}
Ay = mxGetPr(plhs[0]);
/***** main body *****/
if (options==0) {
if (!isspA) {
for (j=0; j<n1; j++){
jm1 = j*m1;
tmp = y[j];
if (tmp !=0) {
saxpymat(A,jm1,0,m1,tmp,Ay,0); }
}
} else {
for (j=0; j<n1; j++){
tmp = y[j];
if (tmp != 0) {
istart = jcA[j]; iend = jcA[j+1];
for (i=istart; i<iend; i++) {
r = irA[i];
Ay[r] += tmp*A[i]; }
}
}
}
} else {
if (!isspA) {
for (j=0; j<n1; j++){
jm1 = j*m1;
Ay[j] = realdotde(A,jm1,y,m1);
}
} else {
for (j=0; j<n1; j++){
istart = jcA[j]; iend = jcA[j+1];
tmp = 0;
for (i=istart; i<iend; i++) {
r = irA[i];
tmp += y[r]*A[i]; }
Ay[j] = tmp;
}
}
}
return;
}
/**********************************************************/

Event Timeline

c4science · Help