blob: 5f8fc2e587190ba5144665dcdfe1861e36634ae3 [file] [log] [blame]
Austin Schuh9a24b372018-01-28 16:12:29 -08001/**************************************************************************************************
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
30
31#if defined(LA_REFERENCE)
32
33
34
35void AXPY_LIBSTR(int m, REAL alpha, struct STRVEC *sx, int xi, struct STRVEC *sy, int yi, struct STRVEC *sz, int zi)
36 {
37 if(m<=0)
38 return;
39 int ii;
40 REAL *x = sx->pa + xi;
41 REAL *y = sy->pa + yi;
42 REAL *z = sz->pa + zi;
43 ii = 0;
44 for(; ii<m-3; ii+=4)
45 {
46 z[ii+0] = y[ii+0] + alpha*x[ii+0];
47 z[ii+1] = y[ii+1] + alpha*x[ii+1];
48 z[ii+2] = y[ii+2] + alpha*x[ii+2];
49 z[ii+3] = y[ii+3] + alpha*x[ii+3];
50 }
51 for(; ii<m; ii++)
52 z[ii+0] = y[ii+0] + alpha*x[ii+0];
53 return;
54 }
55
56
57
58// multiply two vectors and compute dot product
59REAL VECMULDOT_LIBSTR(int m, struct STRVEC *sx, int xi, struct STRVEC *sy, int yi, struct STRVEC *sz, int zi)
60 {
61 if(m<=0)
62 return 0.0;
63 REAL *x = sx->pa + xi;
64 REAL *y = sy->pa + yi;
65 REAL *z = sz->pa + zi;
66 int ii;
67 REAL dot = 0.0;
68 ii = 0;
69 for(; ii<m-3; ii+=4)
70 {
71 z[ii+0] = x[ii+0] * y[ii+0];
72 z[ii+1] = x[ii+1] * y[ii+1];
73 z[ii+2] = x[ii+2] * y[ii+2];
74 z[ii+3] = x[ii+3] * y[ii+3];
75 dot += z[ii+0] + z[ii+1] + z[ii+2] + z[ii+3];
76 }
77 for(; ii<m; ii++)
78 {
79 z[ii+0] = x[ii+0] * y[ii+0];
80 dot += z[ii+0];
81 }
82 return dot;
83 }
84
85
86
87// compute dot product of two vectors
88REAL DOT_LIBSTR(int m, struct STRVEC *sx, int xi, struct STRVEC *sy, int yi)
89 {
90 if(m<=0)
91 return 0.0;
92 REAL *x = sx->pa + xi;
93 REAL *y = sy->pa + yi;
94 int ii;
95 REAL dot = 0.0;
96 ii = 0;
97 for(; ii<m-3; ii+=4)
98 {
99 dot += x[ii+0] * y[ii+0];
100 dot += x[ii+1] * y[ii+1];
101 dot += x[ii+2] * y[ii+2];
102 dot += x[ii+3] * y[ii+3];
103 }
104 for(; ii<m; ii++)
105 {
106 dot += x[ii+0] * y[ii+0];
107 }
108 return dot;
109 }
110
111
112
113#elif defined(LA_BLAS)
114
115
116
117void AXPY_LIBSTR(int m, REAL alpha, struct STRVEC *sx, int xi, struct STRVEC *sy, int yi, struct STRVEC *sz, int zi)
118 {
119 if(m<=0)
120 return;
121 int i1 = 1;
122 REAL *x = sx->pa + xi;
123 REAL *y = sy->pa + yi;
124 REAL *z = sz->pa + zi;
125 if(y!=z)
126 COPY(&m, y, &i1, z, &i1);
127 AXPY(&m, &alpha, x, &i1, z, &i1);
128 return;
129 }
130
131
132
133// multiply two vectors and compute dot product
134REAL VECMULDOT_LIBSTR(int m, struct STRVEC *sx, int xi, struct STRVEC *sy, int yi, struct STRVEC *sz, int zi)
135 {
136 if(m<=0)
137 return 0.0;
138 REAL *x = sx->pa + xi;
139 REAL *y = sy->pa + yi;
140 REAL *z = sz->pa + zi;
141 int ii;
142 REAL dot = 0.0;
143 ii = 0;
144 for(; ii<m; ii++)
145 {
146 z[ii+0] = x[ii+0] * y[ii+0];
147 dot += z[ii+0];
148 }
149 return dot;
150 }
151
152
153
154// compute dot product of two vectors
155REAL DOT_LIBSTR(int m, struct STRVEC *sx, int xi, struct STRVEC *sy, int yi)
156 {
157 if(m<=0)
158 return 0.0;
159 REAL *x = sx->pa + xi;
160 REAL *y = sy->pa + yi;
161 int ii;
162 REAL dot = 0.0;
163 ii = 0;
164 for(; ii<m-3; ii+=4)
165 {
166 dot += x[ii+0] * y[ii+0];
167 dot += x[ii+1] * y[ii+1];
168 dot += x[ii+2] * y[ii+2];
169 dot += x[ii+3] * y[ii+3];
170 }
171 for(; ii<m; ii++)
172 {
173 dot += x[ii+0] * y[ii+0];
174 }
175 return dot;
176 }
177
178
179
180#else
181
182#error : wrong LA choice
183
184#endif
185
186