linalg: Matrix Inverse

doc/specs/stdlib_linalg.md

@@ -777,7 +777,7 @@ 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 `min(m,n)`, returning the list of singular values `s(i)>=cond*maxval(s)` from the internal SVD, 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 +879,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 +989,121 @@ 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  [[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`.
+
+If an exception occurred on input errors, or singular matrix, NaNs will be returned.
+
+
+### Example
+
+```fortran
+{!example/linalg/example_inverse_operator.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`, that contains the diagonal pivot indices on return. It is an `intent(inout)` argument. 
+
+`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
+
+### Return value
+
+Computes the inverse of the matrix \( A \), \(A^{-1}\, and returns it either in \( A \) or in another matrix.
+
+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_inverse_inplace.f90!}
+```
+
+```fortran
+{!example/linalg/example_inverse_subroutine.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. It 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_inverse_function.f90!}
+```

example/linalg/CMakeLists.txt

@@ -12,6 +12,10 @@ ADD_EXAMPLE(is_skew_symmetric)
 ADD_EXAMPLE(is_square)
 ADD_EXAMPLE(is_symmetric)
 ADD_EXAMPLE(is_triangular)
+ADD_EXAMPLE(inverse_operator)
+ADD_EXAMPLE(inverse_function)
+ADD_EXAMPLE(inverse_inplace)
+ADD_EXAMPLE(inverse_subroutine)
 ADD_EXAMPLE(outer_product)
 ADD_EXAMPLE(trace)
 ADD_EXAMPLE(state1)

example/linalg/example_inverse_function.f90

@@ -0,0 +1,22 @@
+! Matrix inversion example: function interface
+program example_inverse_function
+  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_inverse_function

example/linalg/example_inverse_inplace.f90

@@ -0,0 +1,23 @@
+! Matrix inversion example: in-place inversion
+program example_inverse_inplace
+  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_inverse_inplace

example/linalg/example_inverse_operator.f90

@@ -0,0 +1,22 @@
+! Matrix inversion example: operator interface
+program example_inverse_operator
+  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_inverse_operator

example/linalg/example_inverse_subroutine.f90

@@ -0,0 +1,22 @@
+! Matrix inversion example: subroutine interface
+program example_inverse_subroutine
+  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] ))
+
+  ! Invert matrix
+  call invert(A,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_inverse_subroutine

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,122 @@ 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`. On input errors or singular matrix, 
+    !! NaNs will be returned.
+    !!
+    !!@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