linalg: Matrix Inverse #828

perazz · 2024-05-23T14:02:37Z

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).

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 · 2024-05-27T08:17:30Z

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 · 2024-05-27T10:38:50Z

Great idea @loiseaujc, thanks!

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

e10e4c4

perazz · 2024-06-04T17:07:58Z

From Fortran Monthly call:

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

jalvesz · 2024-06-08T14:17:06Z

example/linalg/example_inverse3.f90

+  Am1 = A
+
+  ! Invert matrix
+  call invert(Am1)


Why not use this example to show how to use the subroutine interface with the optional output matrix:

Suggested change

Am1 = A

! Invert matrix

call invert(Am1)

! Invert matrix without modifying the original matrix

Am1 = A

call invert(Am1)

! Which would be equivalent to: call invert(A,Am1)

jalvesz · 2024-06-08T14:19:03Z

example/linalg/example_inverse1.f90

@@ -0,0 +1,22 @@
+! Matrix inversion example: operator interface
+program example_inverse1


I would propose to give explicit names to the example files such as:
example_inverse_operator, example_inverse_function, example_inverse_subroutine

jvdp1

Thank you @perazz . On overall LGTM. Here are some suggestions.

jvdp1 · 2024-06-09T17:38:42Z

doc/specs/stdlib_linalg.md

