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

Revision 12581, 11.1 KB checked in by xavier, 5 weeks ago (diff)

Small change: Improvement of the vectorization for the Blue Gene/Q.

  • 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#ifdef VEC_SCAL_LD
71          register VEC_TYPE wj = VEC_SCAL_LD(w + j);
72#else     
73          register VEC_TYPE wj = VEC_SCAL(w[j]);
74#endif
75          int indexj = (index[j] + i)<<ldf;
76          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
77          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
78          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
79          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
80          a4 = VEC_FMA(wj, LOAD(fi + indexj + 4*VEC_SIZE + k), a4);
81          a5 = VEC_FMA(wj, LOAD(fi + indexj + 5*VEC_SIZE + k), a5);
82          a6 = VEC_FMA(wj, LOAD(fi + indexj + 6*VEC_SIZE + k), a6);
83          a7 = VEC_FMA(wj, LOAD(fi + indexj + 7*VEC_SIZE + k), a7);
84          a8 = VEC_FMA(wj, LOAD(fi + indexj + 8*VEC_SIZE + k), a8);
85          a9 = VEC_FMA(wj, LOAD(fi + indexj + 9*VEC_SIZE + k), a9);
86          aa = VEC_FMA(wj, LOAD(fi + indexj + 10*VEC_SIZE + k), aa);
87          ab = VEC_FMA(wj, LOAD(fi + indexj + 11*VEC_SIZE + k), ab);
88          ac = VEC_FMA(wj, LOAD(fi + indexj + 12*VEC_SIZE + k), ac);
89          ad = VEC_FMA(wj, LOAD(fi + indexj + 13*VEC_SIZE + k), ad);
90          ae = VEC_FMA(wj, LOAD(fi + indexj + 14*VEC_SIZE + k), ae);
91          af = VEC_FMA(wj, LOAD(fi + indexj + 15*VEC_SIZE + k), af);
92          b0 = VEC_FMA(wj, LOAD(fi + indexj + 16*VEC_SIZE + k), b0);
93          b1 = VEC_FMA(wj, LOAD(fi + indexj + 17*VEC_SIZE + k), b1);
94          b2 = VEC_FMA(wj, LOAD(fi + indexj + 18*VEC_SIZE + k), b2);
95          b3 = VEC_FMA(wj, LOAD(fi + indexj + 19*VEC_SIZE + k), b3);
96          b4 = VEC_FMA(wj, LOAD(fi + indexj + 20*VEC_SIZE + k), b4);
97          b5 = VEC_FMA(wj, LOAD(fi + indexj + 21*VEC_SIZE + k), b5);
98          b6 = VEC_FMA(wj, LOAD(fi + indexj + 22*VEC_SIZE + k), b6);
99          b7 = VEC_FMA(wj, LOAD(fi + indexj + 23*VEC_SIZE + k), b7);
100          b8 = VEC_FMA(wj, LOAD(fi + indexj + 24*VEC_SIZE + k), b8);
101          b9 = VEC_FMA(wj, LOAD(fi + indexj + 25*VEC_SIZE + k), b9);
102          ba = VEC_FMA(wj, LOAD(fi + indexj + 26*VEC_SIZE + k), ba);
103          bb = VEC_FMA(wj, LOAD(fi + indexj + 27*VEC_SIZE + k), bb);
104          bc = VEC_FMA(wj, LOAD(fi + indexj + 28*VEC_SIZE + k), bc);
105          bd = VEC_FMA(wj, LOAD(fi + indexj + 29*VEC_SIZE + k), bd);
106          be = VEC_FMA(wj, LOAD(fi + indexj + 30*VEC_SIZE + k), be);
107          bf = VEC_FMA(wj, LOAD(fi + indexj + 31*VEC_SIZE + k), bf);
108
109        }
110        STORE(fo + (i<<ldf)              + k, a0);
111        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
112        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
113        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
114        STORE(fo + (i<<ldf) + 4*VEC_SIZE + k, a4);
115        STORE(fo + (i<<ldf) + 5*VEC_SIZE + k, a5);
116        STORE(fo + (i<<ldf) + 6*VEC_SIZE + k, a6);
117        STORE(fo + (i<<ldf) + 7*VEC_SIZE + k, a7);
118        STORE(fo + (i<<ldf) + 8*VEC_SIZE + k, a8);
119        STORE(fo + (i<<ldf) + 9*VEC_SIZE + k, a9);
120        STORE(fo + (i<<ldf) + 10*VEC_SIZE + k, aa);
121        STORE(fo + (i<<ldf) + 11*VEC_SIZE + k, ab);
122        STORE(fo + (i<<ldf) + 12*VEC_SIZE + k, ac);
123        STORE(fo + (i<<ldf) + 13*VEC_SIZE + k, ad);
124        STORE(fo + (i<<ldf) + 14*VEC_SIZE + k, ae);
125        STORE(fo + (i<<ldf) + 15*VEC_SIZE + k, af);
126        STORE(fo + (i<<ldf) + 16*VEC_SIZE + k, b0);
127        STORE(fo + (i<<ldf) + 17*VEC_SIZE + k, b1);
128        STORE(fo + (i<<ldf) + 18*VEC_SIZE + k, b2);
129        STORE(fo + (i<<ldf) + 19*VEC_SIZE + k, b3);
130        STORE(fo + (i<<ldf) + 20*VEC_SIZE + k, b4);
131        STORE(fo + (i<<ldf) + 21*VEC_SIZE + k, b5);
132        STORE(fo + (i<<ldf) + 22*VEC_SIZE + k, b6);
133        STORE(fo + (i<<ldf) + 23*VEC_SIZE + k, b7);
134        STORE(fo + (i<<ldf) + 24*VEC_SIZE + k, b8);
135        STORE(fo + (i<<ldf) + 25*VEC_SIZE + k, b9);
136        STORE(fo + (i<<ldf) + 26*VEC_SIZE + k, ba);
137        STORE(fo + (i<<ldf) + 27*VEC_SIZE + k, bb);
138        STORE(fo + (i<<ldf) + 28*VEC_SIZE + k, bc);
139        STORE(fo + (i<<ldf) + 29*VEC_SIZE + k, bd);
140        STORE(fo + (i<<ldf) + 30*VEC_SIZE + k, be);
141        STORE(fo + (i<<ldf) + 31*VEC_SIZE + k, bf);
142      }
143    }
144#endif
145
146#if DEPTH >= 16
147    for (; i < (rimap_inv_max[l] - unroll16 + 1) ; i += unroll16){
148      int k;
149      for(k = 0; k < (1<<ldf); k += 16*VEC_SIZE){
150        register VEC_TYPE a0, a1, a2, a3;
151        register VEC_TYPE a4, a5, a6, a7;
152        register VEC_TYPE a8, a9, aa, ab;
153        register VEC_TYPE ac, ad, ae, af;
154
155        a0 = a1 = a2 = a3 = VEC_ZERO;
156        a4 = a5 = a6 = a7 = VEC_ZERO;
157        a8 = a9 = aa = ab = VEC_ZERO;
158        ac = ad = ae = af = VEC_ZERO;
159
160        for(j = 0; j < n; j++) {
161#ifdef VEC_SCAL_LD
162          register VEC_TYPE wj = VEC_SCAL_LD(w + j);
163#else
164          register VEC_TYPE wj = VEC_SCAL(w[j]);
165#endif
166          int indexj = (index[j] + i)<<ldf;
167          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
168          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
169          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
170          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
171          a4 = VEC_FMA(wj, LOAD(fi + indexj + 4*VEC_SIZE + k), a4);
172          a5 = VEC_FMA(wj, LOAD(fi + indexj + 5*VEC_SIZE + k), a5);
173          a6 = VEC_FMA(wj, LOAD(fi + indexj + 6*VEC_SIZE + k), a6);
174          a7 = VEC_FMA(wj, LOAD(fi + indexj + 7*VEC_SIZE + k), a7);
175          a8 = VEC_FMA(wj, LOAD(fi + indexj + 8*VEC_SIZE + k), a8);
176          a9 = VEC_FMA(wj, LOAD(fi + indexj + 9*VEC_SIZE + k), a9);
177          aa = VEC_FMA(wj, LOAD(fi + indexj + 10*VEC_SIZE + k), aa);
178          ab = VEC_FMA(wj, LOAD(fi + indexj + 11*VEC_SIZE + k), ab);
179          ac = VEC_FMA(wj, LOAD(fi + indexj + 12*VEC_SIZE + k), ac);
180          ad = VEC_FMA(wj, LOAD(fi + indexj + 13*VEC_SIZE + k), ad);
181          ae = VEC_FMA(wj, LOAD(fi + indexj + 14*VEC_SIZE + k), ae);
182          af = VEC_FMA(wj, LOAD(fi + indexj + 15*VEC_SIZE + k), af);
183        }
184        STORE(fo + (i<<ldf)              + k, a0);
185        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
186        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
187        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
188        STORE(fo + (i<<ldf) + 4*VEC_SIZE + k, a4);
189        STORE(fo + (i<<ldf) + 5*VEC_SIZE + k, a5);
190        STORE(fo + (i<<ldf) + 6*VEC_SIZE + k, a6);
191        STORE(fo + (i<<ldf) + 7*VEC_SIZE + k, a7);
192        STORE(fo + (i<<ldf) + 8*VEC_SIZE + k, a8);
193        STORE(fo + (i<<ldf) + 9*VEC_SIZE + k, a9);
194        STORE(fo + (i<<ldf) + 10*VEC_SIZE + k, aa);
195        STORE(fo + (i<<ldf) + 11*VEC_SIZE + k, ab);
196        STORE(fo + (i<<ldf) + 12*VEC_SIZE + k, ac);
197        STORE(fo + (i<<ldf) + 13*VEC_SIZE + k, ad);
198        STORE(fo + (i<<ldf) + 14*VEC_SIZE + k, ae);
199        STORE(fo + (i<<ldf) + 15*VEC_SIZE + k, af);
200      }
201    }
202#endif
203
204#if DEPTH >= 8
205    for (; i < (rimap_inv_max[l] - unroll8 + 1) ; i += unroll8){
206      int k;
207      for(k = 0; k < (1<<ldf); k += 8*VEC_SIZE){
208        register VEC_TYPE a0, a1, a2, a3;
209        register VEC_TYPE a4, a5, a6, a7;
210
211        a0 = a1 = a2 = a3 = VEC_ZERO;
212        a4 = a5 = a6 = a7 = VEC_ZERO;
213
214        for(j = 0; j < n; j++) {
215#ifdef VEC_SCAL_LD
216          register VEC_TYPE wj = VEC_SCAL_LD(w + j);
217#else
218          register VEC_TYPE wj = VEC_SCAL(w[j]);
219#endif
220          int indexj = (index[j] + i)<<ldf;
221          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
222          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
223          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
224          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
225          a4 = VEC_FMA(wj, LOAD(fi + indexj + 4*VEC_SIZE + k), a4);
226          a5 = VEC_FMA(wj, LOAD(fi + indexj + 5*VEC_SIZE + k), a5);
227          a6 = VEC_FMA(wj, LOAD(fi + indexj + 6*VEC_SIZE + k), a6);
228          a7 = VEC_FMA(wj, LOAD(fi + indexj + 7*VEC_SIZE + k), a7);
229        }
230        STORE(fo + (i<<ldf)              + k, a0);
231        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
232        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
233        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
234        STORE(fo + (i<<ldf) + 4*VEC_SIZE + k, a4);
235        STORE(fo + (i<<ldf) + 5*VEC_SIZE + k, a5);
236        STORE(fo + (i<<ldf) + 6*VEC_SIZE + k, a6);
237        STORE(fo + (i<<ldf) + 7*VEC_SIZE + k, a7);
238      }
239    }
240#endif
241
242#if DEPTH >= 4
243    for (; i < (rimap_inv_max[l] - unroll4 + 1) ; i += unroll4){
244      int k;
245      for(k = 0; k < (1<<ldf); k += 4*VEC_SIZE){
246        register VEC_TYPE a0, a1, a2, a3;
247
248        a0 = a1 = a2 = a3 = VEC_ZERO;
249
250        for(j = 0; j < n; j++) {
251#ifdef VEC_SCAL_LD
252          register VEC_TYPE wj = VEC_SCAL_LD(w + j);
253#else
254          register VEC_TYPE wj = VEC_SCAL(w[j]);
255#endif
256          int indexj = (index[j] + i)<<ldf;
257          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
258          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
259          a2 = VEC_FMA(wj, LOAD(fi + indexj + 2*VEC_SIZE + k), a2);
260          a3 = VEC_FMA(wj, LOAD(fi + indexj + 3*VEC_SIZE + k), a3);
261        }
262        STORE(fo + (i<<ldf)              + k, a0);
263        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
264        STORE(fo + (i<<ldf) + 2*VEC_SIZE + k, a2);
265        STORE(fo + (i<<ldf) + 3*VEC_SIZE + k, a3);
266      }
267    }
268#endif
269
270#if DEPTH >= 2
271    for (; i < (rimap_inv_max[l] - unroll2 + 1) ; i += unroll2){
272      int k;
273      for(k = 0; k < (1<<ldf); k += 2*VEC_SIZE){
274        register VEC_TYPE a0, a1;
275
276        a0 = a1 = VEC_ZERO;
277
278        for(j = 0; j < n; j++) {
279          register VEC_TYPE wj = VEC_SCAL(w[j]);
280          int indexj = (index[j] + i)<<ldf;
281          a0 = VEC_FMA(wj, LOAD(fi + indexj              + k), a0);
282          a1 = VEC_FMA(wj, LOAD(fi + indexj + 1*VEC_SIZE + k), a1);
283        }
284        STORE(fo + (i<<ldf)              + k, a0);
285        STORE(fo + (i<<ldf) + 1*VEC_SIZE + k, a1);
286      }
287    }
288#endif
289
290#if DEPTH >= 1
291    for (; i < (rimap_inv_max[l] - unroll1 + 1) ; i += unroll1){
292      int k;
293      for(k = 0; k < (1<<ldf); k += VEC_SIZE){
294        register VEC_TYPE a0 = VEC_ZERO;
295        for(j = 0; j < n; j++) {
296#ifdef VEC_SCAL_LD
297          register VEC_TYPE wj = VEC_SCAL_LD(w + j);
298#else
299          register VEC_TYPE wj = VEC_SCAL(w[j]);
300#endif
301          int indexj = (index[j] + i)<<ldf;
302          a0 = VEC_FMA(wj, LOAD(fi + indexj + k), a0);
303        }
304        STORE(fo + (i<<ldf) + k, a0);
305      }
306    }
307#endif
308
309#if VEC_SIZE > 1
310    for (; i < rimap_inv_max[l]; i++){
311      int k;
312      for(k = 0; k < (1<<ldf); k++){
313        double a = 0.0;
314        for(j = 0; j < n; j++) a += w[j]*fi[((index[j] + i)<<ldf) + k];
315        fo[(i<<ldf) + k] = a;
316      }
317    }
318#endif
319
320  } /* l */
321
322}
323
324#undef LOAD
325#undef STORE
Note: See TracBrowser for help on using the repository browser.