Solving a linear system with a symmetric positive definite matrix

[aamrindersingh]

Resolves #1069

This PR adds two interfaces to stdlib_linalg for solving linear systems with symmetric positive definite (SPD) matrices. The solve_chol interface is for two step workflows where you factorize once and solve multiple times , common in iterative algorithms like quadratic programming. Usage: call cholesky(A, L, lower=.true.) then call solve_chol(L, b, x, lower=.true.). It wraps LAPACK's POTRS. The cholesky_solve interface is a one shot convenience that factorizes and solves in a single call using LAPACK's POSV: call cholesky_solve(A, b, x). This is simpler for one time solves where convenience matters more than reusing factorizations.

The key design decision is making lower required in solve_chol to prevent silent failures when it doesn't match the factorization step. If someone calls cholesky(A, L, lower=.false.) but forgets to specify lower=.false. in solve_chol, they would get wrong results silently. Making it required forces explicit matching. In contrast, cholesky_solve keeps lower optional (defaults to .true.) since factorization happens internally. Both interfaces support real(sp/dp/qp/xdp) and complex(sp/dp/qp/xdp) when available, following the existing solve_lu naming pattern for consistency.

Testing includes basic SPD solves, multiple right-hand sides, upper/lower triangular forms, real and complex types, error handling, residual verification (||A*x - b|| / ||b|| < sqrt(epsilon)), and comparison with solve_lu. All tests pass with GFortran 13.2 and Intel ifx 2024.0.

[codecov]

Codecov Report

❌ Patch coverage is 0% with 42 lines in your changes missing coverage. Please review.
✅ Project coverage is 68.31%. Comparing base (0ede301) to head (f429a98).
⚠️ Report is 38 commits behind head on master.

Files with missing lines	Patch %	Lines
example/linalg/example_solve_lower_chol.f90	0.00%	16 Missing ⚠️
example/linalg/example_solve_upper_chol.f90	0.00%	16 Missing ⚠️
example/linalg/example_solve_chol.f90	0.00%	10 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1093      +/-   ##
==========================================
- Coverage   68.55%   68.31%   -0.25%     
==========================================
  Files         396      399       +3     
  Lines       12746    12788      +42     
  Branches     1376     1383       +7     
==========================================
- Hits         8738     8736       -2     
- Misses       4008     4052      +44     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[loiseaujc]

Thanks for your contribution. It looks pretty good. A couple of general things:

• Whenever you make a call to a lapack function returning an info flag, it is best to have this call followed by call handle_whateverlapackfunction_info(this, ..., err0) to parse the returned info flag. These handling functions are available (or need to be implemented) from the stdlib_linalg_lapack_aux module.

• Could you add the specs of the different functions to the docs ?

• By default, the function cholesky_solve overwrites A with its Cholesky factorization. While this is beneficial for the sake of performances (less allocation or data movement), users may want to keep A untouched. Could you extend the interface and the internal logic with a optional boolean flag overwrite_a (default .false.) similar to what is being done with solve_lu ?

Moreover, I see that you have created a new submodule. It's fine, but since we already have an stdlib_linalg_solve submodule, and solve_cholesky is just a variation around this theme, I think it would be better to have your implementations there. What do you think @perazz ?

[aamrindersingh]

Thanks @loiseaujc

I have addressed all the points:-
Added handle_potrs_info and handle_posv_info to stdlib_linalg_lapack_aux for proper LAPACK info flag handling
Added overwrite_a flag to cholesky_solve (default .false. to preserve A, matching solve_lu pattern)
Moved implementations into stdlib_linalg_solve.fypp and deleted the separate submodule
Added specs to doc/specs/stdlib_linalg.md
Updated tests to cover the new overwrite_a behavior

Let me know if anything else needs adjustment.

[aamrindersingh]

Solved all merge conficts!
ready for review
Thank you

[Copilot]

Pull request overview

This PR adds two interfaces to stdlib_linalg for solving linear systems with symmetric positive definite (SPD) matrices, addressing issue #1069. The implementation wraps LAPACK's POTRS and POSV routines to provide both a two-step workflow (cholesky + solve_chol) for repeated solves and a one-shot convenience interface (cholesky_solve) for single solves.

