linalg: Matrix Inverse

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
+al(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` argument.
 
 `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `A` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
 
@@ -879,15 +880,15 @@ This interface is equivalent to the `pure` version of determinant [[stdlib_linal
 
 ### Syntax
 
-`c = ` [[stdlib_linalg(module):operator(.det.)(interface)]] `(a)`
+`c = ` [[stdlib_linalg(module):operator(.det.)(interface)]] `a`
 
 ### Arguments
 
 `a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument.
 
 ### Return value
 
-Returns a real scalar value that represents the determinnt of the matrix.
+Returns a real scalar value that represents the determinant of the matrix.
 
 Raises `LINALG_ERROR` if the matrix is singular.
 Raises `LINALG_VALUE_ERROR` if the matrix is non-square.
@@ -989,4 +990,116 @@ Exceptions trigger an `error stop`, unless argument `err` is present.
 ```fortran
 {!example/linalg/example_svdvals.f90!}
 ```
+## `.inv.` - Inverse operator of a square matrix
 
+### Status
+
+Experimental
+
+### Description
+
+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)]]. 
+
+### Syntax
+
+`b = ` [[stdlib_linalg(module):operator(.inv.)(interface)]] `a`
+
+### Arguments
+
+`a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument.
+
+### Return value
+
+Returns a rank-2 square array with the same type, kind and rank as `a`, that contains the inverse of `a`.
+
+Exceptions always trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_inverse1.f90!}
+```
+
+## `invert` - Inversion of a square matrix. 
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine inverts a square `real` or `complex` matrix in-place.
+The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
+
+On return, the input matrix `a` is replaced by its inverse.
+The solver is based on LAPACK's `*GETRF` and `*GETRI` backends.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):invert(interface)]] `(a, [,inva] [, pivot] [, err])`
+
+### Arguments
+
+`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. 
+If `inva` is provided, it is an `intent(in)` argument.
+If `inva` is not provided, it is an `intent(inout)` argument: on output, it is replaced by the inverse of `a`. 
+
+`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. 
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
+
+### Return value
+
+Replaces matrix \( A \) with its inverse, \(A^{-1}\).
+
+Raises `LINALG_ERROR` if the matrix is singular or has invalid size.
+Raises `LINALG_VALUE_ERROR` if `inva` and `a` do not have the same size.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_inverse3.f90!}
+```
+
+## `inv` - Inverse of a square matrix. 
+
+### Status
+
+Experimental
+
+### Description
+
+This function returns the inverse of a square `real` or `complex` matrix in-place.
+The inverse, \( A^{-1} \), is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
+
+The solver is based on LAPACK's `*GETRF` and `*GETRI` backends.
+
+### Syntax
+
+`b ` [[stdlib_linalg(module):inv(interface)]] `(a, [, err])`
+
+### Arguments
+
+`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.
+
+### Return value
+
+Returns an array value of the same type, kind and rank as `a`, that contains the inverse matrix \(A^{-1}\).
+
+Raises `LINALG_ERROR` if the matrix is singular or has invalid size.
+If `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_inverse2.f90!}
+```

example/linalg/CMakeLists.txt

@@ -12,6 +12,9 @@ ADD_EXAMPLE(is_skew_symmetric)
 ADD_EXAMPLE(is_square)
 ADD_EXAMPLE(is_symmetric)
 ADD_EXAMPLE(is_triangular)
+ADD_EXAMPLE(inverse1)
+ADD_EXAMPLE(inverse2)
+ADD_EXAMPLE(inverse3)
 ADD_EXAMPLE(outer_product)
 ADD_EXAMPLE(trace)
 ADD_EXAMPLE(state1)

example/linalg/example_inverse1.f90

@@ -0,0 +1,22 @@
+! Matrix inversion example: operator interface
+program example_inverse1
+ use stdlib_linalg_constants, only: dp
+ use stdlib_linalg, only: operator(.inv.),eye
+ implicit none
+
+ real(dp) :: A(2,2), Am1(2,2)
+
+ ! Input matrix (NB Fortran is column major! input columns then transpose)
+ A = transpose(reshape( [4, 3, &
+ 3, 2], [2,2] ))
+
+ ! Invert matrix
+ Am1 = .inv.A
+
+ print *, ' |',Am1(1,:),'|' ! | -2 3 |
+ print *, ' inv(A)= |',Am1(2,:),'|' ! | 3 -4 |
+
+ ! Final check 
+ print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)
+
+end program example_inverse1

example/linalg/example_inverse2.f90

@@ -0,0 +1,22 @@
+! Matrix inversion example: function interface
+program example_inverse2
+ use stdlib_linalg_constants, only: dp
+ use stdlib_linalg, only: inv,eye
+ implicit none
+
+ real(dp) :: A(2,2), Am1(2,2)
+
+ ! Input matrix (NB Fortran is column major! input columns then transpose)
+ A = transpose(reshape( [4, 3, &
+ 3, 2], [2,2] ))
+
+ ! Invert matrix
+ Am1 = inv(A)
+
+ print *, ' |',Am1(1,:),'|' ! | -2 3 |
+ print *, ' inv(A)= |',Am1(2,:),'|' ! | 3 -4 |
+
+ ! Final check 
+ print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)
+
+end program example_inverse2

example/linalg/example_inverse3.f90

@@ -0,0 +1,23 @@
+! Matrix inversion example: in-place inversion
+program example_inverse3
+ use stdlib_linalg_constants, only: dp
+ use stdlib_linalg, only: invert,eye
+ implicit none
+
+ real(dp) :: A(2,2), Am1(2,2)
+
+ ! Input matrix (NB Fortran is column major! input columns then transpose)
+ A = transpose(reshape( [4, 3, &
+ 3, 2], [2,2] ))
+ Am1 = A
+
+ ! Invert matrix
+ call invert(Am1)
+
+ print *, ' |',Am1(1,:),'|' ! | -2 3 |
+ print *, ' inv(A)= |',Am1(2,:),'|' ! | 3 -4 |
+
+ ! Final check 
+ print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)
+
+end program example_inverse3

src/CMakeLists.txt

@@ -29,6 +29,7 @@ set(fppFiles
  stdlib_linalg_cross_product.fypp
  stdlib_linalg_solve.fypp
  stdlib_linalg_determinant.fypp
+ stdlib_linalg_inverse.fypp
  stdlib_linalg_state.fypp
  stdlib_linalg_svd.fypp 
  stdlib_optval.fypp

src/stdlib_linalg.fypp

@@ -20,6 +20,9 @@ module stdlib_linalg
  public :: operator(.det.)
  public :: diag
  public :: eye
+ public :: inv
+ public :: invert
+ public :: operator(.inv.)
  public :: lstsq
  public :: lstsq_space
  public :: solve
@@ -554,6 +557,121 @@ module stdlib_linalg
  #:endfor
  end interface 
 
+ ! Function interface
+ interface inv
+ !! version: experimental 
+ !!
+ !! Inverse of a square matrix
+ !! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-of-a-square-matrix))
+ !!
+ !!### Summary
+ !! This interface provides methods for computing the inverse of a square `real` or `complex` matrix.
+ !! The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
+ !!
+ !!### Description
+ !! 
+ !! This function interface provides methods that return the inverse of a square matrix. 
+ !! Supported data types include `real` and `complex`. 
+ !! The inverse matrix \( A^{-1} \) is returned as a function result. 
+ !! Exceptions are raised in case of singular matrix or invalid size, and trigger an `error stop` if 
+ !! the state flag `err` is not provided. 
+ !!
+ !!@note The provided functions are intended for square matrices.
+ !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
+ !! 
+ #:for rk,rt,ri in RC_KINDS_TYPES
+ #:if rk!="xdp" 
+ module function stdlib_linalg_inverse_${ri}$(a,err) result(inva)
+ !> Input matrix a[n,n]
+ ${rt}$, intent(in) :: a(:,:)
+ !> Output matrix inverse
+ ${rt}$, allocatable :: inva(:,:)
+ !> [optional] state return flag. On error if not requested, the code will stop
+ type(linalg_state_type), optional, intent(out) :: err
+ end function stdlib_linalg_inverse_${ri}$
+ #:endif
+ #:endfor
+ end interface inv
+
+ ! Subroutine interface: in-place factorization
+ interface invert
+ !! version: experimental 
+ !!
+ !! Inversion of a square matrix
+ !! ([Specification](../page/specs/stdlib_linalg.html#invert-inversion-of-a-square-matrix))
+ !!
+ !!### Summary
+ !! This interface provides methods for inverting a square `real` or `complex` matrix in-place.
+ !! The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \).
+ !!
+ !!### Description
+ !! 
+ !! This subroutine interface provides a way to compute the inverse of a matrix. 
+ !! Supported data types include `real` and `complex`. 
+ !! The user may provide a unique matrix argument `a`. In this case, `a` is replaced by the inverse matrix.
+ !! on output. Otherwise, one may provide two separate arguments: an input matrix `a` and an output matrix 
+ !! `inva`. In this case, `a` will not be modified, and the inverse is returned in `inva`.
+ !! Pre-allocated storage may be provided for the array of pivot indices, `pivot`. If all pre-allocated 
+ !! work spaces are provided, no internal memory allocations take place when using this interface. 
+ !!
+ !!@note The provided subroutines are intended for square matrices.
+ !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
+ !! 
+ #:for rk,rt,ri in RC_KINDS_TYPES
+ #:if rk!="xdp" 
+ module subroutine stdlib_linalg_invert_inplace_${ri}$(a,pivot,err)
+ !> Input matrix a[n,n]
+ ${rt}$, intent(inout) :: a(:,:)
+ !> [optional] Storage array for the diagonal pivot indices
+ integer(ilp), optional, intent(inout), target :: pivot(:) 
+ !> [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_invert_inplace_${ri}$
+ ! Compute the square matrix inverse of a
+ module subroutine stdlib_linalg_invert_split_${ri}$(a,inva,pivot,err)
+ !> Input matrix a[n,n].
+ ${rt}$, intent(in) :: a(:,:)
+ !> Inverse matrix a[n,n]. 
+ ${rt}$, intent(out) :: inva(:,:) 
+ !> [optional] Storage array for the diagonal pivot indices
+ integer(ilp), optional, intent(inout), target :: pivot(:)
+ !> [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_invert_split_${ri}$
+ #:endif
+ #:endfor
+ end interface invert
+
+ ! Operator interface
+ interface operator(.inv.)
+ !! version: experimental 
+ !!
+ !! Inverse operator of a square matrix
+ !! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-operator-of-a-square-matrix))
+ !!
+ !!### Summary
+ !! Operator interface for computing the inverse of a square `real` or `complex` matrix.
+ !!
+ !!### Description
+ !! 
+ !! This operator interface provides a convenient way to compute the inverse of a matrix.
+ !! Supported data types include `real` and `complex`.
+ !!
+ !!@note The provided functions are intended for square matrices.
+ !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
+ !!
+ #:for rk,rt,ri in RC_KINDS_TYPES
+ #:if rk!="xdp" 
+ module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva)
+ !> Input matrix a[n,n]
+ ${rt}$, intent(in) :: a(:,:)
+ !> Result matrix
+ ${rt}$, allocatable :: inva(:,:) 
+ end function stdlib_linalg_inverse_${ri}$_operator
+ #:endif
+ #:endfor
+ end interface operator(.inv.)
+
  ! Singular value decomposition 
  interface svd 
  !! version: experimental