BUG: ensure empty cholesky upper does not hang.

numpy/linalg/tests/test_regression.py

@@ -2,6 +2,8 @@
 """
 import warnings
 
+import pytest
+
 import numpy as np
 from numpy import linalg, arange, float64, array, dot, transpose
 from numpy.testing import (
@@ -146,3 +148,9 @@ def test_lstsq_complex_larger_rhs(self):
         u_lstsq, res, rank, sv = linalg.lstsq(G, b, rcond=None)
         # check results just in case
         assert_array_almost_equal(u_lstsq, u)
+
+    @pytest.mark.parametrize("upper", [True, False])
+    def test_cholesky_empty_array(self, upper):
+        # gh-25840 - upper=True hung before.
+        res = np.linalg.cholesky(np.zeros((0, 0)), upper=upper)
+        assert res.size == 0

numpy/linalg/umath_linalg.cpp

@@ -983,36 +983,6 @@ identity_matrix(typ *matrix, size_t n)
     }
 }
 
-         /* zero the undefined part in a upper/lower triangular matrix */
-          /* Note: matrix from fortran routine, so column-major order */
-
-template<typename typ>
-static inline void
-triu_matrix(typ *matrix, size_t n)
-{
-    size_t i, j;
-    for (i = 0; i < n-1; ++i) {
-        for (j = i+1; j < n; ++j) {
-            matrix[j] = numeric_limits<typ>::zero;
-        }
-        matrix += n;
-    }
-}
-
-template<typename typ>
-static inline void
-tril_matrix(typ *matrix, size_t n)
-{
-    size_t i, j;
-    matrix += n;
-    for (i = 1; i < n; ++i) {
-        for (j = 0; j < i; ++j) {
-            matrix[j] = numeric_limits<typ>::zero;
-        }
-        matrix += n;
-    }
-}
-
 /* -------------------------------------------------------------------------- */
                           /* Determinants */
 
@@ -1884,6 +1854,40 @@ struct POTR_PARAMS_t
 };
 
 
+         /* zero the undefined part in a upper/lower triangular matrix */
+          /* Note: matrix from fortran routine, so column-major order */
+
+template<typename typ>
+static inline void
+zero_lower_triangle(POTR_PARAMS_t<typ> *params)
+{
+    fortran_int n = params->N;
+    typ *matrix = params->A;
+    fortran_int i, j;
+    for (i = 0; i < n-1; ++i) {
+        for (j = i+1; j < n; ++j) {
+            matrix[j] = numeric_limits<typ>::zero;
+        }
+        matrix += n;
+    }
+}
+
+template<typename typ>
+static inline void
+zero_upper_triangle(POTR_PARAMS_t<typ> *params)
+{
+    fortran_int n = params->N;
+    typ *matrix = params->A;
+    fortran_int i, j;
+    matrix += n;
+    for (i = 1; i < n; ++i) {
+        for (j = 0; j < i; ++j) {
+            matrix[j] = numeric_limits<typ>::zero;
+        }
+        matrix += n;
+    }
+}
+
 static inline fortran_int
 call_potrf(POTR_PARAMS_t<fortran_real> *params)
 {
@@ -1983,10 +1987,10 @@ cholesky(char uplo, char **args, npy_intp const *dimensions, npy_intp const *ste
             not_ok = call_potrf(&params);
             if (!not_ok) {
                 if (uplo == 'L') {
-                    tril_matrix(params.A, params.N);
+                    zero_upper_triangle(&params);
                 }
                 else {
-                    triu_matrix(params.A, params.N);
+                    zero_lower_triangle(&params);
                 }
                 delinearize_matrix((ftyp*)args[1], params.A, &r_out);
             } else {