Line data Source code
1 : /**
2 : * @file lin_alg.c
3 : * @author Ben Hermans
4 : * @brief Linear algebra with vectors.
5 : * @details Common operations, such as vector products, infinity norm, elementwise
6 : * add/product/division/max etc. are included in this file.
7 : */
8 :
9 : #include "lin_alg.h"
10 :
11 282 : c_float* vec_copy(const c_float *a, size_t n) {
12 : c_float *b;
13 : size_t i;
14 :
15 282 : b = qpalm_malloc(n * sizeof(c_float));
16 :
17 1439 : for (i = 0; i < n; i++) {
18 1157 : b[i] = a[i];
19 : }
20 :
21 282 : return b;
22 : }
23 :
24 4629 : void prea_vec_copy(const c_float *a, c_float *b, size_t n) {
25 : size_t i;
26 :
27 25468 : for (i = 0; i < n; i++) {
28 20839 : b[i] = a[i];
29 : }
30 4629 : }
31 :
32 430 : void prea_int_vec_copy(const c_int *a, c_int *b, size_t n) {
33 : size_t i;
34 :
35 2657 : for (i = 0; i < n; i++) {
36 2227 : b[i] = a[i];
37 : }
38 430 : }
39 :
40 852 : void vec_set_scalar(c_float *a, c_float sc, size_t n) {
41 : size_t i;
42 :
43 4563 : for (i = 0; i < n; i++) {
44 3711 : a[i] = sc;
45 : }
46 852 : }
47 :
48 187 : void vec_set_scalar_int(c_int *a, c_int sc, size_t n) {
49 : size_t i;
50 :
51 1017 : for (i = 0; i < n; i++) {
52 830 : a[i] = sc;
53 : }
54 187 : }
55 :
56 2255 : void vec_self_mult_scalar(c_float *a, c_float sc, size_t n) {
57 : size_t i;
58 :
59 13844 : for (i = 0; i < n; i++) {
60 11589 : a[i] *= sc;
61 : }
62 2255 : }
63 :
64 6 : void vec_mult_scalar(const c_float *a, c_float sc, c_float *b, size_t n) {
65 : size_t i;
66 :
67 30 : for (i = 0; i < n; i++) {
68 24 : b[i] = sc*a[i];
69 : }
70 6 : }
71 :
72 2013 : c_float vec_prod(const c_float *a, const c_float *b, size_t n) {
73 2013 : c_float prod = 0.0;
74 2013 : size_t i = 0;
75 :
76 2013 : if(n >= 4) {
77 2564 : for (; i <= n-4; i+=4) {
78 1432 : prod += (a[i]*b[i] + a[i+1]*b[i+1] + a[i+2]*b[i+2] + a[i+3]*b[i+3]);
79 : }
80 : }
81 4401 : for (; i < n; i++) {
82 2388 : prod += a[i] * b[i];
83 : }
84 :
85 2013 : return prod;
86 : }
87 :
88 21 : c_float vec_norm_two(const c_float *a, size_t n) {
89 21 : return c_sqrt(vec_prod(a, a, n));
90 : }
91 :
92 10542 : void vec_ew_prod(const c_float *a, const c_float *b, c_float *c, size_t n) {
93 : size_t i;
94 :
95 60119 : for (i = 0; i < n; i++) {
96 49577 : c[i] = a[i] * b[i];
97 : }
98 10542 : }
99 :
100 :
101 1321 : void vec_ew_div(const c_float *a, const c_float *b, c_float *c, size_t n) {
102 : size_t i;
103 :
104 10307 : for (i = 0; i < n; i++) {
105 8986 : c[i] = a[i] / b[i];
106 : }
107 1321 : }
108 :
109 :
110 11492 : void vec_add_scaled(const c_float *a, const c_float *b, c_float *c, c_float sc, size_t n) {
111 : size_t i;
112 :
113 62865 : for (i = 0; i < n; i++) {
114 51373 : c[i] = a[i] + sc * b[i];
115 : }
116 11492 : }
117 :
118 411 : void vec_mult_add_scaled(c_float *a, const c_float *b, c_float sc1, c_float sc2, size_t n) {
119 : size_t i;
120 :
121 2391 : for (i = 0; i < n; i++) {
122 1980 : a[i] = sc1*a[i] + sc2*b[i];
123 : }
124 411 : }
125 :
126 9018 : c_float vec_norm_inf(const c_float *a, size_t n) {
127 9018 : register size_t j = 0;
128 9018 : register c_float s0 = 0.;
129 9018 : register c_float s1 = 0.;
130 9018 : register c_float s2 = 0.;
131 9018 : register c_float s3 = 0.;
132 9018 : register c_float max0 = 0.;
133 9018 : register c_float max1 = 0.;
134 9018 : register c_float max2 = 0.;
135 9018 : register c_float max3 = 0.;
136 9018 : const size_t block_size = 4;
137 9018 : const size_t block_len = n >> 2;
138 9018 : const size_t remaining = n % block_size; /*Initializing four blocks for
139 : * efficient implementation of
140 : * inner product using registers */
141 :
142 15930 : while (j < block_len * block_size) {
143 6912 : s0 = c_absval(a[j]); max0 = s0 > max0 ? s0 : max0;
144 6912 : s1 = c_absval(a[j+1]); max1 = s1 > max1 ? s1 : max1;
145 6912 : s2 = c_absval(a[j+2]); max2 = s2 > max2 ? s2 : max2;
146 6912 : s3 = c_absval(a[j+3]); max3 = s3 > max3 ? s3 : max3;
147 6912 : j+=4;
148 : }
149 :
150 9018 : max0 = max0 > max1 ? max0 : max1;
151 9018 : max0 = max0 > max2 ? max0 : max2;
152 9018 : max0 = max0 > max3 ? max0 : max3;
153 9018 : j = block_size * block_len;
154 : /*Taking contribution from the last terms
155 : * that were not included in the block*/
156 9018 : switch (remaining) {
157 2031 : case 3:
158 2031 : max0 = max0 > c_absval(a[j + 2]) ? max0 : c_absval(a[j + 2]);
159 2031 : max0 = max0 > c_absval(a[j + 1]) ? max0 : c_absval(a[j + 1]);
160 2031 : max0 = max0 > c_absval(a[j]) ? max0 : c_absval(a[j]);
161 2031 : break;
162 2543 : case 2:
163 2543 : max0 = max0 > c_absval(a[j + 1]) ? max0 : c_absval(a[j + 1]);
164 2543 : max0 = max0 > c_absval(a[j]) ? max0 : c_absval(a[j]);
165 2543 : break;
166 1877 : case 1:
167 1877 : max0 = max0 > c_absval(a[j]) ? max0 : c_absval(a[j]);
168 1877 : break;
169 2567 : case 0:
170 2567 : break;
171 : }
172 9018 : return max0;
173 :
174 : }
175 :
176 1139 : void vec_ew_recipr(const c_float *a, c_float *b, size_t n) {
177 : size_t i;
178 :
179 6346 : for (i = 0; i < n; i++) {
180 5207 : b[i] = (c_float)1.0 / a[i];
181 : }
182 1139 : }
183 :
184 1 : void vec_ew_max_vec(const c_float *a, const c_float *b, c_float *c, size_t n) {
185 : size_t i;
186 :
187 4 : for (i = 0; i < n; i++) {
188 3 : c[i] = c_max(a[i], b[i]);
189 : }
190 1 : }
191 :
192 1 : void vec_ew_min_vec(const c_float *a, const c_float *b, c_float *c, size_t n) {
193 : size_t i;
194 :
195 4 : for (i = 0; i < n; i++) {
196 3 : c[i] = c_min(a[i], b[i]);
197 : }
198 1 : }
199 :
200 892 : void vec_ew_mid_vec(const c_float *a, const c_float *bmin, const c_float *bmax, c_float *c, size_t n) {
201 : size_t i;
202 :
203 5096 : for (i = 0; i < n; i++) {
204 4204 : c[i] = c_max(bmin[i], c_min(a[i], bmax[i]));
205 : }
206 892 : }
207 :
208 1009 : void vec_ew_sqrt(const c_float *a, c_float *b, size_t n) {
209 : size_t i;
210 :
211 5667 : for (i = 0; i < n; i++) {
212 4658 : b[i] = c_sqrt(a[i]);
213 : }
214 1009 : }
215 :
216 : /* MATRIX FUNCTIONS ----------------------------------------------------------*/
217 :
218 : /* Moved to solver_interface.c*/
219 :
|
