comparison Problems/private/rowlincomb.c @ 10:207a6ae9a76f version1.0

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