linalg: Matrix Inverse

[perazz]

Compute the multiplicative inverse of a real or complex square matrix: $A \cdot A^{-1} = A^{-1} \cdot A = I_n$ . Based on LAPACK General factorization (*GETRF) and inversion (*GETRI).

• [x] base implementation
• [x] tests
• [x] exclude unsupported xdp
• [x] documentation
• [x] submodule
• [x] examples
• [x] subroutine interface

Prior art

• Numpy: B = inv(A)
• Scipy: inv(A, overwrite_a=False, check_finite=True)
• IMSL: .i.A

Proposed implementation

• B = inv(A [, err]) function interface
• call invert(A [, pivot] [, err]) in-place (no-allocation) subroutine interface (optional pivot array)
• .inv.A operator interface

cc: @jalvesz @jvdp1 @loiseaujc @fortran-lang/stdlib

[loiseaujc]

As far as this is my first review for stdlib, it looks pretty good to me. I was jut wondering about the tests. Wouldn't including a test with a known singular matrix be useful to make sure the whole error handling is working correctly?

[perazz]

Great idea @loiseaujc, thanks!

• [x] test for singular matrix returning LINALG_ERROR
• [x] noticed that complex tests were not running -> add them

e10e4c4

[perazz]

From Fortran Monthly call:

• [x] in subroutine invert, add interface to not operate in-place but return inverse to another variable (i.e., call invert(a, am1, [, pivot] [, err])
• [x] in .inv.A operator, with singular matrix, do not return all zeros, but stop the program

[perazz]

Thanks a lot @loiseaujc @jvdp1 @jalvesz for the reviews - I believe I've addressed all your comments (see linked commits in each comment). In d0af9be I've introduced a tiny behavior change:

On exceptions raised when the operator interface (.inv.A) is used, instead of error stopping the program, I think it is more reasonable to instead return an array of NaNs (which are initialized from ieee_value(0.0,ieee_quiet_nan). The reason I think this is better is:

• operators can be chained, I don't think we want the program to stop in between. For that, there are the subroutine and function interfaces
• other Fortran operators (+, -, but also matmul etc. ) would behave exactly in the same way: when something goes wrong, they typically return NaNs.

So please do let me know if you also like this update.

[jalvesz]

Given your idea to enable the operator interface to run through the computation without stopping in case of singular matrices I just suggested a small addition to the documentation. This change seems reasonable to me as even a scalar division would not stop the program but also simply return a NaN. Otherwise LGTM!!

[perazz]

Ok @jvdp1 @jalvesz I've restored the xdp precision. Please let me know if you have further comments. This is probably ready to be merged sometime next week.

[perazz]

Thank you all, I will merge this now.