ENH: linalg: add Python wrapper of ?gtcon (scipy#21328)

* ENH: linalg: add Python wrapper of ?gtcon

* MAINT: linalg: wrap both signatures of gtcon; fix test

* TST: linalg: accommodate old versions of NumPy

Author: mdhaber

scipy/linalg/flapack_gen_tri.pyf.src

@@ -66,6 +66,57 @@ subroutine <prefix>gttrs(trans, n, nrhs, dl, d, du, du2, ipiv, b, ldb, info)
 end subroutine gttrs
 
 
+subroutine <prefix2>gtcon(norm,n,dl,d,du,du2,ipiv,anorm,rcond,work,iwork,info)
+    ! ?GTCON estimates the reciprocal of the condition number of a real
+    ! tridiagonal matrix A using the LU factorization as computed by
+    ! ?GTTRF.
+    ! An estimate is obtained for norm(inv(A)), and the reciprocal of the
+    ! condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
+    threadsafe
+    callstatement (*f2py_func)(norm,&n,dl,d,du,du2,ipiv,&anorm,&rcond,work,iwork,&info)
+    callprotoargument char*, F_INT*, <ctype2>*, <ctype2>*, <ctype2>*, <ctype2>*, F_INT*, <ctype2>*, <ctype2>*, <ctype2>*, F_INT*, F_INT*
+
+    character optional, intent(in) :: norm = '1'
+    integer intent(hide), depend(d) :: n = max(3, len(d))
+    <ftype2> intent(in), depend(n), dimension(n - 1) :: dl
+    <ftype2> intent(in), dimension(n) :: d
+    <ftype2> intent(in), depend(n), dimension(n - 1) :: du
+    <ftype2> intent(in), depend(n), dimension(n - 2) :: du2
+    integer intent(in), depend(n), dimension(n) :: ipiv
+    <ftype2> intent(in) :: anorm
+    <ftype2> intent(out) :: rcond
+    <ftype2> intent(hide, cache), dimension(2*n), depend(n) :: work
+    integer intent(hide, cache), dimension(n), depend(n) :: iwork
+    integer intent(out) :: info
+
+end subroutine <prefix2>gtcon
+
+
+subroutine <prefix2c>gtcon(norm,n,dl,d,du,du2,ipiv,anorm,rcond,work,info)
+    ! ?GTCON estimates the reciprocal of the condition number of a real
+    ! tridiagonal matrix A using the LU factorization as computed by
+    ! ?GTTRF.
+    ! An estimate is obtained for norm(inv(A)), and the reciprocal of the
+    ! condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
+    threadsafe
+    callstatement (*f2py_func)(norm,&n,dl,d,du,du2,ipiv,&anorm,&rcond,work,&info)
+    callprotoargument char*, F_INT*, <ctype2c>*, <ctype2c>*, <ctype2c>*, <ctype2c>*, F_INT*, <ctype2>*, <ctype2>*, <ctype2c>*, F_INT*
+
+    character optional, intent(in) :: norm = '1'
+    integer intent(hide), depend(d) :: n = max(3, len(d))
+    <ftype2c> intent(in), depend(n), dimension(n - 1) :: dl
+    <ftype2c> intent(in), dimension(n) :: d
+    <ftype2c> intent(in), depend(n), dimension(n - 1) :: du
+    <ftype2c> intent(in), depend(n), dimension(n - 2) :: du2
+    integer intent(in), depend(n), dimension(n) :: ipiv
+    <ftype2> intent(in) :: anorm
+    <ftype2> intent(out) :: rcond
+    <ftype2c> intent(hide, cache), dimension(2*n), depend(n) :: work
+    integer intent(out) :: info
+
+end subroutine <prefix2c>gtcon
+
+
 subroutine <prefix2>gtsvx(fact,trans,n,nrhs,dl,d,du,dlf,df,duf,du2,ipiv,b,ldb,x,ldx,rcond,ferr,berr,work,iwork,info)
     ! ?GTSVX uses the LU factorization to compute the solution to a real
     ! system of linear equations A * X = B or A**T * X = B,

scipy/linalg/lapack.py

@@ -793,6 +793,11 @@
    cgttrs
    zgttrs
 
+   sgtcon
+   dgtcon
+   cgtcon
+   zgtcon
+
    stpqrt
    dtpqrt
    ctpqrt

scipy/linalg/tests/test_lapack.py

@@ -2150,6 +2150,34 @@ def test_gttrf_gttrs_NAG_f07cdf_f07cef_f07crf_f07csf(du, d, dl, du_exp, d_exp,
     assert_allclose(x_gttrs, x)
 
 
+@pytest.mark.parametrize('dtype', DTYPES)
+@pytest.mark.parametrize('norm', ['1', 'I', 'O'])
+@pytest.mark.parametrize('n', [3, 10])
+def test_gtcon(dtype, norm, n):
+    rng = np.random.default_rng(23498324)
+
+    d = rng.random(n) + rng.random(n)*1j
+    dl = rng.random(n - 1) + rng.random(n - 1)*1j
+    du = rng.random(n - 1) + rng.random(n - 1)*1j
+    A = np.diag(d) + np.diag(dl, -1) + np.diag(du, 1)
+    if np.issubdtype(dtype, np.floating):
+        A, d, dl, du = A.real, d.real, dl.real, du.real
+    A, d, dl, du = A.astype(dtype), d.astype(dtype), dl.astype(dtype), du.astype(dtype)
+
+    anorm = np.abs(A).sum(axis=0).max()
+
+    gttrf, gtcon = get_lapack_funcs(('gttrf', 'gtcon'), (A,))
+    dl, d, du, du2, ipiv, info = gttrf(dl, d, du)
+    res, _ = gtcon(dl, d, du, du2, ipiv, anorm, norm=norm)
+
+    gecon, getrf = get_lapack_funcs(('gecon', 'getrf'), (A,))
+    lu, ipvt, info = getrf(A)
+    ref, _ = gecon(lu, anorm, norm=norm)
+
+    rtol = np.finfo(dtype).eps**0.75
+    assert_allclose(res, ref, rtol=rtol)
+
+
 @pytest.mark.parametrize('dtype', DTYPES)
 @pytest.mark.parametrize('shape', [(3, 7), (7, 3), (2**18, 2**18)])
 def test_geqrfp_lwork(dtype, shape):