linalg: solve

[perazz]

Linear system solution $Ax=b$, for both 1-rhs (b(:)) or n-rhs (b(:,:)) cases, based on LAPACK *GESV functions.

• base implementation
• tests
• exclude unsupported xdp
• documentation
• submodule
• examples
• pure subroutine interface

Prior art:

• NumPy: solve(a, b)
• Scipy: solve(a, b, lower=False, overwrite_a=False, overwrite_b=False, check_finite=True, assume_a='gen', transposed=False)
• IMSL: lu_solve(a, b, transpose=False)
• Add solve based on OpenBLAS #710 cc @zoziha

Proposed implementation:

• pure (only intent(in) arguments)
• "expert" interface (intent(out) error handler, a can be overwritten)

Example calls:

• x = solve(a,b) -> simplest pure call
• x = solve(a,b,overwrite_a=.false.,err) -> expert call:
  • logical(lk), optional, intent(in) :: overwrite_a = option to avoid internal allocation, but destroy input matrix a. Default: .false.
  • type(linalg_state_type), intent(out) :: err = return state variable

cc: @jvdp1 @zoziha @jalvesz @fortran-lang/stdlib: I believe this is ready for consideration.

[zoziha]

LGTM. @perazz, thanks for your efforts, stdlib gradually has a user-friendly, high-level BLAS API.

[jvdp1]

Thank you. Here are a few suggestions/comments.

[perazz]

@jvdp1 @zoziha I've implemented a pure subroutine interface, that looks like this:

call solve_lu(a,b,x,pivot=ipiv,overwrite_a=.true.,err=err)

• Is solve_lu a good name? I've chosen this because it follows the same pattern as solve_lstsq in the least-squares interface
• Both the subroutine and the function interfaces are now pure

Please, let me know what you think of this approach, I will complete the documentation for this after we've defined the final version. Thanks!

[jvdp1]

@jvdp1 @zoziha I've implemented a pure subroutine interface

Great! Thanks @perazz. LGTM

• Is solve_lu a good name? I've chosen this because it follows the same pattern as solve_lstsq in the least-squares interface

LGTM. it makes sense and gives some info on how it is solves.

[jvdp1]

Thank you @perazz . LGTM, and it is almost ready to be merged. Nice addition!

[perazz]

Dear all @jvdp1 @jalvesz @zoziha I've documented the subroutine interface and added a no-allocation example.
I believe this should be ready now. pivot is the only pre-allocated array needed in this case. I suggest we wait until we have the full picture (what pre-allocated arrays do svd, eig, etc. need? ) so we can create a linalg_storage type a posteriori

[jvdp1]

Thank you @perazz . I will merge it.
As discussed here, I will prepare a new release today/tomorrow.

[perazz]

Thanks a lot @jvdp1