Solving a linear system with a symmetric positive definite matrix

[loiseaujc]

Motivation

I'm currently doing some more modernization work on the QuadProg solver for solving strictly convex quadratic programs. One of the critical steps is solving a linear system P x = q where P is a symmetric positive definite matrix. The matrix is being factorized using its Cholesky decomposition (originally using dpotrf from lapack, now replaced with cholesky from stdlib). The linear solve currently makes an explicit call to dpotrs from lapack but, for the sake of code clarity, I'd like to replace it with an equivalent chol_solve which does not currently exist in stdlib.

As of today, stdlib has a solve_lu subroutine whose interface reads

  interface solve_lu        
     #:for nd,ndsuf,nde in ALL_RHS
     #:for rk,rt,ri in RC_KINDS_TYPES
     pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err)     
         !> Input matrix a[n,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$
         !> Result array/matrix x[n] or x[n,nrhs]     
         ${rt}$, intent(inout), contiguous, target :: x${nd}$
         !> [optional] Storage array for the diagonal pivot indices
         integer(ilp), optional, intent(inout), target :: pivot(:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
     end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$
     #:endfor
     #:endfor    
  end interface solve_lu     

We could write an equivalent routine for chol_solve (or solve_chol if we want to be consistent with the LU call). Almost all of the necessary routines are already present in stdlib_lapack_solve_chol.fypp and stdlib_lapack_solve_chol_comp.fypp. Solving a linear system with a symmetric positive definite matrix could then be done using two successives calls, e.g. (omitting optional arguments)

!> Factorize the matrix (in-place).
call cholesky(P)
!> Solve the linear system.
call chol_solve(P, q, x)

Prior Art

• scipy.linalg has the function cho_solve whose interface reads cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True) (see here).
• In Julia, you would do F = cholesky(P) followed by x = F \ b where F is the equivalent of a Fortran derived-type specialized for the Cholesky factorization (see Derived type for matrix factorizations #1040 for a discussion with the QR example).

Additional Information

ping: @perazz, @jvdp1, @jalvesz