BUG: linalg: avoid a segfault in ?gesdd

linalg.svd with too large matrices may segfault if m*n > int_max on large-memory machines (on smaller memory machines it may fail with a MemoryError instead). The root cause is an integer overflow in indexing 2D arrays, deep in the LAPACK code. Thus, detect a possible error condition in the f2py wrapper.
scipy · Mar 28, 2024 · 16512f2 · 16512f2
1 parent 01ef73e
commit 16512f2
Show file tree

Hide file tree

Showing 3 changed files with 12 additions and 2 deletions.
diff --git a/scipy/linalg/flapack.pyf.src b/scipy/linalg/flapack.pyf.src
@@ -26,7 +26,14 @@
 
 python module _flapack
     usercode '''
+#ifdef HAVE_ILP64
+#define F_INT npy_int64
+#define F_INT_MAX NPY_MAX_INT64
+#else
 #define F_INT int
+#define F_INT_MAX INT_MAX
+#endif
+
 '''
 
 interface

diff --git a/scipy/linalg/flapack_64.pyf.src b/scipy/linalg/flapack_64.pyf.src
@@ -20,6 +20,7 @@ python module _flapack_64
 #include "blas64-prefix-defines.h"
 #endif
 #define F_INT npy_int64
+#define F_INT_MAX NPY_MAX_INT64
 '''
 
 interface

diff --git a/scipy/linalg/flapack_gen.pyf.src b/scipy/linalg/flapack_gen.pyf.src
@@ -307,7 +307,8 @@ subroutine <prefix2>gesdd(m,n,minmn,u0,u1,vt0,vt1,a,compute_uv,full_matrices,u,s
     integer intent(in),optional,check(full_matrices==0||full_matrices==1):: full_matrices = 1
     integer intent(hide),depend(a):: m = shape(a,0)
     integer intent(hide),depend(a):: n = shape(a,1)
-    integer intent(hide),depend(m,n):: minmn = MIN(m,n)
+    check(m <= sqrt(F_INT_MAX)/n)
+    integer intent(hide),depend(m,n) :: minmn = MIN(m,n)
     integer intent(hide),depend(compute_uv,minmn) :: u0 = (compute_uv?m:1)
     integer intent(hide),depend(compute_uv,minmn, full_matrices) :: u1 = (compute_uv?(full_matrices?m:minmn):1)
     integer intent(hide),depend(compute_uv,minmn, full_matrices) :: vt0 = (compute_uv?(full_matrices?n:minmn):1)
@@ -335,7 +336,8 @@ subroutine <prefix2>gesdd_lwork(m,n,minmn,u0,vt0,a,compute_uv,full_matrices,u,s,
     integer intent(in),optional,check(full_matrices==0||full_matrices==1):: full_matrices = 1
     integer intent(in) :: m
     integer intent(in) :: n
-    integer intent(hide),depend(m,n):: minmn = MIN(m,n)
+    check(m <= sqrt(F_INT_MAX)/n)
+    integer intent(hide),depend(m,n) :: minmn = MIN(m,n)
     integer intent(hide),depend(compute_uv,minmn) :: u0 = (compute_uv?m:1)
     integer intent(hide),depend(compute_uv,minmn, full_matrices) :: vt0 = (compute_uv?(full_matrices?n:minmn):1)
     <ftype2> intent(hide) :: a