source: trunk/src/grid/operate_inc.c @ 11591

Revision 11591, 10.7 KB checked in by joseba, 4 months ago (diff)

Added SVN Id to (almost) all the missing files, with svn propset svn:keywords "Id" FILE

  • Property svn:keywords set to Id
Line 
1/*
2 Copyright (C) 2006 X. Andrade
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
17 02110-1301, USA.
18
19 $Id$
20*/
21
22#ifdef ALIGNED
23#define LOAD VEC_LD
24#define STORE VEC_ST
25#else
26#define LOAD VEC_LDU
27#define STORE VEC_STU
28#endif
29
30{
31  const int n = opn[0];
32  const int nri = opnri[0];
33  const int unroll32 = max1(32*VEC_SIZE >> ldf);
34  const int unroll16 = max1(16*VEC_SIZE >> ldf);
35  const int unroll8  = max1(8*VEC_SIZE >> ldf);
36  const int unroll4  = max1(4*VEC_SIZE >> ldf);
37  const int unroll2  = max1(2*VEC_SIZE >> ldf);
38  const int unroll1  = max1(1*VEC_SIZE >> ldf);
39
40  int l, i, j;
41  const int * restrict index;
42   
43  for (l = 0; l < nri ; l++) {
44    index = opri + n * l;
45    i = rimap_inv[l];
46
47#if DEPTH >= 32
48    for (; i < (rimap_inv_max[l] - unroll32 + 1) ; i += unroll32){
49      int k;
50      for(k = 0; k < (1<<ldf); k += 32*VEC_SIZE){
51        register VEC_TYPE a0, a1, a2, a3;
52        register VEC_TYPE a4, a5, a6, a7;
53        register VEC_TYPE a8, a9, aa, ab;
54        register VEC_TYPE ac, ad, ae, af;
55        register VEC_TYPE b0, b1, b2, b3;
56        register VEC_TYPE b4, b5, b6, b7;
57        register VEC_TYPE b8, b9, ba, bb;
58        register VEC_TYPE bc, bd, be, bf;
59
60        a0 = a1 = a2 = a3 = VEC_ZERO;
61        a4 = a5 = a6 = a7 = VEC_ZERO;
62        a8 = a9 = aa = ab = VEC_ZERO;
63        ac = ad = ae = af = VEC_ZERO;
64        b0 = b1 = b2 = b3 = VEC_ZERO;
65        b4 = b5 = b6 = b7 = VEC_ZERO;
66        b8 = b9 = ba = bb = VEC_ZERO;
67        bc = bd = be = bf = VEC_ZERO;
68
69        for(j = 0; j < n; j++) {
70          register VEC_TYPE wj = VEC_SCAL(w[j]);
71          int indexj = (index[j] + i)<<ldf;
72          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
73          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
74          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
75          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
76          a4 = VEC_FMA(wj, LOAD(fi + indexj + 4*VEC_SIZE + k), a4);
77          a5 = VEC_FMA(wj, LOAD(fi + indexj + 5*VEC_SIZE + k), a5);
78          a6 = VEC_FMA(wj, LOAD(fi + indexj + 6*VEC_SIZE + k), a6);
79          a7 = VEC_FMA(wj, LOAD(fi + indexj + 7*VEC_SIZE + k), a7);
80          a8 = VEC_FMA(wj, LOAD(fi + indexj + 8*VEC_SIZE + k), a8);
81          a9 = VEC_FMA(wj, LOAD(fi + indexj + 9*VEC_SIZE + k), a9);
82          aa = VEC_FMA(wj, LOAD(fi + indexj + 10*VEC_SIZE + k), aa);
83          ab = VEC_FMA(wj, LOAD(fi + indexj + 11*VEC_SIZE + k), ab);
84          ac = VEC_FMA(wj, LOAD(fi + indexj + 12*VEC_SIZE + k), ac);
85          ad = VEC_FMA(wj, LOAD(fi + indexj + 13*VEC_SIZE + k), ad);
86          ae = VEC_FMA(wj, LOAD(fi + indexj + 14*VEC_SIZE + k), ae);
87          af = VEC_FMA(wj, LOAD(fi + indexj + 15*VEC_SIZE + k), af);
88          b0 = VEC_FMA(wj, LOAD(fi + indexj + 16*VEC_SIZE + k), b0);
89          b1 = VEC_FMA(wj, LOAD(fi + indexj + 17*VEC_SIZE + k), b1);
90          b2 = VEC_FMA(wj, LOAD(fi + indexj + 18*VEC_SIZE + k), b2);
91          b3 = VEC_FMA(wj, LOAD(fi + indexj + 19*VEC_SIZE + k), b3);
92          b4 = VEC_FMA(wj, LOAD(fi + indexj + 20*VEC_SIZE + k), b4);
93          b5 = VEC_FMA(wj, LOAD(fi + indexj + 21*VEC_SIZE + k), b5);
94          b6 = VEC_FMA(wj, LOAD(fi + indexj + 22*VEC_SIZE + k), b6);
95          b7 = VEC_FMA(wj, LOAD(fi + indexj + 23*VEC_SIZE + k), b7);
96          b8 = VEC_FMA(wj, LOAD(fi + indexj + 24*VEC_SIZE + k), b8);
97          b9 = VEC_FMA(wj, LOAD(fi + indexj + 25*VEC_SIZE + k), b9);
98          ba = VEC_FMA(wj, LOAD(fi + indexj + 26*VEC_SIZE + k), ba);
99          bb = VEC_FMA(wj, LOAD(fi + indexj + 27*VEC_SIZE + k), bb);
100          bc = VEC_FMA(wj, LOAD(fi + indexj + 28*VEC_SIZE + k), bc);
101          bd = VEC_FMA(wj, LOAD(fi + indexj + 29*VEC_SIZE + k), bd);
102          be = VEC_FMA(wj, LOAD(fi + indexj + 30*VEC_SIZE + k), be);
103          bf = VEC_FMA(wj, LOAD(fi + indexj + 31*VEC_SIZE + k), bf);
104
105        }
106        STORE(fo + (i<<ldf)              + k, a0);
107        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
108        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
109        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
110        STORE(fo + (i<<ldf) + 4*VEC_SIZE + k, a4);
111        STORE(fo + (i<<ldf) + 5*VEC_SIZE + k, a5);
112        STORE(fo + (i<<ldf) + 6*VEC_SIZE + k, a6);
113        STORE(fo + (i<<ldf) + 7*VEC_SIZE + k, a7);
114        STORE(fo + (i<<ldf) + 8*VEC_SIZE + k, a8);
115        STORE(fo + (i<<ldf) + 9*VEC_SIZE + k, a9);
116        STORE(fo + (i<<ldf) + 10*VEC_SIZE + k, aa);
117        STORE(fo + (i<<ldf) + 11*VEC_SIZE + k, ab);
118        STORE(fo + (i<<ldf) + 12*VEC_SIZE + k, ac);
119        STORE(fo + (i<<ldf) + 13*VEC_SIZE + k, ad);
120        STORE(fo + (i<<ldf) + 14*VEC_SIZE + k, ae);
121        STORE(fo + (i<<ldf) + 15*VEC_SIZE + k, af);
122        STORE(fo + (i<<ldf) + 16*VEC_SIZE + k, b0);
123        STORE(fo + (i<<ldf) + 17*VEC_SIZE + k, b1);
124        STORE(fo + (i<<ldf) + 18*VEC_SIZE + k, b2);
125        STORE(fo + (i<<ldf) + 19*VEC_SIZE + k, b3);
126        STORE(fo + (i<<ldf) + 20*VEC_SIZE + k, b4);
127        STORE(fo + (i<<ldf) + 21*VEC_SIZE + k, b5);
128        STORE(fo + (i<<ldf) + 22*VEC_SIZE + k, b6);
129        STORE(fo + (i<<ldf) + 23*VEC_SIZE + k, b7);
130        STORE(fo + (i<<ldf) + 24*VEC_SIZE + k, b8);
131        STORE(fo + (i<<ldf) + 25*VEC_SIZE + k, b9);
132        STORE(fo + (i<<ldf) + 26*VEC_SIZE + k, ba);
133        STORE(fo + (i<<ldf) + 27*VEC_SIZE + k, bb);
134        STORE(fo + (i<<ldf) + 28*VEC_SIZE + k, bc);
135        STORE(fo + (i<<ldf) + 29*VEC_SIZE + k, bd);
136        STORE(fo + (i<<ldf) + 30*VEC_SIZE + k, be);
137        STORE(fo + (i<<ldf) + 31*VEC_SIZE + k, bf);
138      }
139    }
140#endif
141
142#if DEPTH >= 16
143    for (; i < (rimap_inv_max[l] - unroll16 + 1) ; i += unroll16){
144      int k;
145      for(k = 0; k < (1<<ldf); k += 16*VEC_SIZE){
146        register VEC_TYPE a0, a1, a2, a3;
147        register VEC_TYPE a4, a5, a6, a7;
148        register VEC_TYPE a8, a9, aa, ab;
149        register VEC_TYPE ac, ad, ae, af;
150
151        a0 = a1 = a2 = a3 = VEC_ZERO;
152        a4 = a5 = a6 = a7 = VEC_ZERO;
153        a8 = a9 = aa = ab = VEC_ZERO;
154        ac = ad = ae = af = VEC_ZERO;
155
156        for(j = 0; j < n; j++) {
157          register VEC_TYPE wj = VEC_SCAL(w[j]);
158          int indexj = (index[j] + i)<<ldf;
159          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
160          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
161          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
162          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
163          a4 = VEC_FMA(wj, LOAD(fi + indexj + 4*VEC_SIZE + k), a4);
164          a5 = VEC_FMA(wj, LOAD(fi + indexj + 5*VEC_SIZE + k), a5);
165          a6 = VEC_FMA(wj, LOAD(fi + indexj + 6*VEC_SIZE + k), a6);
166          a7 = VEC_FMA(wj, LOAD(fi + indexj + 7*VEC_SIZE + k), a7);
167          a8 = VEC_FMA(wj, LOAD(fi + indexj + 8*VEC_SIZE + k), a8);
168          a9 = VEC_FMA(wj, LOAD(fi + indexj + 9*VEC_SIZE + k), a9);
169          aa = VEC_FMA(wj, LOAD(fi + indexj + 10*VEC_SIZE + k), aa);
170          ab = VEC_FMA(wj, LOAD(fi + indexj + 11*VEC_SIZE + k), ab);
171          ac = VEC_FMA(wj, LOAD(fi + indexj + 12*VEC_SIZE + k), ac);
172          ad = VEC_FMA(wj, LOAD(fi + indexj + 13*VEC_SIZE + k), ad);
173          ae = VEC_FMA(wj, LOAD(fi + indexj + 14*VEC_SIZE + k), ae);
174          af = VEC_FMA(wj, LOAD(fi + indexj + 15*VEC_SIZE + k), af);
175        }
176        STORE(fo + (i<<ldf)              + k, a0);
177        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
178        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
179        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
180        STORE(fo + (i<<ldf) + 4*VEC_SIZE + k, a4);
181        STORE(fo + (i<<ldf) + 5*VEC_SIZE + k, a5);
182        STORE(fo + (i<<ldf) + 6*VEC_SIZE + k, a6);
183        STORE(fo + (i<<ldf) + 7*VEC_SIZE + k, a7);
184        STORE(fo + (i<<ldf) + 8*VEC_SIZE + k, a8);
185        STORE(fo + (i<<ldf) + 9*VEC_SIZE + k, a9);
186        STORE(fo + (i<<ldf) + 10*VEC_SIZE + k, aa);
187        STORE(fo + (i<<ldf) + 11*VEC_SIZE + k, ab);
188        STORE(fo + (i<<ldf) + 12*VEC_SIZE + k, ac);
189        STORE(fo + (i<<ldf) + 13*VEC_SIZE + k, ad);
190        STORE(fo + (i<<ldf) + 14*VEC_SIZE + k, ae);
191        STORE(fo + (i<<ldf) + 15*VEC_SIZE + k, af);
192      }
193    }
194#endif
195
196#if DEPTH >= 8
197    for (; i < (rimap_inv_max[l] - unroll8 + 1) ; i += unroll8){
198      int k;
199      for(k = 0; k < (1<<ldf); k += 8*VEC_SIZE){
200        register VEC_TYPE a0, a1, a2, a3;
201        register VEC_TYPE a4, a5, a6, a7;
202
203        a0 = a1 = a2 = a3 = VEC_ZERO;
204        a4 = a5 = a6 = a7 = VEC_ZERO;
205
206        for(j = 0; j < n; j++) {
207          register VEC_TYPE wj = VEC_SCAL(w[j]);
208          int indexj = (index[j] + i)<<ldf;
209          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
210          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
211          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
212          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
213          a4 = VEC_FMA(wj, LOAD(fi + indexj + 4*VEC_SIZE + k), a4);
214          a5 = VEC_FMA(wj, LOAD(fi + indexj + 5*VEC_SIZE + k), a5);
215          a6 = VEC_FMA(wj, LOAD(fi + indexj + 6*VEC_SIZE + k), a6);
216          a7 = VEC_FMA(wj, LOAD(fi + indexj + 7*VEC_SIZE + k), a7);
217        }
218        STORE(fo + (i<<ldf)              + k, a0);
219        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
220        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
221        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
222        STORE(fo + (i<<ldf) + 4*VEC_SIZE + k, a4);
223        STORE(fo + (i<<ldf) + 5*VEC_SIZE + k, a5);
224        STORE(fo + (i<<ldf) + 6*VEC_SIZE + k, a6);
225        STORE(fo + (i<<ldf) + 7*VEC_SIZE + k, a7);
226      }
227    }
228#endif
229
230#if DEPTH >= 4
231    for (; i < (rimap_inv_max[l] - unroll4 + 1) ; i += unroll4){
232      int k;
233      for(k = 0; k < (1<<ldf); k += 4*VEC_SIZE){
234        register VEC_TYPE a0, a1, a2, a3;
235
236        a0 = a1 = a2 = a3 = VEC_ZERO;
237
238        for(j = 0; j < n; j++) {
239          register VEC_TYPE wj = VEC_SCAL(w[j]);
240          int indexj = (index[j] + i)<<ldf;
241          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
242          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
243          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
244          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
245        }
246        STORE(fo + (i<<ldf)              + k, a0);
247        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
248        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
249        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
250      }
251    }
252#endif
253
254#if DEPTH >= 2
255    for (; i < (rimap_inv_max[l] - unroll2 + 1) ; i += unroll2){
256      int k;
257      for(k = 0; k < (1<<ldf); k += 2*VEC_SIZE){
258        register VEC_TYPE a0, a1;
259
260        a0 = a1 = VEC_ZERO;
261
262        for(j = 0; j < n; j++) {
263          register VEC_TYPE wj = VEC_SCAL(w[j]);
264          int indexj = (index[j] + i)<<ldf;
265          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
266          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
267        }
268        STORE(fo + (i<<ldf)              + k, a0);
269        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
270      }
271    }
272#endif
273
274#if DEPTH >= 1
275    for (; i < (rimap_inv_max[l] - unroll1 + 1) ; i += unroll1){
276      int k;
277      for(k = 0; k < (1<<ldf); k += VEC_SIZE){
278        register VEC_TYPE a0 = VEC_ZERO;
279        for(j = 0; j < n; j++) {
280          register VEC_TYPE wj = VEC_SCAL(w[j]);
281          int indexj = (index[j] + i)<<ldf;
282          a0 = VEC_FMA(wj, LOAD(fi + indexj + k), a0);
283        }
284        STORE(fo + (i<<ldf) + k, a0);
285      }
286    }
287#endif
288
289#if VEC_SIZE > 1
290    for (; i < rimap_inv_max[l]; i++){
291      int k;
292      for(k = 0; k < (1<<ldf); k++){
293        double a = 0.0;
294        for(j = 0; j < n; j++) a += w[j]*fi[((index[j] + i)<<ldf) + k];
295        fo[(i<<ldf) + k] = a;
296      }
297    }
298#endif
299
300  } /* l */
301
302}
303
304#undef LOAD
305#undef STORE
Note: See TracBrowser for help on using the repository browser.