fix nasty lda bug in lange (#278)

* fix nasty lda bug in lange * switch to using blas amax instead
cvxgrp · Jun 27, 2024 · 7fc2a88 · 7fc2a88
1 parent 45988c9
commit 7fc2a88
Show file tree

Hide file tree

Showing 2 changed files with 21 additions and 3 deletions.
diff --git a/include/scs_blas.h b/include/scs_blas.h
@@ -17,8 +17,10 @@ extern "C" {
 /* single or double precision */
 #ifndef SFLOAT
 #define BLAS(x) d##x
+#define BLASI(x) id##x
 #else
 #define BLAS(x) s##x
+#define BLASI(x) is##x
 #endif
 #else
 /* this extra indirection is needed for BLASSUFFIX to work correctly as a
@@ -28,8 +30,10 @@ extern "C" {
 /* single or double precision */
 #ifndef SFLOAT
 #define BLAS(x) stitch__(d, x, BLASSUFFIX)
+#define BLASI(x) stitch__(id, x, BLASSUFFIX)
 #else
 #define BLAS(x) stitch__(s, x, BLASSUFFIX)
+#define BLASI(x) stitch__(is, x, BLASSUFFIX)
 #endif
 #endif
 

diff --git a/src/linalg.c b/src/linalg.c
@@ -102,12 +102,16 @@ extern "C" {
 scs_float BLAS(nrm2)(blas_int *n, const scs_float *x, blas_int *incx);
 scs_float BLAS(dot)(const blas_int *n, const scs_float *x, const blas_int *incx,
                     const scs_float *y, const blas_int *incy);
-scs_float BLAS(lange)(const char *norm, const blas_int *m, const blas_int *n,
-                      const scs_float *a, blas_int *lda, scs_float *work);
 void BLAS(axpy)(blas_int *n, const scs_float *a, const scs_float *x,
                 blas_int *incx, scs_float *y, blas_int *incy);
 void BLAS(scal)(const blas_int *n, const scs_float *sa, scs_float *sx,
                 const blas_int *incx);
+blas_int BLASI(amax)(blas_int *n, const scs_float *x, blas_int *incx);
+
+/* Possibly not working correctly on all platforms.
+scs_float BLAS(lange)(const char *norm, const blas_int *m, const blas_int *n,
+                      const scs_float *a, blas_int *lda, scs_float *work);
+*/
 
 #ifdef __cplusplus
 }
@@ -140,10 +144,20 @@ scs_float SCS(norm_2)(const scs_float *v, scs_int len) {
   return BLAS(nrm2)(&blen, v, &bone);
 }
 
+/* Possibly not working correctly on all platforms.
+scs_float SCS(norm_inf)(const scs_float *a, scs_int len) {
+  blas_int bone = 1;
+  blas_int blen = (blas_int)len;
+  return BLAS(lange)("Max", &blen, &bone, a, &blen, SCS_NULL);
+}
+*/
+
 scs_float SCS(norm_inf)(const scs_float *a, scs_int len) {
   blas_int bone = 1;
   blas_int blen = (blas_int)len;
-  return BLAS(lange)("Max", &blen, &bone, a, &bone, SCS_NULL);
+  scs_int idx = (scs_int)BLASI(amax)(&blen, a, &bone);
+  /* Returned idx is 1-based. */
+  return ABS(a[idx - 1]);
 }
 
 /* axpy a += sc*b */