Repository: scipy/scipy
Pull request: ENH:linalg: Add LAPACK wrappers pptrf/pptrs/ppsv/pptri/ppcon
As title mentions.

Closes #12431

Problem statement:
Request for Python wrapper for LAPACK's ?pptrf (Cholesky factorization in packed format)
**Is your feature request related to a problem? Please describe.**

scipy.linalg.lapack already has wrappers for ?tpsv (triangular solve in packed format).  ?pptrf (Cholesky factorization in packed format) seems to be a common enough function that can be used to generate the matrix passed to ?tpsv, but it is not wrapped at the Python level.

**Describe the solution you'd like**

Add a wrapper for ?pptrf.

**Describe alternatives you've considered**

Manually converting between unpacked and packed formats, but this is not exactly fast.

Hints from issue discussion:
@anntzer consider adding items for all wrappers you want?

- [ ] `?pptrf`
- [ ] `?ZZZZ`
etc. in that way one can have more in one.
(Right now that's the only one I need, but I will add more as needed :))

Commit message: ENH:linalg: Add wrappers pptrf/pptrs/ppsv

Patch to review:
diff --git a/scipy/linalg/flapack_other.pyf.src b/scipy/linalg/flapack_other.pyf.src
index 39c3ac7a2365..c1567b6c49c1 100644
--- a/scipy/linalg/flapack_other.pyf.src
+++ b/scipy/linalg/flapack_other.pyf.src
@@ -1170,6 +1170,71 @@ subroutine <prefix>gglse_lwork(m,n,p,a,lda,b,ldb,c,d,x,work,lwork,info)
 
 end subroutine <prefix>gglse_lwork
 
+subroutine <prefix>pptrf(lower,n,ap,info,L)
+    ! ?PPTRF computes the Cholesky factorization of a symmetric/hermitian
+    ! positive definite matrix A stored in packed format.
+    ! 
+    ! The factorization has the form
+    !    A = U**T * U,  if UPLO = 'U', or
+    !    A = L  * L**T,  if UPLO = 'L',
+    ! where U is an upper triangular matrix and L is lower triangular.
+
+    callstatement (*f2py_func)((lower?"L":"U"),&n,ap,&info)
+    callprotoargument char*,int*,<ctype>*,int*
+    
+    integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
+    integer intent(int),check(n>=0) :: n
+    <ftype> dimension(L),depend(n),intent(in,out,copy,out=ul) :: ap
+    integer intent(hide),depend(ap,n),check(L>=(n*(n+1)/2)) :: L = len(ap)
+    integer intent(out) :: info
+
+end subroutine<prefix>pptrf
+
+subroutine <prefix>pptrs(lower,n,nrhs,ap,b,ldb,info,L)
+    ! DPPTRS solves a system of linear equations A*X = B with a symmetric
+    ! positive definite matrix A in packed storage using the Cholesky
+    ! factorization A = U**T*U or A = L*L**T computed by DPPTRF.
+
+    callstatement (*f2py_func)((lower?"L":"U"),&n,&nrhs,ap,b,&ldb,&info)
+    callprotoargument char*,int*,int*,<ctype>*,<ctype>*,int*,int*
+    
+    integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
+    integer intent(int),check(n>=0) :: n
+    integer intent(hide),depend(b) :: ldb = max(1, shape(b,0))
+    integer intent(hide),depend(b) :: nrhs = shape(b,1)
+    <ftype> dimension(L),depend(n),intent(in) :: ap
+    integer intent(hide),depend(ap,n),check(L>=(n*(n+1)/2)) :: L = len(ap)
+    <ftype> dimension(ldb,nrhs),intent(in,out,copy,out=x) :: b
+    integer intent(out) :: info
+
+end subroutine <prefix>pptrs
+
+subroutine <prefix>ppsv(lower,n,nrhs,ap,b,ldb,info,L)
+    ! DPPSV computes the solution to a real system of linear equations
+    !    A * X = B,
+    ! where A is an N-by-N symmetric positive definite matrix stored in
+    ! packed format and X and B are N-by-NRHS matrices.
+    !
+    ! The Cholesky decomposition is used to factor A as
+    !    A = U**T* U,  if UPLO = 'U', or
+    !    A = L * L**T,  if UPLO = 'L',
+    ! where U is an upper triangular matrix and L is a lower triangular
+    ! matrix.  The factored form of A is then used to solve the system of
+    ! equations A * X = B.
+    callstatement (*f2py_func)((lower?"L":"U"),&n,&nrhs,ap,b,&ldb,&info)
+    callprotoargument char*,int*,int*,<ctype>*,<ctype>*,int*,int*
+    
+    integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
+    integer intent(int),check(n>=0) :: n
+    integer intent(hide),depend(b) :: ldb = max(1, shape(b,0))
+    integer intent(hide),depend(b) :: nrhs = shape(b,1)
+    <ftype> dimension(L),depend(n),intent(in) :: ap
+    integer intent(hide),depend(ap,n),check(L>=(n*(n+1)/2)) :: L = len(ap)
+    <ftype> dimension(ldb,nrhs),intent(in,out,copy,out=x) :: b
+    integer intent(out) :: info
+
+end subroutine <prefix>ppsv
+
 subroutine <prefix2>sbev(ab,compute_v,lower,n,ldab,kd,w,z,ldz,work,info)
     ! in :Band:dsbev.f
     ! principally <s,d>sbevd does the same, and are recommended for use.
diff --git a/scipy/linalg/lapack.py b/scipy/linalg/lapack.py
index ecaedee393c1..0cf656a9572b 100644
--- a/scipy/linalg/lapack.py
+++ b/scipy/linalg/lapack.py
@@ -513,6 +513,21 @@
    cpteqr
    zpteqr
 
+   spptrf
+   dpptrf
+   cpptrf
+   zpptrf
+
+   spptrs
+   dpptrs
+   cpptrs
+   zpptrs
+
+   sppsv
+   dppsv
+   cppsv
+   zppsv
+
    crot
    zrot
 
diff --git a/scipy/linalg/tests/test_lapack.py b/scipy/linalg/tests/test_lapack.py
index 3e1db62c8b3c..78b2f41796eb 100644
--- a/scipy/linalg/tests/test_lapack.py
+++ b/scipy/linalg/tests/test_lapack.py
@@ -19,9 +19,9 @@
 
 from numpy.random import rand, randint, seed
 
-from scipy.linalg import _flapack as flapack, lapack
-from scipy.linalg import inv, svd, cholesky, solve, ldl, norm, block_diag, qr
-from scipy.linalg import eigh
+from scipy.linalg import (_flapack as flapack, lapack, inv, svd, cholesky,
+                          solve, ldl, norm, block_diag, qr, eigh)
+
 from scipy.linalg.lapack import _compute_lwork
 from scipy.stats import ortho_group, unitary_group
 
@@ -1748,7 +1748,8 @@ def test_syequb():
         assert_equal(np.log2(s).astype(int), desired_log2s)
 
 
-@pytest.mark.skipif(True, reason="Failing on some OpenBLAS version, see gh-12276")
+@pytest.mark.skipif(True,
+                    reason="Failing on some OpenBLAS version, see gh-12276")
 def test_heequb():
     # zheequb has a bug for versions =< LAPACK 3.9.0
     # See Reference-LAPACK gh-61 and gh-408
@@ -2923,3 +2924,29 @@ def test_ptsvx_NAG(d, e, b, x):
     df, ef, x_ptsvx, rcond, ferr, berr, info = ptsvx(d, e, b)
     # determine ptsvx's solution and x are the same.
     assert_array_almost_equal(x, x_ptsvx)
+
+
+@pytest.mark.parametrize('dtype', DTYPES)
+def test_pptrs_pptrf_ppsv(dtype):
+    seed(1234)
+    # Manual conversion to/from packed format is feasible here.
+    n, nrhs = 10, 4
+    a = generate_random_dtype_array([n, n], dtype=dtype)
+    b = generate_random_dtype_array([n, nrhs], dtype=dtype)
+    a = a.conj().T + a + np.eye(n, dtype=dtype) * dtype(5.)
+    inds = ([x for y in range(1, n+1) for x in range(y)],
+            [y-1 for y in range(1, n+1) for x in range(y)])
+    ap = a[inds]
+    ppsv, pptrf, pptrs = get_lapack_funcs(('ppsv', 'pptrf', 'pptrs'),
+                                          dtype=dtype)
+
+    ul, info = pptrf(n, ap)
+    aul = cholesky(a)[inds]
+    assert_allclose(ul, aul, rtol=np.finfo(dtype).eps*50)
+
+    x, info = pptrs(n, ul, b)
+    bx = solve(a, b)
+    assert_allclose(x, bx, rtol=np.finfo(dtype).eps*50)
+
+    xv, info = ppsv(n, ap, b)
+    assert_allclose(xv, bx, rtol=np.finfo(dtype).eps*50)


Code review comments:
diff_hunk: @@ -2923,3 +2924,29 @@ def test_ptsvx_NAG(d, e, b, x):
     df, ef, x_ptsvx, rcond, ferr, berr, info = ptsvx(d, e, b)
     # determine ptsvx's solution and x are the same.
     assert_array_almost_equal(x, x_ptsvx)
+
+
+@pytest.mark.parametrize('dtype', DTYPES)
path: scipy/linalg/tests/test_lapack.py
text: @user1:
Should test UPLO?

@author:
Bah, I forgot lower keyword. Will add.

diff_hunk: @@ -2923,3 +2924,29 @@ def test_ptsvx_NAG(d, e, b, x):
     df, ef, x_ptsvx, rcond, ferr, berr, info = ptsvx(d, e, b)
     # determine ptsvx's solution and x are the same.
     assert_array_almost_equal(x, x_ptsvx)
+
+
+@pytest.mark.parametrize('dtype', DTYPES)
+def test_pptrs_pptrf_ppsv(dtype):
+    seed(1234)
+    # Manual conversion to/from packed format is feasible here.
+    n, nrhs = 10, 4
+    a = generate_random_dtype_array([n, n], dtype=dtype)
+    b = generate_random_dtype_array([n, nrhs], dtype=dtype)
+    a = a.conj().T + a + np.eye(n, dtype=dtype) * dtype(5.)
+    inds = ([x for y in range(1, n+1) for x in range(y)],
+            [y-1 for y in range(1, n+1) for x in range(y)])
+    ap = a[inds]
+    ppsv, pptrf, pptrs = get_lapack_funcs(('ppsv', 'pptrf', 'pptrs'),
+                                          dtype=dtype)
+
+    ul, info = pptrf(n, ap)
+    aul = cholesky(a)[inds]
+    assert_allclose(ul, aul, rtol=np.finfo(dtype).eps*50)
+
+    x, info = pptrs(n, ul, b)
path: scipy/linalg/tests/test_lapack.py
text: @user1:
Shouldn't the routine be able to figure out `n`?

@author:
See #12333 for a similar discussion

@user1:
By the time I got there, I forgot it was in a packed format.

diff_hunk: @@ -2923,3 +2924,50 @@ def test_ptsvx_NAG(d, e, b, x):
     df, ef, x_ptsvx, rcond, ferr, berr, info = ptsvx(d, e, b)
     # determine ptsvx's solution and x are the same.
     assert_array_almost_equal(x, x_ptsvx)
+
+
+@pytest.mark.parametrize('lower', [False, True])
+@pytest.mark.parametrize('dtype', DTYPES)
+def test_pptrs_pptri_pptrf_ppsv_ppcon(dtype, lower):
+    seed(1234)
+    atol = np.finfo(dtype).eps*100
+    # Manual conversion to/from packed format is feasible here.
+    n, nrhs = 10, 4
+    a = generate_random_dtype_array([n, n], dtype=dtype)
+    b = generate_random_dtype_array([n, nrhs], dtype=dtype)
+
+    a = a.conj().T + a + np.eye(n, dtype=dtype) * dtype(5.)
+    if lower:
+        inds = ([x for y in range(n) for x in range(y, n)],
+                [y for y in range(n) for x in range(y, n)])
+    else:
+        inds = ([x for y in range(1, n+1) for x in range(y)],
+                [y-1 for y in range(1, n+1) for x in range(y)])
+    ap = a[inds]
+    ppsv, pptrf, pptrs, pptri, ppcon = get_lapack_funcs(('ppsv', 'pptrf',
+                                                         'pptrs', 'pptri',
+                                                         'ppcon'), dtype=dtype)
path: scipy/linalg/tests/test_lapack.py
text: @user1:
If you add `ilp64="preferred"` here, also that path gets tested on the CI. (Or, parametrize the test; `scipy.linalg.blas.HAS_ILP64` is the relevant enable flag.)

@author:
Ah sure. This will take some time for me to get used to this.