annotate Problems/private/collincomb.c @ 61:42fcbcfca132

(none)
author idamnjanovic
date Tue, 15 Mar 2011 12:21:31 +0000
parents 207a6ae9a76f
children
rev   line source
idamnjanovic@10 1 /**************************************************************************
idamnjanovic@10 2 *
idamnjanovic@10 3 * File name: collincomb.c
idamnjanovic@10 4 *
idamnjanovic@10 5 * Ron Rubinstein
idamnjanovic@10 6 * Computer Science Department
idamnjanovic@10 7 * Technion, Haifa 32000 Israel
idamnjanovic@10 8 * ronrubin@cs
idamnjanovic@10 9 *
idamnjanovic@10 10 * Last Updated: 21.5.2009
idamnjanovic@10 11 *
idamnjanovic@10 12 *************************************************************************/
idamnjanovic@10 13
idamnjanovic@10 14
idamnjanovic@10 15 #include "mex.h"
idamnjanovic@10 16
idamnjanovic@10 17
idamnjanovic@10 18 /* Input Arguments */
idamnjanovic@10 19
idamnjanovic@10 20 #define A_IN prhs[0]
idamnjanovic@10 21 #define ROWS_IN prhs[1]
idamnjanovic@10 22 #define COLS_IN1 prhs[1]
idamnjanovic@10 23 #define COLS_IN2 prhs[2]
idamnjanovic@10 24 #define X_IN1 prhs[2]
idamnjanovic@10 25 #define X_IN2 prhs[3]
idamnjanovic@10 26
idamnjanovic@10 27
idamnjanovic@10 28 /* Output Arguments */
idamnjanovic@10 29
idamnjanovic@10 30 #define Y_OUT plhs[0]
idamnjanovic@10 31
idamnjanovic@10 32
idamnjanovic@10 33 void mexFunction(int nlhs, mxArray *plhs[],
idamnjanovic@10 34 int nrhs, const mxArray*prhs[])
idamnjanovic@10 35
idamnjanovic@10 36 {
idamnjanovic@10 37 double *A, *x, *y, *rows, *cols;
idamnjanovic@10 38 mwSize m,n,m1,n1,m2,n2,rownum,colnum;
idamnjanovic@10 39 mwIndex col_id,*row_ids,i,j;
idamnjanovic@10 40 int rownumspecified=0;
idamnjanovic@10 41
idamnjanovic@10 42
idamnjanovic@10 43 /* Check for proper number of arguments */
idamnjanovic@10 44
idamnjanovic@10 45 if (nrhs!=3 && nrhs!=4) {
idamnjanovic@10 46 mexErrMsgTxt("Invalid number of arguments.");
idamnjanovic@10 47 } else if (nlhs > 1) {
idamnjanovic@10 48 mexErrMsgTxt("Too many output arguments.");
idamnjanovic@10 49 }
idamnjanovic@10 50
idamnjanovic@10 51
idamnjanovic@10 52 /* Check the input dimensions */
idamnjanovic@10 53
idamnjanovic@10 54 m = mxGetM(A_IN);
idamnjanovic@10 55 n = mxGetN(A_IN);
idamnjanovic@10 56 if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) {
idamnjanovic@10 57 mexErrMsgTxt("COLLINCOMB requires that A be a double matrix.");
idamnjanovic@10 58 }
idamnjanovic@10 59
idamnjanovic@10 60 if (nrhs==3) {
idamnjanovic@10 61
idamnjanovic@10 62 m1 = mxGetM(COLS_IN1);
idamnjanovic@10 63 n1 = mxGetN(COLS_IN1);
idamnjanovic@10 64 if (!mxIsDouble(COLS_IN1) || mxIsComplex(COLS_IN1) || (m1!=1 && n1!=1)) {
idamnjanovic@10 65 mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double.");
idamnjanovic@10 66 }
idamnjanovic@10 67 colnum = (m1 > n1) ? m1 : n1; /* the number of columns in the linear combination */
idamnjanovic@10 68
idamnjanovic@10 69 m2 = mxGetM(X_IN1);
idamnjanovic@10 70 n2 = mxGetN(X_IN1);
idamnjanovic@10 71 if (!mxIsDouble(X_IN1) || mxIsComplex(X_IN1) || (m2!=1 && n2!=1)) {
idamnjanovic@10 72 mexErrMsgTxt("COLLINCOMB requires that X be a double vector.");
idamnjanovic@10 73 }
idamnjanovic@10 74
idamnjanovic@10 75 if (m2!=colnum && n2!=colnum) {
idamnjanovic@10 76 mexErrMsgTxt("The length of X does not match the number of columns in COLS.");
idamnjanovic@10 77 }
idamnjanovic@10 78
idamnjanovic@10 79 rows = 0;
idamnjanovic@10 80 Y_OUT = mxCreateDoubleMatrix(m, 1, mxREAL);
idamnjanovic@10 81 cols = mxGetPr(COLS_IN1);
idamnjanovic@10 82 x = mxGetPr(X_IN1);
idamnjanovic@10 83 }
idamnjanovic@10 84 else {
idamnjanovic@10 85
idamnjanovic@10 86 m1 = mxGetM(ROWS_IN);
idamnjanovic@10 87 n1 = mxGetN(ROWS_IN);
idamnjanovic@10 88 if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) {
idamnjanovic@10 89 mexErrMsgTxt("COLLINCOMB requires that ROWS be an index vector of type double.");
idamnjanovic@10 90 }
idamnjanovic@10 91 rownum = (m1 > n1) ? m1 : n1; /* the number of rows in the linear combination */
idamnjanovic@10 92 rownumspecified = 1;
idamnjanovic@10 93 rows = mxGetPr(ROWS_IN);
idamnjanovic@10 94
idamnjanovic@10 95 m1 = mxGetM(COLS_IN2);
idamnjanovic@10 96 n1 = mxGetN(COLS_IN2);
idamnjanovic@10 97 if (!mxIsDouble(COLS_IN2) || mxIsComplex(COLS_IN2) || (m1!=1 && n1!=1)) {
idamnjanovic@10 98 mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double.");
idamnjanovic@10 99 }
idamnjanovic@10 100 colnum = (m1 > n1) ? m1 : n1; /* the number of columns in the linear combination */
idamnjanovic@10 101
idamnjanovic@10 102 m2 = mxGetM(X_IN2);
idamnjanovic@10 103 n2 = mxGetN(X_IN2);
idamnjanovic@10 104 if (!mxIsDouble(X_IN2) || mxIsComplex(X_IN2) || (m2!=1 && n2!=1)) {
idamnjanovic@10 105 mexErrMsgTxt("COLLINCOMB requires that X be a double vector.");
idamnjanovic@10 106 }
idamnjanovic@10 107
idamnjanovic@10 108 if (m2!=colnum && n2!=colnum) {
idamnjanovic@10 109 mexErrMsgTxt("The length of X does not match the number of columns in COLS.");
idamnjanovic@10 110 }
idamnjanovic@10 111
idamnjanovic@10 112 Y_OUT = mxCreateDoubleMatrix(rownum, 1, mxREAL);
idamnjanovic@10 113 cols = mxGetPr(COLS_IN2);
idamnjanovic@10 114 x = mxGetPr(X_IN2);
idamnjanovic@10 115 }
idamnjanovic@10 116
idamnjanovic@10 117
idamnjanovic@10 118 /* Assign pointers to the various parameters */
idamnjanovic@10 119 A = mxGetPr(A_IN);
idamnjanovic@10 120 y = mxGetPr(Y_OUT);
idamnjanovic@10 121
idamnjanovic@10 122
idamnjanovic@10 123 if (rownumspecified) {
idamnjanovic@10 124
idamnjanovic@10 125 /* check row indices */
idamnjanovic@10 126
idamnjanovic@10 127 row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex));
idamnjanovic@10 128
idamnjanovic@10 129 for (i=0; i<rownum; ++i) {
idamnjanovic@10 130 row_ids[i] = (mwIndex)(rows[i]+0.1)-1;
idamnjanovic@10 131 if (row_ids[i]<0 || row_ids[i]>=m) {
idamnjanovic@10 132 mexErrMsgTxt("Row index in ROWS is out of range.");
idamnjanovic@10 133 }
idamnjanovic@10 134 }
idamnjanovic@10 135
idamnjanovic@10 136 /* Do the actual computation */
idamnjanovic@10 137 for (i=0; i<colnum; ++i) {
idamnjanovic@10 138 col_id = (mwIndex)(cols[i]+0.1)-1;
idamnjanovic@10 139 if (col_id<0 || col_id>=n) {
idamnjanovic@10 140 mexErrMsgTxt("Column index in COLS is out of range.");
idamnjanovic@10 141 }
idamnjanovic@10 142 for (j=0; j<rownum; ++j) {
idamnjanovic@10 143 y[j] += A[m*col_id+row_ids[j]]*x[i];
idamnjanovic@10 144 }
idamnjanovic@10 145 }
idamnjanovic@10 146
idamnjanovic@10 147 mxFree(row_ids);
idamnjanovic@10 148 }
idamnjanovic@10 149
idamnjanovic@10 150 else {
idamnjanovic@10 151
idamnjanovic@10 152 /* Do the actual computation */
idamnjanovic@10 153 for (i=0; i<colnum; ++i) {
idamnjanovic@10 154 col_id = (mwIndex)(cols[i]+0.1)-1;
idamnjanovic@10 155 if (col_id<0 || col_id>=n) {
idamnjanovic@10 156 mexErrMsgTxt("Column index in COLS is out of range.");
idamnjanovic@10 157 }
idamnjanovic@10 158 for (j=0; j<m; ++j) {
idamnjanovic@10 159 y[j] += A[m*col_id+j]*x[i];
idamnjanovic@10 160 }
idamnjanovic@10 161 }
idamnjanovic@10 162 }
idamnjanovic@10 163
idamnjanovic@10 164 return;
idamnjanovic@10 165 }