diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 4e8fe3b19..5a696b7ae 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -104,7 +104,8 @@ end interface axpy Note that the 128-bit functions are only provided by `stdlib` and always point to the internal implementation. Because 128-bit precision is identified as [stdlib_kinds(module):qp], initials for 128-bit procedures were labelled as `q` (quadruple-precision reals) and `w` ("wide" or quadruple-precision complex numbers). -Extended precision ([stdlib_kinds(module):xdp]) calculations are currently not supported. +Extended precision ([stdlib_kinds(module):xdp]) calculations are labelled as `x` (extended-precision reals). +and `y` (extended-precision complex numbers). ### Example @@ -779,7 +780,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. @@ -881,7 +882,7 @@ 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 @@ -889,7 +890,7 @@ This interface is equivalent to the `pure` version of determinant [[stdlib_linal ### 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. @@ -1165,3 +1166,125 @@ 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, `NaN`s will be returned. +For fine-grained error control in case of singular matrices prefer the `subroutine` and the `function` +interfaces. + + +### 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!} +``` + diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 3606f0a7e..ecbb1a9cc 100644 --- a/example/linalg/CMakeLists.txt +++ b/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(eig) ADD_EXAMPLE(eigh) diff --git a/example/linalg/example_inverse_function.f90 b/example/linalg/example_inverse_function.f90 new file mode 100644 index 000000000..a59b0ebdb --- /dev/null +++ b/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 diff --git a/example/linalg/example_inverse_inplace.f90 b/example/linalg/example_inverse_inplace.f90 new file mode 100644 index 000000000..9a3816a34 --- /dev/null +++ b/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 diff --git a/example/linalg/example_inverse_operator.f90 b/example/linalg/example_inverse_operator.f90 new file mode 100644 index 000000000..a4f71c75d --- /dev/null +++ b/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 diff --git a/example/linalg/example_inverse_subroutine.f90 b/example/linalg/example_inverse_subroutine.f90 new file mode 100644 index 000000000..1c75d04f0 --- /dev/null +++ b/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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index fa1831e98..579b70c72 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -30,6 +30,7 @@ set(fppFiles stdlib_linalg_eigenvalues.fypp stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp + stdlib_linalg_inverse.fypp stdlib_linalg_state.fypp stdlib_linalg_svd.fypp stdlib_optval.fypp diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index fe8fb0bc0..9cb7a13a7 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -24,6 +24,9 @@ module stdlib_linalg public :: eigvals public :: eigvalsh public :: eye + public :: inv + public :: invert + public :: operator(.inv.) public :: lstsq public :: lstsq_space public :: solve @@ -544,6 +547,114 @@ module stdlib_linalg #:endfor end interface + ! Matrix Inverse: 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. + !! + #:for rk,rt,ri in RC_KINDS_TYPES + 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}$ + #:endfor + end interface inv + + ! Matrix Inverse: Subroutine interface - in-place inversion + 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. + !! + #:for rk,rt,ri in RC_KINDS_TYPES + 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}$ + #:endfor + end interface invert + + ! Matrix Inverse: 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. + !! + #:for rk,rt,ri in RC_KINDS_TYPES + 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 + #:endfor + end interface operator(.inv.) + + ! Eigendecomposition of a square matrix: eigenvalues, and optionally eigenvectors interface eig !! version: experimental @@ -747,7 +858,6 @@ module stdlib_linalg #:endfor end interface eigvalsh - ! Singular value decomposition interface svd !! version: experimental diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp new file mode 100644 index 000000000..38ad3dda7 --- /dev/null +++ b/src/stdlib_linalg_inverse.fypp @@ -0,0 +1,179 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +submodule (stdlib_linalg) stdlib_linalg_inverse +!! Compute inverse of a square matrix + use stdlib_linalg_constants + use stdlib_linalg_lapack, only: getri,getrf,stdlib_ilaenv + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan + implicit none(type,external) + + character(*), parameter :: this = 'inverse' + + contains + + elemental subroutine handle_getri_info(info,lda,n,err) + integer(ilp), intent(in) :: info,lda,n + type(linalg_state_type), intent(out) :: err + + ! Process output + select case (info) + case (0) + ! Success + case (:-1) + err = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=',[lda,n]) + case (1:) + ! Matrix is singular + err = linalg_state_type(this,LINALG_ERROR,'singular matrix') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select + + end subroutine handle_getri_info + + #:for rk,rt,ri in RC_KINDS_TYPES + ! Compute the in-place square matrix inverse of a + module subroutine stdlib_linalg_invert_inplace_${ri}$(a,pivot,err) + !> Input matrix a[n,n]. On return, A is destroyed and replaced by the inverse + ${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 + + !> Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: lda,n,info,nb,lwork,npiv + integer(ilp), pointer :: ipiv(:) + ${rt}$, allocatable :: work(:) + + !> Problem sizes + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + + ! Has a pre-allocated pivots storage array been provided? + if (present(pivot)) then + ipiv => pivot + else + allocate(ipiv(n)) + endif + npiv = size(ipiv,kind=ilp) + + if (lda<1 .or. n<1 .or. lda/=n .or. npiv 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 + + type(linalg_state_type) :: err0 + integer(ilp) :: sa(2),sinva(2) + + sa = shape(a,kind=ilp) + sinva = shape(inva,kind=ilp) + + if (any(sa/=sinva)) then + + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',sa,' inva=',sinva) + + else + + !> Copy data in + inva = a + + !> Compute matrix inverse + call stdlib_linalg_invert_inplace_${ri}$(inva,err=err0) + + end if + + ! Process output and return + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_invert_split_${ri}$ + + ! Invert matrix in place + 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 + + !> Allocate with copy + allocate(inva,source=a) + + !> Compute matrix inverse + call stdlib_linalg_invert_inplace_${ri}$(inva,err=err) + + end function stdlib_linalg_inverse_${ri}$ + + ! Inverse matrix operator + module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva) + !> Input matrix a[n,n] + ${rt}$, intent(in) :: a(:,:) + !> Result matrix + ${rt}$, allocatable :: inva(:,:) + + type(linalg_state_type) :: err + + ! Provide an error handler to return NaNs on issues + inva = stdlib_linalg_inverse_${ri}$(a,err=err) + + ! Return NaN on issues + if (err%error()) then + if (allocated(inva)) deallocate(inva) + allocate(inva(size(a,1,kind=ilp),size(a,2,kind=ilp))) + + #:if rt.startswith('complex') + inva = ieee_value(1.0_${rk}$,ieee_quiet_nan) + #:else + inva = cmplx(ieee_value(1.0_${rk}$,ieee_quiet_nan), & + ieee_value(1.0_${rk}$,ieee_quiet_nan), kind=${rk}$) + #:endif + endif + + end function stdlib_linalg_inverse_${ri}$_operator + + #:endfor + +end submodule stdlib_linalg_inverse diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index cff60532d..3add198ff 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -4,6 +4,7 @@ set( "test_blas_lapack.fypp" "test_linalg_eigenvalues.fypp" "test_linalg_solve.fypp" + "test_linalg_inverse.fypp" "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" "test_linalg_svd.fypp" @@ -15,6 +16,7 @@ ADDTEST(linalg) ADDTEST(linalg_determinant) ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) +ADDTEST(linalg_inverse) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(linalg_svd) diff --git a/test/linalg/test_linalg_inverse.fypp b/test/linalg/test_linalg_inverse.fypp new file mode 100644 index 000000000..679a586c5 --- /dev/null +++ b/test/linalg/test_linalg_inverse.fypp @@ -0,0 +1,321 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test inverse matrix operator +module test_linalg_inverse + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: inv,invert,operator(.inv.),eye + use stdlib_linalg_state, only: linalg_state_type,LINALG_ERROR + + implicit none (type,external) + private + + public :: test_inverse_matrix + + contains + + !> Matrix inversion tests + subroutine test_inverse_matrix(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in RC_KINDS_TYPES + tests = [tests,new_unittest("${ri}$_eye_inverse",test_${ri}$_eye_inverse), & + new_unittest("${ri}$_singular_inverse",test_${ri}$_singular_inverse), & + new_unittest("${ri}$_random_spd_inverse",test_${ri}$_random_spd_inverse)] + #:endfor + + end subroutine test_inverse_matrix + + #:for rk,rt,ri in REAL_KINDS_TYPES + !> Invert real identity matrix + subroutine test_${ri}$_eye_inverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp), parameter :: n = 25_ilp + ${rt}$ :: a(n,n),inva(n,n) + + a = eye(n) + + !> Inverse function + inva = inv(a,err=state) + call check(error,state%ok(),'inverse_${ri}$_eye (function): '//state%print()) + if (allocated(error)) return + + call check(error,all(abs(a-inva) Inverse subroutine: split + call invert(a,inva,err=state) + + call check(error,state%ok(),'inverse_${ri}$_eye (subroutine): '//state%print()) + if (allocated(error)) return + + call check(error,all(abs(a-inva) Inverse subroutine in-place + call invert(a,err=state) + + call check(error,state%ok(),'inverse_${ri}$_eye (in-place): '//state%print()) + if (allocated(error)) return + + call check(error,all(abs(a-inva) Invert singular matrix + subroutine test_${ri}$_singular_inverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: err + + integer(ilp), parameter :: n = 25_ilp + ${rt}$ :: a(n,n) + + a = eye(n) + + !> Make rank-deficient + a(12,12) = 0 + + !> Inverse function + call invert(a,err=err) + call check(error,err%state==LINALG_ERROR,'singular ${rt}$ inverse returned '//err%print()) + if (allocated(error)) return + + end subroutine test_${ri}$_singular_inverse + + !> Create a random symmetric positive definite matrix + function random_spd_matrix_${ri}$(n) result(A) + integer(ilp), intent(in) :: n + ${rt}$ :: A(n,n) + + ${rt}$, parameter :: one = 1.0_${rk}$ + ${rt}$, parameter :: half = 0.5_${rk}$ + + !> Initialize with randoms + call random_number(A) + + !> Make symmetric + A = half*(A+transpose(A)) + + !> Add diagonally dominant part + A = A + n*eye(n) + + end function random_spd_matrix_${ri}$ + + !> Test random symmetric positive-definite matrix + subroutine test_${ri}$_random_spd_inverse(error) + type(error_type), allocatable, intent(out) :: error + + !> Solution tolerance + ${rt}$, parameter :: tol = sqrt(epsilon(0.0_${rk}$)) + + !> Local variables + integer(ilp), parameter :: n = 5_ilp + type(linalg_state_type) :: state + ${rt}$ :: A(n,n),Am1(n,n) + + !> Generate random SPD matrix + A = random_spd_matrix_${ri}$(n) + + !> Invert matrix + call invert(A,Am1,err=state) + + !> Check result + call check(error,state%ok(),'random SPD matrix (${rk}$): '//state%print()) + if (allocated(error)) return + + call check(error,all(abs(matmul(Am1,A)-eye(n)) Invert complex identity matrix + #:for ck,ct,ci in CMPLX_KINDS_TYPES + subroutine test_${ci}$_eye_inverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp) :: i,j,failed + integer(ilp), parameter :: n = 25_ilp + + ${ct}$ :: a(n,n),copya(n,n),inva(n,n) + + do concurrent (i=1:n,j=1:n) + a(i,j) = merge((1.0_${ck}$,1.0_${ck}$),(0.0_${ck}$,0.0_${ck}$),i==j) + end do + copya = a + + !> The inverse of a complex diagonal matrix has conjg(z_ii)/abs(z_ii)^2 on the diagonal + inva = inv(a,err=state) + + call check(error,state%ok(),'inverse_${ci}$_eye (function): '//state%print()) + if (allocated(error)) return + + failed = 0 + do i=1,n + do j=1,n + if (.not.is_diagonal_inverse(a(i,j),inva(i,j),i,j)) failed = failed+1 + end do + end do + + call check(error,failed==0,'inverse_${ci}$_eye (function): data converged') + if (allocated(error)) return + + !> Inverse subroutine + call invert(copya,err=state) + + call check(error,state%ok(),'inverse_${ci}$_eye (subroutine): '//state%print()) + if (allocated(error)) return + + failed = 0 + do i=1,n + do j=1,n + if (.not.is_diagonal_inverse(a(i,j),copya(i,j),i,j)) failed = failed+1 + end do + end do + + call check(error,failed==0,'inverse_${ci}$_eye (subroutine): data converged') + if (allocated(error)) return + + contains + + elemental logical function is_diagonal_inverse(aij,invaij,i,j) + ${ct}$, intent(in) :: aij,invaij + integer(ilp), intent(in) :: i,j + if (i/=j) then + is_diagonal_inverse = max(abs(aij),abs(invaij)) Create a random symmetric positive definite matrix + function random_spd_matrix_${ci}$(n) result(A) + integer(ilp), intent(in) :: n + ${ct}$ :: A(n,n) + + ${ct}$, parameter :: half = (0.5_${ck}$,0.0_${ck}$) + real(${ck}$) :: reA(n,n),imA(n,n) + integer(ilp) :: i + + !> Initialize with randoms + call random_number(reA) + call random_number(imA) + + A = cmplx(reA,imA,kind=${ck}$) + + !> Make symmetric + A = half*(A+transpose(A)) + + !> Add diagonally dominant part + forall(i=1:n) A(i,i) = A(i,i) + n*(1.0_${ck}$,0.0_${ck}$) + + end function random_spd_matrix_${ci}$ + + !> Test random symmetric positive-definite matrix + subroutine test_${ci}$_random_spd_inverse(error) + type(error_type), allocatable, intent(out) :: error + + !> Local variables + integer(ilp) :: failed,i,j + integer(ilp), parameter :: n = 5_ilp + type(linalg_state_type) :: state + ${ct}$ :: A(n,n),Am1(n,n),AA(n,n) + + !> Generate random SPD matrix + A = random_spd_matrix_${ci}$(n) + + !> Invert matrix + call invert(A,Am1,err=state) + + !> Check result + call check(error,state%ok(),'random complex SPD matrix (${ck}$): '//state%print()) + if (allocated(error)) return + + failed = 0 + AA = matmul(A,Am1) + do i=1,n + do j=1,n + if (.not.is_complex_inverse(AA(i,j),i,j)) failed = failed+1 + end do + end do + + call check(error,failed==0,'inverse_${ci}$_eye (subroutine): data converged') + if (allocated(error)) return + + contains + + elemental logical function is_complex_inverse(aij,i,j) + ${ct}$, intent(in) :: aij + integer(ilp), intent(in) :: i,j + real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) + if (i/=j) then + is_complex_inverse = abs(aij) Invert singular matrix + subroutine test_${ci}$_singular_inverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: err + + integer(ilp), parameter :: n = 25_ilp + ${ct}$ :: a(n,n) + + a = (0.0_${ck}$,0.0_${ck}$) + + !> Inverse function + call invert(a,err=err) + call check(error,err%state==LINALG_ERROR,'singular ${ct}$ inverse returned '//err%print()) + if (allocated(error)) return + + end subroutine test_${ci}$_singular_inverse + + #:endfor + +end module test_linalg_inverse + +program test_inv + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_inverse, only : test_inverse_matrix + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_inverse", test_inverse_matrix) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_inv