Context Navigation

operate_inc.c @ 7019

Revision 7019, 10.7 KB checked in by dstrubbe, 3 years ago (diff)
) Now the configure script checks that the hardware actually supports AVX instructions. ) Removed some unnecessary semicolons.

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., 59 Temple Place - Suite 330, Boston, MA
17	02111-1307, USA.
18
19	$Id: operate_ri_vec.c 2146 2006-05-23 17:36:00Z xavier $
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.

Context Navigation

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

Download in other formats: