1 #include "linalg.h"
2 #include "scs_blas.h"
3 #include <math.h>
4
5 /* these routines do not have BLAS implementations (that I can find at least) */
6
SCS(norm_diff)7 scs_float SCS(norm_diff)(const scs_float *a, const scs_float *b, scs_int len) {
8 scs_float nm_diff = 0.0, tmp;
9 scs_int i;
10 for (i = 0; i < len; ++i) {
11 tmp = (a[i] - b[i]);
12 nm_diff += tmp * tmp;
13 }
14 return SQRTF(nm_diff);
15 }
16
SCS(norm_inf_diff)17 scs_float SCS(norm_inf_diff)(const scs_float *a, const scs_float *b,
18 scs_int len) {
19 scs_float tmp, max = 0.0;
20 scs_int i;
21 for (i = 0; i < len; ++i) {
22 tmp = ABS(a[i] - b[i]);
23 if (tmp > max) {
24 max = tmp;
25 }
26 }
27 return max;
28 }
29
30 #ifndef USE_LAPACK
31 /* Self-rolled basic linear algebra routines */
32
33 /* a *= b */
SCS(scale_array)34 void SCS(scale_array)(scs_float *a, const scs_float b, scs_int len) {
35 scs_int i;
36 for (i = 0; i < len; ++i)
37 a[i] *= b;
38 }
39
40 /* x'*y */
SCS(dot)41 scs_float SCS(dot)(const scs_float *x, const scs_float *y, scs_int len) {
42 scs_int i;
43 scs_float ip = 0.0;
44 for (i = 0; i < len; ++i) {
45 ip += x[i] * y[i];
46 }
47 return ip;
48 }
49
50 /* ||v||_2^2 */
SCS(norm_sq)51 scs_float SCS(norm_sq)(const scs_float *v, scs_int len) {
52 scs_int i;
53 scs_float nmsq = 0.0;
54 for (i = 0; i < len; ++i) {
55 nmsq += v[i] * v[i];
56 }
57 return nmsq;
58 }
59
60 /* ||v||_2 */
SCS(norm_2)61 scs_float SCS(norm_2)(const scs_float *v, scs_int len) {
62 return SQRTF(SCS(norm_sq)(v, len));
63 }
64
SCS(norm_inf)65 scs_float SCS(norm_inf)(const scs_float *a, scs_int len) {
66 scs_float tmp, max = 0.0;
67 scs_int i;
68 for (i = 0; i < len; ++i) {
69 tmp = ABS(a[i]);
70 if (tmp > max) {
71 max = tmp;
72 }
73 }
74 return max;
75 }
76
77 /* axpy a += sc*b */
SCS(add_scaled_array)78 void SCS(add_scaled_array)(scs_float *a, const scs_float *b, scs_int n,
79 const scs_float sc) {
80 scs_int i;
81 for (i = 0; i < n; ++i) {
82 a[i] += sc * b[i];
83 }
84 }
85
86 #else
87 /* If we have BLAS / LAPACK we may as well use them */
88
89 #ifdef __cplusplus
90 extern "C" {
91 #endif
92
93 scs_float BLAS(nrm2)(blas_int *n, const scs_float *x, blas_int *incx);
94 scs_float BLAS(dot)(const blas_int *n, const scs_float *x, const blas_int *incx,
95 const scs_float *y, const blas_int *incy);
96 scs_float BLAS(lange)(const char *norm, const blas_int *m, const blas_int *n,
97 const scs_float *a, blas_int *lda, scs_float *work);
98 void BLAS(axpy)(blas_int *n, const scs_float *a, const scs_float *x,
99 blas_int *incx, scs_float *y, blas_int *incy);
100 void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx,
101 const blas_int *incx);
102
103 #ifdef __cplusplus
104 }
105 #endif
106
107 /* a *= b */
SCS(scale_array)108 void SCS(scale_array)(scs_float *a, const scs_float b, scs_int len) {
109 blas_int bone = 1;
110 blas_int blen = (blas_int)len;
111 BLAS(scal)(&blen, &b, a, &bone);
112 }
113
114 /* x'*y */
SCS(dot)115 scs_float SCS(dot)(const scs_float *x, const scs_float *y, scs_int len) {
116 blas_int bone = 1;
117 blas_int blen = (blas_int)len;
118 return BLAS(dot)(&blen, x, &bone, y, &bone);
119 }
120
121 /* ||v||_2^2 */
SCS(norm_sq)122 scs_float SCS(norm_sq)(const scs_float *v, scs_int len) {
123 scs_float nrm = SCS(norm_2)(v, len);
124 return nrm * nrm;
125 }
126
127 /* ||v||_2 */
SCS(norm_2)128 scs_float SCS(norm_2)(const scs_float *v, scs_int len) {
129 blas_int bone = 1;
130 blas_int blen = (blas_int)len;
131 return BLAS(nrm2)(&blen, v, &bone);
132 }
133
SCS(norm_inf)134 scs_float SCS(norm_inf)(const scs_float *a, scs_int len) {
135 blas_int bone = 1;
136 blas_int blen = (blas_int)len;
137 return BLAS(lange)("Max", &blen, &bone, a, &bone, SCS_NULL);
138 }
139
140 /* axpy a += sc*b */
SCS(add_scaled_array)141 void SCS(add_scaled_array)(scs_float *a, const scs_float *b, scs_int len,
142 const scs_float sc) {
143 blas_int bone = 1;
144 blas_int blen = (blas_int)len;
145 BLAS(axpy)(&blen, &sc, b, &bone, a, &bone);
146 }
147
148 #endif
149