Gitiles
Code Review Sign In
realtimeroboticsgroup.org / RealtimeRoboticsGroup / test / 9a24b370cbc07492fa6c205f1027367fa55b6742 / . / blas / d_blas1_lib4.c
blob: a4155a932be1b9948c985fcab7f6c90d7a7fc89e [file] [log] [blame]
Austin Schuh9a24b372018-01-28 16:12:29 -0800[diff] [blame^]1/**************************************************************************************************
2* *
3* This file is part of BLASFEO. *
4* *
5* BLASFEO -- BLAS For Embedded Optimization. *
6* Copyright (C) 2016-2017 by Gianluca Frison. *
7* Developed at IMTEK (University of Freiburg) under the supervision of Moritz Diehl. *
8* All rights reserved. *
9* *
10* HPMPC is free software; you can redistribute it and/or *
11* modify it under the terms of the GNU Lesser General Public *
12* License as published by the Free Software Foundation; either *
13* version 2.1 of the License, or (at your option) any later version. *
14* *
15* HPMPC is distributed in the hope that it will be useful, *
16* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
18* See the GNU Lesser General Public License for more details. *
19* *
20* You should have received a copy of the GNU Lesser General Public *
21* License along with HPMPC; if not, write to the Free Software *
22* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
23* *
24* Author: Gianluca Frison, giaf (at) dtu.dk *
25* gianluca.frison (at) imtek.uni-freiburg.de *
26* *
27**************************************************************************************************/
28
29#include <stdlib.h>
30#include <stdio.h>
31
32#if defined(TARGET_X64_INTEL_HASWELL) || defined(TARGET_X64_INTEL_SANDY_BRIDGE)
33#include <mmintrin.h>
34#include <xmmintrin.h> // SSE
35#include <emmintrin.h> // SSE2
36#include <pmmintrin.h> // SSE3
37#include <smmintrin.h> // SSE4
38#include <immintrin.h> // AVX
39#endif
40
41#include "../include/blasfeo_common.h"
42#include "../include/blasfeo_d_kernel.h"
43
44
45
46#if defined(LA_HIGH_PERFORMANCE)
47
48
49
50void daxpy_libstr(int m, double alpha, struct d_strvec *sx, int xi, struct d_strvec *sy, int yi, struct d_strvec *sz, int zi)
51 {
52
53 if(m<=0)
54 return;
55
56 double *x = sx->pa + xi;
57 double *y = sy->pa + yi;
58 double *z = sz->pa + zi;
59
60 int ii;
61
62#if defined(TARGET_X64_INTEL_HASWELL) || defined(TARGET_X64_INTEL_SANDY_BRIDGE)
63 __m256d
64 v_alpha, v_tmp,
65 v_x0, v_y0,
66 v_x1, v_y1;
67#endif
68
69 ii = 0;
70#if defined(TARGET_X64_INTEL_HASWELL) || defined(TARGET_X64_INTEL_SANDY_BRIDGE)
71 v_alpha = _mm256_broadcast_sd( &alpha );
72 for( ; ii<m-7; ii+=8)
73 {
74 v_x0 = _mm256_loadu_pd( &x[ii+0] );
75 v_x1 = _mm256_loadu_pd( &x[ii+4] );
76 v_y0 = _mm256_loadu_pd( &y[ii+0] );
77 v_y1 = _mm256_loadu_pd( &y[ii+4] );
78#if defined(TARGET_X64_INTEL_HASWELL)
79 v_y0 = _mm256_fmadd_pd( v_alpha, v_x0, v_y0 );
80 v_y1 = _mm256_fmadd_pd( v_alpha, v_x1, v_y1 );
81#else // sandy bridge
82 v_tmp = _mm256_mul_pd( v_alpha, v_x0 );
83 v_y0 = _mm256_add_pd( v_tmp, v_y0 );
84 v_tmp = _mm256_mul_pd( v_alpha, v_x1 );
85 v_y1 = _mm256_add_pd( v_tmp, v_y1 );
86#endif
87 _mm256_storeu_pd( &z[ii+0], v_y0 );
88 _mm256_storeu_pd( &z[ii+4], v_y1 );
89 }
90 for( ; ii<m-3; ii+=4)
91 {
92 v_x0 = _mm256_loadu_pd( &x[ii] );
93 v_y0 = _mm256_loadu_pd( &y[ii] );
94#if defined(TARGET_X64_INTEL_HASWELL)
95 v_y0 = _mm256_fmadd_pd( v_alpha, v_x0, v_y0 );
96#else // sandy bridge
97 v_tmp = _mm256_mul_pd( v_alpha, v_x0 );
98 v_y0 = _mm256_add_pd( v_tmp, v_y0 );
99#endif
100 _mm256_storeu_pd( &z[ii], v_y0 );
101 }
102#else
103 for( ; ii<m-3; ii+=4)
104 {
105 z[ii+0] = y[ii+0] + alpha*x[ii+0];
106 z[ii+1] = y[ii+1] + alpha*x[ii+1];
107 z[ii+2] = y[ii+2] + alpha*x[ii+2];
108 z[ii+3] = y[ii+3] + alpha*x[ii+3];
109 }
110#endif
111 for( ; ii<m; ii++)
112 {
113 z[ii+0] = y[ii+0] + alpha*x[ii+0];
114 }
115
116 return;
117 }
118
119
120
121// multiply two vectors and compute dot product
122double dvecmuldot_libstr(int m, struct d_strvec *sx, int xi, struct d_strvec *sy, int yi, struct d_strvec *sz, int zi)
123 {
124
125 if(m<=0)
126 return 0.0;
127
128 double *x = sx->pa + xi;
129 double *y = sy->pa + yi;
130 double *z = sz->pa + zi;
131 int ii;
132 double dot = 0.0;
133#if defined(TARGET_X64_INTEL_HASWELL) || defined(TARGET_X64_INTEL_SANDY_BRIDGE)
134 __m128d
135 u_tmp, u_dot;
136 __m256d
137 v_tmp,
138 v_x0, v_y0, v_z0;
139
140 v_tmp = _mm256_setzero_pd();
141#endif
142
143 ii = 0;
144
145#if defined(TARGET_X64_INTEL_HASWELL) || defined(TARGET_X64_INTEL_SANDY_BRIDGE)
146 for(; ii<m-3; ii+=4)
147 {
148 v_x0 = _mm256_loadu_pd( &x[ii+0] );
149 v_y0 = _mm256_loadu_pd( &y[ii+0] );
150 v_z0 = _mm256_mul_pd( v_x0, v_y0 );
151 _mm256_storeu_pd( &z[ii+0], v_z0 );
152 v_tmp = _mm256_add_pd( v_tmp, v_z0 );
153 }
154#endif
155 for(; ii<m; ii++)
156 {
157 z[ii+0] = x[ii+0] * y[ii+0];
158 dot += z[ii+0];
159 }
160#if defined(TARGET_X64_INTEL_HASWELL) || defined(TARGET_X64_INTEL_SANDY_BRIDGE)
161 // dot product
162 u_tmp = _mm_add_pd( _mm256_castpd256_pd128( v_tmp ), _mm256_extractf128_pd( v_tmp, 0x1 ) );
163 u_tmp = _mm_hadd_pd( u_tmp, u_tmp);
164 u_dot = _mm_load_sd( &dot );
165 u_dot = _mm_add_sd( u_dot, u_tmp );
166 _mm_store_sd( &dot, u_dot );
167#endif
168 return dot;
169 }
170
171
172
173// compute dot product of two vectors
174double ddot_libstr(int m, struct d_strvec *sx, int xi, struct d_strvec *sy, int yi)
175 {
176
177 if(m<=0)
178 return 0.0;
179
180 double *x = sx->pa + xi;
181 double *y = sy->pa + yi;
182 int ii;
183 double dot = 0.0;
184
185#if defined(TARGET_X64_INTEL_HASWELL) || defined(TARGET_X64_INTEL_SANDY_BRIDGE)
186 __m128d
187 u_dot0, u_x0, u_y0, u_tmp;
188 __m256d
189 v_dot0, v_dot1, v_x0, v_x1, v_y0, v_y1, v_tmp;
190
191 v_dot0 = _mm256_setzero_pd();
192 v_dot1 = _mm256_setzero_pd();
193 u_dot0 = _mm_setzero_pd();
194
195 ii = 0;
196 for(; ii<m-7; ii+=8)
197 {
198 v_x0 = _mm256_loadu_pd( &x[ii+0] );
199 v_x1 = _mm256_loadu_pd( &x[ii+4] );
200 v_y0 = _mm256_loadu_pd( &y[ii+0] );
201 v_y1 = _mm256_loadu_pd( &y[ii+4] );
202#if defined(TARGET_X64_INTEL_HASWELL)
203 v_dot0 = _mm256_fmadd_pd( v_x0, v_y0, v_dot0 );
204 v_dot1 = _mm256_fmadd_pd( v_x1, v_y1, v_dot1 );
205#else // sandy bridge
206 v_tmp = _mm256_mul_pd( v_x0, v_y0 );
207 v_dot0 = _mm256_add_pd( v_dot0, v_tmp );
208 v_tmp = _mm256_mul_pd( v_x1, v_y1 );
209 v_dot1 = _mm256_add_pd( v_dot1, v_tmp );
210#endif
211 }
212 for(; ii<m-3; ii+=4)
213 {
214 v_x0 = _mm256_loadu_pd( &x[ii+0] );
215 v_y0 = _mm256_loadu_pd( &y[ii+0] );
216#if defined(TARGET_X64_INTEL_HASWELL)
217 v_dot0 = _mm256_fmadd_pd( v_x0, v_y0, v_dot0 );
218#else // sandy bridge
219 v_tmp = _mm256_mul_pd( v_x0, v_y0 );
220 v_dot0 = _mm256_add_pd( v_dot0, v_tmp );
221#endif
222 }
223 for(; ii<m; ii++)
224 {
225 u_x0 = _mm_load_sd( &x[ii+0] );
226 u_y0 = _mm_load_sd( &y[ii+0] );
227#if defined(TARGET_X64_INTEL_HASWELL)
228 u_dot0 = _mm_fmadd_sd( u_x0, u_y0, u_dot0 );
229#else // sandy bridge
230 u_tmp = _mm_mul_sd( u_x0, u_y0 );
231 u_dot0 = _mm_add_sd( u_dot0, u_tmp );
232#endif
233 }
234 // reduce
235 v_dot0 = _mm256_add_pd( v_dot0, v_dot1 );
236 u_tmp = _mm_add_pd( _mm256_castpd256_pd128( v_dot0 ), _mm256_extractf128_pd( v_dot0, 0x1 ) );
237 u_tmp = _mm_hadd_pd( u_tmp, u_tmp);
238 u_dot0 = _mm_add_sd( u_dot0, u_tmp );
239 _mm_store_sd( &dot, u_dot0 );
240#else // no haswell, no sandy bridge
241 ii = 0;
242 for(; ii<m-3; ii+=4)
243 {
244 dot += x[ii+0] * y[ii+0];
245 dot += x[ii+1] * y[ii+1];
246 dot += x[ii+2] * y[ii+2];
247 dot += x[ii+3] * y[ii+3];
248 }
249 for(; ii<m; ii++)
250 {
251 dot += x[ii+0] * y[ii+0];
252 }
253#endif // haswell, sandy bridge
254 return dot;
255 }
256
257
258
259#else
260
261#error : wrong LA choice
262
263#endif