@@ -777,7 +777,8 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \(

 `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument.

-`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `minval(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` argument.
+`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `m


Is this a typo?

jvdp1 · 2024-06-09T17:39:32Z

doc/specs/stdlib_linalg.md

+This operator returns the inverse of a `real` or `complex` square matrix \( A \).
+The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
+
+This interface is equivalent to the function version of  [[stdlib_linalg(module):inv(interface)]]. 


Suggested change

This interface is equivalent to the function version of [[stdlib_linalg(module):inv(interface)]].

This interface is equivalent to the function [[stdlib_linalg(module):inv(interface)]].

jvdp1 · 2024-06-09T17:40:28Z

doc/specs/stdlib_linalg.md

+{!example/linalg/example_inverse1.f90!}
+```
+
+## `invert` - Inversion of a square matrix. 


Suggested change

## `invert` - Inversion of a square matrix.

## `invert` - Inversion of a square matrix

jvdp1 · 2024-06-09T17:42:57Z

doc/specs/stdlib_linalg.md

+`inva` (optional): Shall be a rank-2, square, `real` or `complex` array with the same size, and kind as `a`. 
+On output, it contains the inverse of `a`.
+
+`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices. 


Suggested change

`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices.

`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, that contains the diagonal pivot indices on return. It is an `intent(inout)` argument.

jvdp1 · 2024-06-09T17:44:23Z

doc/specs/stdlib_linalg.md

+
+### Return value
+
+Replaces matrix \( A \) with its inverse, \(A^{-1}\).


That is not always the case

Suggested change

Replaces matrix $ A $ with its inverse, $A^{-1}$.

Computes the inverse of the matrix $ A $, $A^{-1}\, and returns it either in \( A $ or in another matrix.

jvdp1 · 2024-06-09T17:45:15Z

doc/specs/stdlib_linalg.md

+
+`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument.
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.


Suggested change

`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.

`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.

jvdp1 · 2024-06-09T17:50:56Z

src/stdlib_linalg_inverse.fypp

+
+         ! Has a pre-allocated pivots storage array been provided? 
+         if (present(pivot)) then 
+            ipiv => pivot


Maybe to allow larger vectors, could something like this be useful:

Suggested change

ipiv => pivot

ipiv => pivot(1:n)

jvdp1 · 2024-06-09T18:11:14Z

test/linalg/test_linalg_inverse.fypp

+
+        #:for rk,rt,ri in RC_KINDS_TYPES
+        #:if rk!="xdp"
+        tests = [tests,new_unittest("${ri}$_eye_inverse",test_${ri}$_eye_inverse), &


I propose to add a test for random SPD matrices. In my own implementation I have something like that:

do i = 2, ndim orimat = 0 call random_number(mat(1:i, 1:i)) orimat(1:i, 1:i) = 0.5_real64 * matmul(mat(1:i, 1:i), transpose(mat(1:i, 1:))) mat = orimat call inv(mat) !identity matrix identity(1:i, 1:i) = reshape([(merge(1._real64, 0._real64, j/i.eq.mod(j,i)), j = 0, i**2 - 1)], [i, i]) mat(1:i, 1:i) = matmul(mat(1:i, 1:i), orimat(1:i, 1:i)) call check(error , all(abs(mat(1:i, 1:i) - identity(1:i, 1:i)) < tol_real64) & , 'inv_pd: wrong inverse for '//value2char(i)) if(allocated(error))return enddo

jvdp1 · 2024-06-09T18:20:59Z

src/stdlib_linalg_inverse.fypp

+
+            ! Get optimal worksize (returned in work(1)) (inflate by a 5% safety margin)
+            nb = stdlib_ilaenv(1,'${ri}$getri',' ',n,-1,-1,-1)
+            lwork = ceiling(1.05*n*nb,kind=ilp)


We should be careful with integer overflow in this case (I have been bitten several times ;). It might be good to have a check for lwork

perazz added 10 commits May 23, 2024 15:22

add source file

f39bf52

working implementation

c0da712

reorganize as submodule

92b1ac5

add tests

e9a3f5b

option for pre-allocated pivot indices

b233156

add examples

dc0127b

document .inv.A

75a5a23

document invert

bf0dfd3

document inv

db90b37

typo

6a1c397

perazz marked this pull request as ready for review May 23, 2024 15:50

test for singular matrix; activate complex matrix tests

e10e4c4

Merge branch 'master' into matrix_inverse

db714bb

perazz added 4 commits June 5, 2024 09:45

subroutine version, non in-place

5f320f9

test subroutine version, non in-place

47ee61a

.inv. operator: error stop on issues

e2159bc

Merge branch 'master' into matrix_inverse

b23c811

jalvesz suggested changes Jun 8, 2024

View reviewed changes

jvdp1 reviewed Jun 9, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

linalg: Matrix Inverse #828

linalg: Matrix Inverse #828

perazz commented May 23, 2024 •

edited

loiseaujc commented May 27, 2024

perazz commented May 27, 2024 •

edited

perazz commented Jun 4, 2024 •

edited

jalvesz Jun 8, 2024

jalvesz Jun 8, 2024

jvdp1 left a comment

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

jvdp1 Jun 9, 2024

		@@ -0,0 +1,22 @@
		! Matrix inversion example: operator interface
		program example_inverse1

	This interface is equivalent to the function version of [[stdlib_linalg(module):inv(interface)]].
	This interface is equivalent to the function [[stdlib_linalg(module):inv(interface)]].

	## `invert` - Inversion of a square matrix.
	## `invert` - Inversion of a square matrix

	`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices.
	`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, that contains the diagonal pivot indices on return. It is an `intent(inout)` argument.


		### Return value

		Replaces matrix \( A \) with its inverse, \(A^{-1}\).

	Replaces matrix \( A \) with its inverse, \(A^{-1}\).
	Computes the inverse of the matrix \( A \), \(A^{-1}\, and returns it either in \( A \) or in another matrix.


		`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument.

		`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.

linalg: Matrix Inverse #828

Are you sure you want to change the base?

linalg: Matrix Inverse #828

Conversation

perazz commented May 23, 2024 • edited

Prior art

Proposed implementation

loiseaujc commented May 27, 2024

perazz commented May 27, 2024 • edited

perazz commented Jun 4, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jvdp1 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

perazz commented May 23, 2024 •

edited

perazz commented May 27, 2024 •

edited

perazz commented Jun 4, 2024 •

edited