1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2012 - INRIA - Serge STEER
4 *
5 * Copyright (C) 2012 - 2016 - Scilab Enterprises
6 *
7 * This file is hereby licensed under the terms of the GNU GPL v2.0,
8 * pursuant to article 5.3.4 of the CeCILL v.2.1.
9 * This file was originally licensed under the terms of the CeCILL v2.1,
10 * and continues to be available under such terms.
11 * For more information, see the COPYING file which you should have received
12 * along with this program.
13 *
14 */
15 /*--------------------------------------------------------------------------*/
16 #include "machine.h"
17 #include "core_math.h"
18 #include "sci_malloc.h"
19 #include "conv2.h"
20 /*--------------------------------------------------------------------------*/
21 extern double C2F(ddot)(int *n, double *A, int *iA, double *B, int *iB);
22 /*--------------------------------------------------------------------------*/
conv2_separable_R(double * R,int nR,double * C,int mC,double * A,int mA,int nA,double * Out,int mOut,int nOut,int edgM,int edgN,double * T)23 void conv2_separable_R(double *R, int nR, double *C, int mC, double *A, int mA, int nA, double *Out, int mOut, int nOut, int edgM, int edgN, double *T)
24 {
25 int ai = 0, tj = 0, ci = 0, rj = 0; /*current index over A,T,C and R */
26 int i = 0, j = 0; /* loop variables*/
27 int l = 0;
28 int one = 1, minusone = -1;
29
30 for (i = 0; i < mOut; i++ )
31 {
32 /*Compute the 1-D conv A(i,:) and C in T */
33 ai = Max(0, i - edgM) ;
34 ci = mC - 1 - Max(0, edgM - i);
35 l = Min(ci + 1, mA - ai);
36 for (j = 0; j < nA; j++ )
37 {
38 T[j] = C2F(ddot)(&l, A + ai + mA * j, &one, C + ci - l + 1, &minusone);
39 }
40 /*1-D convolution of T and R */
41 for (j = 0; j < nOut; j++ )
42 {
43 rj = nR - 1 - Max(0, edgN - j);
44 tj = Max(0, j - edgN) ;
45 l = Min(rj + 1, nA - tj);
46 Out[i + j * mOut] = C2F(ddot)(&l, T + tj, &one, R + rj - l + 1, &minusone);
47 }
48 }
49 }
50 /*--------------------------------------------------------------------------*/
conv2_separable_C(double * Rr,double * Ri,int nR,double * Cr,double * Ci,int mC,double * Ar,double * Ai,int mA,int nA,double * Outr,double * Outi,int mOut,int nOut,int edgM,int edgN,double * Tr,double * Ti)51 void conv2_separable_C(double *Rr, double *Ri, int nR, double *Cr, double *Ci, int mC, double *Ar, double *Ai, int mA, int nA, double *Outr, double *Outi, int mOut, int nOut, int edgM, int edgN, double *Tr, double *Ti)
52 {
53 int ai = 0, tj = 0, ci = 0, rj = 0; /*current index over A,T,C and R */
54 int i = 0, j = 0; /* loop variables*/
55 int l = 0;
56 int one = 1, minusone = -1;
57
58 for (i = 0; i < mOut; i++ )
59 {
60 /*Compute the 1-D conv A(i,:) and C in T */
61 ai = Max(0, i - edgM) ;
62 ci = mC - 1 - Max(0, edgM - i);
63 l = Min(ci + 1, mA - ai);
64 if (Ai != NULL && Ci != NULL)
65 {
66 for (j = 0; j < nA; j++ )
67 {
68 Tr[j] = C2F(ddot)(&l, Ar + ai + mA * j, &one, Cr + ci - l + 1, &minusone) -
69 C2F(ddot)(&l, Ai + ai + mA * j, &one, Ci + ci - l + 1, &minusone);
70
71 Ti[j] = C2F(ddot)(&l, Ar + ai + mA * j, &one, Ci + ci - l + 1, &minusone) +
72 C2F(ddot)(&l, Ai + ai + mA * j, &one, Cr + ci - l + 1, &minusone);
73 }
74 }
75 else if (Ci != NULL)
76 {
77 for (j = 0; j < nA; j++ )
78 {
79 Tr[j] = C2F(ddot)(&l, Ar + ai + mA * j, &one, Cr + ci - l + 1, &minusone);
80 Ti[j] = C2F(ddot)(&l, Ar + ai + mA * j, &one, Ci + ci - l + 1, &minusone);
81 }
82 }
83 else if (Ai != NULL)
84 {
85 for (j = 0; j < nA; j++ )
86 {
87 Tr[j] = C2F(ddot)(&l, Ar + ai + mA * j, &one, Cr + ci - l + 1, &minusone);
88 Ti[j] = C2F(ddot)(&l, Ai + ai + mA * j, &one, Cr + ci - l + 1, &minusone);
89 }
90 }
91 else
92 {
93 for (j = 0; j < nA; j++ )
94 {
95 Tr[j] = C2F(ddot)(&l, Ar + ai + mA * j, &one, Cr + ci - l + 1, &minusone);
96 Ti[j] = 0.0;
97 }
98 }
99 /*1-D convolution of T and R */
100 for (j = 0; j < nOut; j++ )
101 {
102 rj = nR - 1 - Max(0, edgN - j);
103 tj = Max(0, j - edgN) ;
104 l = Min(rj + 1, nA - tj);
105 Outr[i + j * mOut] = C2F(ddot)(&l, Tr + tj, &one, Rr + rj - l + 1, &minusone);
106 Outi[i + j * mOut] = C2F(ddot)(&l, Ti + tj, &one, Rr + rj - l + 1, &minusone);
107 if (Ri != NULL)
108 {
109 Outr[i + j * mOut] -= C2F(ddot)(&l, Ti + tj, &one, Ri + rj - l + 1, &minusone);
110 Outi[i + j * mOut] = C2F(ddot)(&l, Tr + tj, &one, Ri + rj - l + 1, &minusone);
111 }
112 }
113 }
114 }
115 /*--------------------------------------------------------------------------*/
conv2_R(double * A,int mA,int nA,double * B,int mB,int nB,double * Out,int mOut,int nOut,int edgM,int edgN)116 void conv2_R(double *A, int mA, int nA, double *B, int mB, int nB, double *Out, int mOut, int nOut, int edgM, int edgN)
117 {
118 int ai = 0, aj = 0, bi = 0, bj = 0; /*current index over A and B */
119 int i = 0, j = 0; /* loop variables*/
120 int l = 0;
121 int one = 1, minusone = -1;
122 double sum = 0;
123 if (nOut == 1)
124 {
125 /* A and B are column vectors nA=nB=nOut=1 */
126 for (i = 0; i < mOut; i++ )
127 {
128 bi = mB - 1 - Max(0, edgM - i);
129 ai = Max(0, i - edgM);
130 l = Min(bi + 1, mA - ai);
131 Out[i] = C2F(ddot)(&l, A + ai, &one, B + bi - l + 1, &minusone);
132 }
133 }
134 else if (mOut == 1)
135 {
136 /* A and B are row vectors mA=mB=mOut=1 */
137 for (j = 0; j < nOut; j++ )
138 {
139 bj = nB - 1 - Max(0, edgN - j);
140 aj = Max(0, j - edgN);
141 l = Min(bj + 1, nA - aj);
142 Out[j] = C2F(ddot)(&l, A + aj, &one, B + bj - l + 1, &minusone);
143 }
144 }
145 else
146 {
147 /* general array case */
148 for (i = 0; i < mOut; i++ )
149 {
150 bi = mB - 1 - Max(0, edgM - i);
151 ai = Max(0, i - edgM);
152 for (j = 0; j < nOut; j++ )
153 {
154 sum = 0;
155 for (bj = nB - 1 - Max(0, edgN - j), aj = Max(0, j - edgN); bj >= 0 && aj < nA; bj--, aj++)
156 {
157 l = Min(bi + 1, mA - ai);
158 sum += C2F(ddot)(&l, A + ai + mA * aj, &one, B + bi - l + 1 + mB * bj, &minusone);
159 }
160 Out[i + j * mOut] = sum;
161 }
162 }
163 }
164 }
165 /*--------------------------------------------------------------------------*/
conv2_C(double * Ar,double * Ai,int mA,int nA,double * Br,double * Bi,int mB,int nB,double * Outr,double * Outi,int mOut,int nOut,int edgM,int edgN)166 void conv2_C(double *Ar, double *Ai, int mA, int nA, double *Br, double *Bi, int mB, int nB, double *Outr, double *Outi, int mOut, int nOut, int edgM, int edgN)
167 {
168 int ai = 0, aj = 0, bi = 0, bj = 0; /*current index over A and B */
169 int i = 0, j = 0; /* loop variables*/
170 int l = 0;
171 int one = 1, minusone = -1;
172 double sumr = 0, sumi = 0;
173
174 if (Bi != NULL && Ai != NULL)
175 {
176 if (nOut == 1)
177 {
178 /* fastest code */
179 for (i = 0; i < mOut; i++ )
180 {
181 bi = mB - 1 - Max(0, edgM - i);
182 ai = Max(0, i - edgM);
183 l = Min(bi + 1, mA - ai);
184 Outr[i] = C2F(ddot)(&l, Ar + ai, &one, Br + bi - l + 1, &minusone) -
185 C2F(ddot)(&l, Ai + ai, &one, Bi + bi - l + 1, &minusone);
186
187 Outi[i] = C2F(ddot)(&l, Ar + ai, &one, Bi + bi - l + 1, &minusone) +
188 C2F(ddot)(&l, Ai + ai, &one, Br + bi - l + 1, &minusone);
189 }
190 }
191 else
192 {
193 for (i = 0; i < mOut; i++ )
194 {
195 bi = mB - 1 - Max(0, edgM - i);
196 ai = Max(0, i - edgM);
197 for (j = 0; j < nOut; j++ )
198 {
199 sumr = 0;
200 sumi = 0;
201 for (bj = nB - 1 - Max(0, edgN - j), aj = Max(0, j - edgN); bj >= 0 && aj < nA; bj--, aj++)
202 {
203 l = Min(bi + 1, mA - ai);
204 sumr += C2F(ddot)(&l, Ar + ai + mA * aj, &one, Br + bi - l + 1 + mB * bj, &minusone)
205 - C2F(ddot)(&l, Ai + ai + mA * aj, &one, Bi + bi - l + 1 + mB * bj, &minusone);
206
207 sumi += C2F(ddot)(&l, Ar + ai + mA * aj, &one, Bi + bi - l + 1 + mB * bj, &minusone) +
208 C2F(ddot)(&l, Ai + ai + mA * aj, &one, Br + bi - l + 1 + mB * bj, &minusone);
209 }
210 Outr[i + j * mOut] = sumr;
211 Outi[i + j * mOut] = sumi;
212 }
213 }
214 }
215 }
216 else if (Ai != NULL)
217 {
218 if (nOut == 1)
219 {
220 /* fastest code */
221 for (i = 0; i < mOut; i++ )
222 {
223 bi = mB - 1 - Max(0, edgM - i);
224 ai = Max(0, i - edgM);
225 l = Min(bi + 1, mA - ai);
226 Outr[i] = C2F(ddot)(&l, Ar + ai, &one, Br + bi - l + 1, &minusone);
227 Outi[i] = C2F(ddot)(&l, Ai + ai, &one, Br + bi - l + 1, &minusone);
228 }
229 }
230 else
231 {
232 for (i = 0; i < mOut; i++ )
233 {
234 bi = mB - 1 - Max(0, edgM - i);
235 ai = Max(0, i - edgM);
236 for (j = 0; j < nOut; j++ )
237 {
238 sumr = 0;
239 sumi = 0;
240 for (bj = nB - 1 - Max(0, edgN - j), aj = Max(0, j - edgN); bj >= 0 && aj < nA; bj--, aj++)
241 {
242 l = Min(bi + 1, mA - ai);
243 sumr += C2F(ddot)(&l, Ar + ai + mA * aj, &one, Br + bi - l + 1 + mB * bj, &minusone);
244 sumi += C2F(ddot)(&l, Ai + ai + mA * aj, &one, Br + bi - l + 1 + mB * bj, &minusone);
245 }
246 Outr[i + j * mOut] = sumr;
247 Outi[i + j * mOut] = sumi;
248 }
249 }
250 }
251 }
252 else if (Bi != NULL)
253 {
254 if (nOut == 1)
255 {
256 /* fastest code */
257 for (i = 0; i < mOut; i++ )
258 {
259 bi = mB - 1 - Max(0, edgM - i);
260 ai = Max(0, i - edgM);
261 l = Min(bi + 1, mA - ai);
262 Outr[i] = C2F(ddot)(&l, Ar + ai, &one, Br + bi - l + 1, &minusone);
263 Outi[i] = C2F(ddot)(&l, Ar + ai, &one, Bi + bi - l + 1, &minusone);
264 }
265 }
266 else
267 {
268 for (i = 0; i < mOut; i++ )
269 {
270 bi = mB - 1 - Max(0, edgM - i);
271 ai = Max(0, i - edgM);
272 for (j = 0; j < nOut; j++ )
273 {
274 sumr = 0;
275 sumi = 0;
276 for (bj = nB - 1 - Max(0, edgN - j), aj = Max(0, j - edgN); bj >= 0 && aj < nA; bj--, aj++)
277 {
278 l = Min(bi + 1, mA - ai);
279 sumr += C2F(ddot)(&l, Ar + ai + mA * aj, &one, Br + bi - l + 1 + mB * bj, &minusone);
280 sumi += C2F(ddot)(&l, Ar + ai + mA * aj, &one, Bi + bi - l + 1 + mB * bj, &minusone);
281 }
282 Outr[i + j * mOut] = sumr;
283 Outi[i + j * mOut] = sumi;
284 }
285 }
286 }
287 }
288 }
289 /*--------------------------------------------------------------------------*/
290