comparison DL/RLS-DLA/private/collincomb.c @ 60:ad36f80e2ccf

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