Repository: scipy/scipy
Pull request: ENH: linalg.lapack: wrap ?sytri and ?hetri
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
https://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes #21764
#### What does this implement/fix?
<!--Please explain your changes.-->
As requested in the above issue wraps the functions from https://www.netlib.org/lapack/explore-html/da/dfa/group__hetri.html
#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
ENH: linalg.lapack: add symmetric solvers
### Describe your issue.

While the positive definite variants have all 3 functions present (trf, trs, tri), the hermitian and symmetric variants are missing 2 functions (trs, tri)

This is true for both the double precision and single precision versions.

### Reproducing Code Example

```python
from scipy.linalg import lapack

for name in ['zpotrf', 'zpotrs', 'zpotri', 'zhetrf', 'zhetrs', 'zhetri', 'dsytrf', 'dsytrs', 'dsytri']:
    print(name, hasattr(lapack, name))
```


### Error message

```shell
zpotrf True
zpotrs True
zpotri True
zhetrf True
zhetrs False
zhetri False
dsytrf True
dsytrs False
dsytri False
```
```


### SciPy/NumPy/Python version and system information

```shell
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.0.11",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.16.1",
      "include directory": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\pip-build-env-w8i5lw_j\\overlay\\Lib\\site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "scipy-openblas",
      "found": true,
      "version": "0.3.27.dev",
      "detection method": "pkgconfig",
      "include directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/include",
      "lib directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/lib",
      "openblas configuration": "OpenBLAS 0.3.27.dev DYNAMIC_ARCH NO_AFFINITY Zen MAX_THREADS=24",
      "pc file directory": "D:/a/scipy/scipy"
    },
    "lapack": {
      "name": "scipy-openblas",
      "found": true,
      "version": "0.3.27.dev",
      "detection method": "pkgconfig",
      "include directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/include",
      "lib directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/lib",
      "openblas configuration": "OpenBLAS 0.3.27.dev DYNAMIC_ARCH NO_AFFINITY Zen MAX_THREADS=24",
      "pc file directory": "D:/a/scipy/scipy"
    },
    "pybind11": {
      "name": "pybind11",
      "version": "2.12.0",
      "detection method": "config-tool",
      "include directory": "unknown"
    }
  },
  "Python Information": {
    "path": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\cibw-run-_mbvtoth\\cp310-win_amd64\\build\\venv\\Scripts\\python.exe",
    "version": "3.10"
  }
}
```

Hints from issue discussion:

Commit message: ENH: linalg: wrap ?hetri

Patch to review:
diff --git a/scipy/linalg/flapack_sym_herm.pyf.src b/scipy/linalg/flapack_sym_herm.pyf.src
index d0decd33c92e..c3d0db143451 100644
--- a/scipy/linalg/flapack_sym_herm.pyf.src
+++ b/scipy/linalg/flapack_sym_herm.pyf.src
@@ -1410,3 +1410,40 @@ subroutine <prefix2c>heequb(lower,n,a,lda,s,scond,amax,work,info)
     integer intent(out) :: info
 
 end subroutine <prefix2c>heequb
+
+
+subroutine <prefix>sytri(lower,n,a,lda,ipiv,work,info)
+    ! ?SYTRI computes the inverse of a symmetric indefinite matrix
+    ! A using the factorization A = U*D*U**T or A = L*D*L**T computed by
+    ! ?SYTRF.
+    threadsafe
+    callstatement (*f2py_func)((lower?"L":"U"),&n,a,&lda,ipiv,work,&info)
+    callprotoargument char*, F_INT*, <ctype>*, F_INT*, F_INT*, <ctype>*, F_INT*
+
+    integer optional intent(in), check(lower==0||lower==1) :: lower = 0
+    integer intent(hide), depend(a) :: n=shape(a,1)
+    <ftype> dimension(n,n),intent(in,out,copy,out=inv_a) :: a
+    integer depend(n),intent(hide) :: lda = max(n,1)
+    integer depend(n),dimension(n),intent(in) :: ipiv
+    <ftype> depend(n),dimension(n),intent(hide,cache) :: work
+    integer intent(out):: info
+end subroutine <prefix>sytri
+
+
+subroutine <prefix2c>hetri(lower,n,a,lda,ipiv,work,info)
+    ! ?HETRI computes the inverse of a complex Hermitian indefinite matrix
+    ! A using the factorization A = U*D*U**H or A = L*D*L**H computed by
+    ! ?HETRF.
+    threadsafe
+    callstatement (*f2py_func)((lower?"L":"U"),&n,a,&lda,ipiv,work,&info)
+    callprotoargument char*, F_INT*, <ctype2c>*, F_INT*, F_INT*, <ctype2c>*, F_INT*
+
+    integer optional intent(in), check(lower==0||lower==1) :: lower = 0
+    integer intent(hide), depend(a) :: n=shape(a,1)
+    <ftype2c> dimension(n,n),intent(in,out,copy,out=inv_a) :: a
+    integer depend(n),intent(hide) :: lda = max(n,1)
+    integer depend(n),dimension(n),intent(in) :: ipiv
+    <ftype2c> depend(n),dimension(n),intent(hide,cache) :: work
+    integer intent(out):: info
+end subroutine <prefix2c>hetri
+
diff --git a/scipy/linalg/lapack.py b/scipy/linalg/lapack.py
index ffce33ffd105..6559a2fef677 100644
--- a/scipy/linalg/lapack.py
+++ b/scipy/linalg/lapack.py
@@ -350,6 +350,9 @@
    chetrf_lwork
    zhetrf_lwork
 
+   chetri
+   zhetri
+
    chetrs
    zhetrs
 
@@ -673,6 +676,11 @@
    csytrf_lwork
    zsytrf_lwork
 
+   ssytri
+   dsytri
+   csytri
+   zsytri
+
    ssytrs
    dsytrs
    csytrs
diff --git a/scipy/linalg/tests/test_lapack.py b/scipy/linalg/tests/test_lapack.py
index f32ce213a5aa..48157870aa65 100644
--- a/scipy/linalg/tests/test_lapack.py
+++ b/scipy/linalg/tests/test_lapack.py
@@ -3487,6 +3487,33 @@ def test_sy_hetrs(mtype, dtype, lower):
     assert_allclose(A@x, b, atol=100*n*eps)
 
 
+@pytest.mark.parametrize('mtype', ['sy', 'he'])  # matrix type
+@pytest.mark.parametrize('dtype', DTYPES)
+@pytest.mark.parametrize('lower', (0, 1))
+def test_sy_he_tri(dtype, lower, mtype):
+    if mtype == 'he' and dtype in REAL_DTYPES:
+        pytest.skip("hetri not for real dtypes.")
+    rng = np.random.default_rng(1723059677121834)
+    n = 20
+    A = rng.random((n, n)) + rng.random((n, n))*1j
+    if np.issubdtype(dtype, np.floating):
+        A = A.real
+    A = A.astype(dtype)
+    A = A + A.T if mtype == 'sy' else A + A.conj().T
+    names = f'{mtype}trf', f'{mtype}tri'
+    trf, tri = get_lapack_funcs(names, dtype=dtype)
+    ldu, ipiv, info = trf(A, lower=lower)
+    assert info == 0
+    A_inv, info = tri(a=ldu, ipiv=ipiv, lower=lower)
+    assert info == 0
+    eps = np.finfo(dtype).eps
+    ref = np.linalg.inv(A)
+    if lower:
+        assert_allclose(np.tril(A_inv), np.tril(ref), atol=100*n*eps)
+    else:
+        assert_allclose(np.triu(A_inv), np.triu(ref), atol=100*n*eps)
+
+
 @pytest.mark.parametrize('norm', list('Mm1OoIiFfEe'))
 @pytest.mark.parametrize('uplo, m, n', [('U', 5, 10), ('U', 10, 10),
                                         ('L', 10, 5), ('L', 10, 10)])


Code review comments:
diff_hunk: @@ -1410,3 +1410,40 @@ subroutine <prefix2c>heequb(lower,n,a,lda,s,scond,amax,work,info)
     integer intent(out) :: info
 
 end subroutine <prefix2c>heequb
+
+
+subroutine <prefix>sytri(lower,n,a,lda,ipiv,work,info)
+    ! ?SYTRI computes the inverse of a symmetric indefinite matrix
+    ! A using the factorization A = U*D*U**T or A = L*D*L**T computed by
+    ! ?SYTRF.
+    threadsafe
+    callstatement (*f2py_func)((lower?"L":"U"),&n,a,&lda,ipiv,work,&info)
+    callprotoargument char*, F_INT*, <ctype>*, F_INT*, F_INT*, <ctype>*, F_INT*
+
+    integer optional intent(in), check(lower==0||lower==1) :: lower = 0
+    integer intent(hide), depend(a) :: n=shape(a,1)
+    <ftype> dimension(n,n),intent(in,out,copy,out=inv_a) :: a
+    integer depend(n),intent(hide) :: lda = max(n,1)
path: scipy/linalg/flapack_sym_herm.pyf.src
text: @user1:
In LAPACK function conventions, lda is the number of rows and n is the number of columns. So typically they are not dependent. This is to allow for passing 20x10 array and be able to solve 5x5 corner of it where n=5 defines the size of the problem and lda tells you how many items to skip for the next column. Hence a size should indeed be `(lda, n)` and `lda = max(1, shape(a, 0))`.

@author:
Thanks for the explanation! Have changed it