Changes:

• Added solve_chol interface for solving systems using pre-computed Cholesky factors (POTRS wrapper)
• Added cholesky_solve interface for one-shot factorization and solve (POSV wrapper)
• Implemented error handlers handle_potrs_info and handle_posv_info for LAPACK return codes

Reviewed changes

Copilot reviewed 9 out of 9 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
src/linalg/stdlib_linalg_solve.fypp	Implementation of solve_chol and cholesky_solve routines following existing solve_lu patterns
src/linalg/stdlib_linalg.fypp	Public interface declarations with comprehensive documentation comments
src/lapack/stdlib_linalg_lapack_aux.fypp	LAPACK error handlers for POTRS and POSV routines
test/linalg/test_linalg_solve_chol.fypp	Test suite covering basic solves, multiple RHS, overwrite behavior, and all supported types
test/linalg/CMakeLists.txt	Integration of new test suite into build system
example/linalg/example_solve_chol.f90	Example demonstrating two-step workflow with pre-computed factors
example/linalg/example_cholesky_solve.f90	Example demonstrating one-shot solve interface
example/linalg/CMakeLists.txt	Registration of example programs
doc/specs/stdlib_linalg.md	Comprehensive documentation for both interfaces with syntax, arguments, and examples

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Pull request overview

Copilot reviewed 9 out of 9 changed files in this pull request and generated no new comments.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[aamrindersingh]

@jvdp1 Addressed all Copilot suggestions
CI Build Failure is unrelated.
Thanks

[loiseaujc]

It looks good to me and seems pretty much ready for merging. Could you add two additional tests though:

• A symmetric but indefinite
• A symmetric positive semi-definite

just to make sure the handle_potrs_info and handle_posv_info capture these cases correctly ?

[aamrindersingh]

@loiseaujc
Addressed all review comments and also added tests for indefinite and semi definite matrices to verify error handling.
CI should also pass now.
Thanks

[loiseaujc]

Code-wise, it looks really good. And very thorough unit-testing. I'll go through the specs tomorrow but I think you're very close.

[perazz]

I stand by @loiseaujc's comments as far as the quality of the PR.

I would like to have a discussion on the choice of the API names: I think solve_chol and cholesky_solve are interchangeable, so not very clear from the user's standpoint. Given we have solve, solve_lu, solve_lstsq, solve_constrained_lstsq, ... a solve_cholesky or solve_chol interface looks more in line with the current API to me. in other words, I think cholesky_solve would be confusing because it does not "imply" any difference from solve_chol. Here is a proposal to unify the API:

• solve_chol or solve_cholesky -> replace cholesky_solve i.e. start from matrix "A"; retain the lower argument
• solve_upper_cholesky or solve_upper_chol -> solve with (u, b, x) pre-formed (drop lower argument)
• solve_lower_cholesky or solve_lower_chol -> solve with (l, b, x) pre-formed (drop lower argument)

[aamrindersingh]

Thanks @perazz,
Going with the shorter _chol suffix for consistency. Almost done with the API restructuring will push in a moment after revalidating everything.

[loiseaujc]

The subroutines solve_lower_chol and solve_upper_chol are essentially duplicated code with a single change. Maintenance wise, a better approach (I guess) would be to have a low-level driver solve_factorized_cholesky(a, b, x, uplo, err) with all arguments required, and then have solve_lower_chol and solve_upper_chol simply call this low-level driver with the appropriate uplo input. Note that the low-level driver need not be exported at the module level.

[perazz]

Thank you @aamrindersingh, it's in a pretty good state. I added a couple clarifying edits (please double check no leftovers)

[loiseaujc]

As far as I can see, I think it looks pretty much ready.

[loiseaujc]

Great job @aamrindersingh ! As far as I'm concerned, this PR is ready for merging. I'll give @jvdp1 and @jalvesz a few days if they want to check out the code and have any comments. If they green light it, I guess we'll merge by the end of next week.

[perazz]

With no more outstanding comments, I will go ahead and merge this one.