diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index cb34c0cc0..a7c8c0ed3 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1459,6 +1459,142 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_inverse_function.f90!} ``` +## `pinv` - Moore-Penrose pseudo-inverse of a matrix + +### Status + +Experimental + +### Description + +This function computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix. +The pseudo-inverse, \( A^{+} \), generalizes the matrix inverse and satisfies the conditions: +- \( A \cdot A^{+} \cdot A = A \) +- \( A^{+} \cdot A \cdot A^{+} = A^{+} \) +- \( (A \cdot A^{+})^T = A \cdot A^{+} \) +- \( (A^{+} \cdot A)^T = A^{+} \cdot A \) + +The computation is based on singular value decomposition (SVD). Singular values below a relative +tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest +singular value, are treated as zero. + +### Syntax + +`b =` [[stdlib_linalg(module):pinv(interface)]] `(a, [, rtol, err])` + +### Arguments + +`a`: Shall be a rank-2, `real` or `complex` array of shape `[m, n]` containing the coefficient matrix. +It is an `intent(in)` argument. + +`rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff. +If `rtol` is not provided, the default relative tolerance is \( \text{rtol} = \text{max}(m, n) \cdot \epsilon \), +where \( \epsilon \) is the machine precision for the element type of `a`. It is an `intent(in)` 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` with shape `[n, m]`, that contains the pseudo-inverse matrix \( A^{+} \). + +Raises `LINALG_ERROR` if the underlying SVD did not converge. +Raises `LINALG_VALUE_ERROR` if `a` has invalid size. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_pseudoinverse.f90!} +``` + +## `pseudoinvert` - Moore-Penrose pseudo-inverse of a matrix + +### Status + +Experimental + +### Description + +This subroutine computes the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix. +The pseudo-inverse \( A^{+} \) is a generalization of the matrix inverse and satisfies the following properties: +- \( A \cdot A^{+} \cdot A = A \) +- \( A^{+} \cdot A \cdot A^{+} = A^{+} \) +- \( (A \cdot A^{+})^T = A \cdot A^{+} \) +- \( (A^{+} \cdot A)^T = A^{+} \cdot A \) + +The computation is based on singular value decomposition (SVD). Singular values below a relative +tolerance threshold \( \text{rtol} \cdot \sigma_{\max} \), where \( \sigma_{\max} \) is the largest +singular value, are treated as zero. + +On return, matrix `pinva` `[n, m]` will store the pseudo-inverse of `a` `[m, n]`. + +### Syntax + +`call ` [[stdlib_linalg(module):pseudoinvert(interface)]] `(a, pinva [, rtol] [, err])` + +### Arguments + +`a`: Shall be a rank-2, `real` or `complex` array containing the coefficient matrix. +It is an `intent(in)` argument. + +`pinva`: Shall be a rank-2 array of the same kind as `a`, and size equal to that of `transpose(a)`. +On output, it contains the Moore-Penrose pseudo-inverse of `a`. + +`rtol` (optional): Shall be a scalar `real` value specifying the relative tolerance for singular value cutoff. +If not provided, the default threshold is \( \text{max}(m, n) \cdot \epsilon \), where \( \epsilon \) is the +machine precision for the element type of `a`. + +`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. + +### Return value + +Computes the Moore-Penrose pseudo-inverse of the matrix \( A \), \( A^{+} \), and returns it in matrix `pinva`. + +Raises `LINALG_ERROR` if the underlying SVD did not converge. +Raises `LINALG_VALUE_ERROR` if `pinva` and `a` have degenerate or incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_pseudoinverse.f90!} +``` + +## `.pinv.` - Moore-Penrose Pseudo-Inverse operator + +### Status + +Experimental + +### Description + +This operator returns the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix \( A \). +The pseudo-inverse \( A^{+} \) is computed using Singular Value Decomposition (SVD), and singular values +below a given threshold are treated as zero. + +This interface is equivalent to the function [[stdlib_linalg(module):pinv(interface)]]. + +### Syntax + +`b = ` [[stdlib_linalg(module):operator(.pinv.)(interface)]] `a` + +### Arguments + +`a`: Shall be a rank-2 array of any `real` or `complex` kinds, with arbitrary dimensions \( m \times n \). It is an `intent(in)` argument. + +### Return value + +Returns a rank-2 array with the same type, kind, and rank as `a`, that contains the Moore-Penrose pseudo-inverse of `a`. + +If an exception occurs, or if the input matrix is degenerate (e.g., rank-deficient), the returned matrix will contain `NaN`s. +For more detailed error handling, it is recommended to use the `subroutine` or `function` interfaces. + +### Example + +```fortran +{!example/linalg/example_pseudoinverse.f90!} +``` + ## `get_norm` - Computes the vector norm of a generic-rank array. ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index ac4f0972f..2aac2a484 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -16,6 +16,7 @@ ADD_EXAMPLE(inverse_operator) ADD_EXAMPLE(inverse_function) ADD_EXAMPLE(inverse_inplace) ADD_EXAMPLE(inverse_subroutine) +ADD_EXAMPLE(pseudoinverse) ADD_EXAMPLE(outer_product) ADD_EXAMPLE(eig) ADD_EXAMPLE(eigh) diff --git a/example/linalg/example_pseudoinverse.f90 b/example/linalg/example_pseudoinverse.f90 new file mode 100644 index 000000000..6cf01642b --- /dev/null +++ b/example/linalg/example_pseudoinverse.f90 @@ -0,0 +1,32 @@ +! Matrix pseudo-inversion example: function, subroutine, and operator interfaces +program example_pseudoinverse + use stdlib_linalg, only: pinv, pseudoinvert, operator(.pinv.), linalg_state_type + implicit none(type,external) + + real :: A(15,5), Am1(5,15) + type(linalg_state_type) :: state + integer :: i, j + real, parameter :: tol = sqrt(epsilon(0.0)) + + ! Generate random matrix A (15x15) + call random_number(A) + + ! Pseudo-inverse: Function interfcae + Am1 = pinv(A, err=state) + print *, 'Max error (function) : ', maxval(abs(A-matmul(A, matmul(Am1,A)))) + + ! User threshold + Am1 = pinv(A, rtol=0.001, err=state) + print *, 'Max error (rtol=0.001): ', maxval(abs(A-matmul(A, matmul(Am1,A)))) + + ! Pseudo-inverse: Subroutine interface + call pseudoinvert(A, Am1, err=state) + + print *, 'Max error (subroutine): ', maxval(abs(A-matmul(A, matmul(Am1,A)))) + + ! Operator interface + Am1 = .pinv.A + + print *, 'Max error (operator) : ', maxval(abs(A-matmul(A, matmul(Am1,A)))) + +end program example_pseudoinverse diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ff9f39417..f5c9c96c7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,7 @@ set(fppFiles stdlib_linalg_determinant.fypp stdlib_linalg_qr.fypp stdlib_linalg_inverse.fypp + stdlib_linalg_pinv.fypp stdlib_linalg_norms.fypp stdlib_linalg_state.fypp stdlib_linalg_svd.fypp diff --git a/src/stdlib_io.fypp b/src/stdlib_io.fypp index 556d0281f..2a96f1a61 100644 --- a/src/stdlib_io.fypp +++ b/src/stdlib_io.fypp @@ -167,7 +167,7 @@ contains read (s,*,iostat=ios,iomsg=iomsg) d(i, :) if (ios/=0) then - write(msgout,1) trim(iomsg),i,trim(filename) + write(msgout,1) trim(iomsg),size(d,2),i,trim(filename) call error_stop(msg=trim(msgout)) end if @@ -178,7 +178,7 @@ contains read (s,fmt_,iostat=ios,iomsg=iomsg) d(i, :) if (ios/=0) then - write(msgout,1) trim(iomsg),i,trim(filename) + write(msgout,1) trim(iomsg),size(d,2),i,trim(filename) call error_stop(msg=trim(msgout)) end if @@ -187,7 +187,7 @@ contains close(s) - 1 format('loadtxt: error <',a,'> reading line ',i0,' of ',a,'.') + 1 format('loadtxt: error <',a,'> reading ',i0,' values from line ',i0,' of ',a,'.') end subroutine loadtxt_${t1[0]}$${k1}$ #:endfor @@ -230,14 +230,14 @@ contains iostat=ios,iomsg=iomsg) d(i, :) if (ios/=0) then - write(msgout,1) trim(iomsg),i,trim(filename) + write(msgout,1) trim(iomsg),size(d,2),i,trim(filename) call error_stop(msg=trim(msgout)) end if end do close(s) - 1 format('savetxt: error <',a,'> writing line ',i0,' of ',a,'.') + 1 format('savetxt: error <',a,'> writing ',i0,' values to line ',i0,' of ',a,'.') end subroutine savetxt_${t1[0]}$${k1}$ #:endfor diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ad8ca3fb8..682288c86 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -29,6 +29,9 @@ module stdlib_linalg public :: inv public :: invert public :: operator(.inv.) + public :: pinv + public :: pseudoinvert + public :: operator(.pinv.) public :: lstsq public :: lstsq_space public :: norm @@ -810,6 +813,131 @@ module stdlib_linalg end interface operator(.inv.) + ! Moore-Penrose Pseudo-Inverse: Function interface + interface pinv + !! version: experimental + !! + !! Pseudo-inverse of a matrix + !! ([Specification](../page/specs/stdlib_linalg.html#pinv-moore-penrose-pseudo-inverse-of-a-matrix)) + !! + !!### Summary + !! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a matrix. + !! The pseudo-inverse \( A^{+} \) is a generalization of the matrix inverse, computed for square, singular, + !! or rectangular matrices. It is defined such that it satisfies the conditions: + !! - \( A \cdot A^{+} \cdot A = A \) + !! - \( A^{+} \cdot A \cdot A^{+} = A^{+} \) + !! - \( (A \cdot A^{+})^T = A \cdot A^{+} \) + !! - \( (A^{+} \cdot A)^T = A^{+} \cdot A \) + !! + !!### Description + !! + !! This function interface provides methods that return the Moore-Penrose pseudo-inverse of a matrix. + !! Supported data types include `real` and `complex`. + !! The pseudo-inverse \( A^{+} \) is returned as a function result. The computation is based on the + !! singular value decomposition (SVD). An optional relative tolerance `rtol` is provided to control the + !! inclusion of singular values during inversion. Singular values below \( \text{rtol} \cdot \sigma_{\max} \) + !! are treated as zero, where \( \sigma_{\max} \) is the largest singular value. If `rtol` is not provided, + !! a default threshold is applied. + !! + !! Exceptions are raised in case of computational errors or invalid input, and trigger an `error stop` + !! if the state flag `err` is not provided. + !! + !!@note The provided functions are intended for both rectangular and square matrices. + !! + #:for rk,rt,ri in RC_KINDS_TYPES + module function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Relative tolerance for singular value cutoff + real(${rk}$), optional, intent(in) :: rtol + !> [optional] State return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + !> Output matrix pseudo-inverse [n,m] + ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) + end function stdlib_linalg_pseudoinverse_${ri}$ + #:endfor + end interface pinv + + ! Moore-Penrose Pseudo-Inverse: Subroutine interface + interface pseudoinvert + !! version: experimental + !! + !! Computation of the Moore-Penrose pseudo-inverse + !! ([Specification](../page/specs/stdlib_linalg.html#pseudoinvert-moore-penrose-pseudo-inverse-of-a-matrix)) + !! + !!### Summary + !! This interface provides methods for computing the Moore-Penrose pseudo-inverse of a rectangular + !! or square `real` or `complex` matrix. + !! The pseudo-inverse \( A^{+} \) generalizes the matrix inverse and satisfies the properties: + !! - \( A \cdot A^{+} \cdot A = A \) + !! - \( A^{+} \cdot A \cdot A^{+} = A^{+} \) + !! - \( (A \cdot A^{+})^T = A \cdot A^{+} \) + !! - \( (A^{+} \cdot A)^T = A^{+} \cdot A \) + !! + !!### Description + !! + !! This subroutine interface provides a way to compute the Moore-Penrose pseudo-inverse of a matrix. + !! Supported data types include `real` and `complex`. + !! Users must provide two matrices: the input matrix `a` [m,n] and the output pseudo-inverse `pinva` [n,m]. + !! The input matrix `a` is used to compute the pseudo-inverse and is not modified. The computed + !! pseudo-inverse is stored in `pinva`. The computation is based on the singular value decomposition (SVD). + !! + !! An optional relative tolerance `rtol` is used to control the inclusion of singular values in the + !! computation. Singular values below \( \text{rtol} \cdot \sigma_{\max} \) are treated as zero, + !! where \( \sigma_{\max} \) is the largest singular value. If `rtol` is not provided, a default + !! threshold is applied. + !! + !! Exceptions are raised in case of computational errors or invalid input, and trigger an `error stop` + !! if the state flag `err` is not provided. + !! + !!@note The provided subroutines are intended for both rectangular and square matrices. + !! + #:for rk,rt,ri in RC_KINDS_TYPES + module subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err) + !> Input matrix a[m,n] + ${rt}$, intent(inout) :: a(:,:) + !> Output pseudo-inverse matrix [n,m] + ${rt}$, intent(out) :: pinva(:,:) + !> [optional] Relative tolerance for singular value cutoff + real(${rk}$), optional, intent(in) :: rtol + !> [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_pseudoinvert_${ri}$ + #:endfor + end interface pseudoinvert + + ! Moore-Penrose Pseudo-Inverse: Operator interface + interface operator(.pinv.) + !! version: experimental + !! + !! Pseudo-inverse operator of a matrix + !! ([Specification](../page/specs/stdlib_linalg.html#pinv-moore-penrose-pseudo-inverse-operator)) + !! + !!### Summary + !! Operator interface for computing the Moore-Penrose pseudo-inverse of a `real` or `complex` matrix. + !! + !!### Description + !! + !! This operator interface provides a convenient way to compute the Moore-Penrose pseudo-inverse + !! of a matrix. Supported data types include `real` and `complex`. The pseudo-inverse \( A^{+} \) + !! is computed using singular value decomposition (SVD), with singular values below an internal + !! threshold treated as zero. + !! + !! For computational errors or invalid input, the function may return a matrix filled with NaNs. + !! + !!@note The provided functions are intended for both rectangular and square matrices. + !! + #:for rk,rt,ri in RC_KINDS_TYPES + module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Result pseudo-inverse matrix + ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) + end function stdlib_linalg_pinv_${ri}$_operator + #:endfor + end interface operator(.pinv.) + + ! Eigendecomposition of a square matrix: eigenvalues, and optionally eigenvectors interface eig !! version: experimental diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp new file mode 100644 index 000000000..e1d41f0a1 --- /dev/null +++ b/src/stdlib_linalg_pinv.fypp @@ -0,0 +1,132 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +submodule(stdlib_linalg) stdlib_linalg_pseudoinverse +!! Compute the (Moore-Penrose) pseudo-inverse of a matrix. + use stdlib_linalg_constants + use stdlib_linalg_blas + use stdlib_linalg_lapack + use stdlib_linalg_state + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan + implicit none(type,external) + + character(*), parameter :: this = 'pseudo-inverse' + + contains + + #:for rk,rt,ri in RC_KINDS_TYPES + + ! Compute the in-place pseudo-inverse of matrix a + module subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err) + !> Input matrix a[m,n] + ${rt}$, intent(inout) :: a(:,:) + !> Output pseudo-inverse matrix + ${rt}$, intent(out) :: pinva(:,:) + !> [optional] Relative tolerance for singular value cutoff + real(${rk}$), optional, intent(in) :: rtol + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + ! Local variables + real(${rk}$) :: tolerance,cutoff + real(${rk}$), allocatable :: s(:) + ${rt}$, allocatable :: u(:,:),vt(:,:) + type(linalg_state_type) :: err0 + integer(ilp) :: m,n,k,i,j + ${rt}$, parameter :: alpha = 1.0_${rk}$, beta = 0.0_${rk}$ + character, parameter :: H = #{if rt.startswith('complex')}# 'C' #{else}# 'T' #{endif}# + + ! Problem size + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + if (m<1 .or. n<1) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + call linalg_error_handling(err0,err) + return + end if + + if (any(shape(pinva,kind=ilp)/=[n,m])) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid pinv size:',shape(pinva),'should be',[n,m]) + call linalg_error_handling(err0,err) + return + end if + + ! Singular value threshold + tolerance = max(m,n)*epsilon(0.0_${rk}$) + + ! User threshold: fallback to default if <=0 + if (present(rtol)) then + if (rtol>0.0_${rk}$) tolerance = rtol + end if + + allocate(s(k),u(m,k),vt(k,n)) + call svd(a,s,u,vt,overwrite_a=.false.,full_matrices=.false.,err=err0) + if (err0%error()) then + err0 = linalg_state_type(this,LINALG_ERROR,'svd failure -',err0%message) + call linalg_error_handling(err0,err) + return + endif + + !> Discard singular values + cutoff = tolerance*maxval(s) + s = merge(1.0_${rk}$/s,0.0_${rk}$,s>cutoff) + + ! Get pseudo-inverse: A_pinv = V * (diag(1/s) * U^H) = V * (U * diag(1/s))^H + + ! 1) compute (U * diag(1/s)) in-place + do concurrent (i=1:m,j=1:k); u(i,j) = s(j)*u(i,j); end do + + ! 2) commutate matmul: A_pinv = V * (U * diag(1/s))^H = ((U * diag(1/s)) * V^H)^H. + ! This avoids one matrix transpose + call gemm(H, H, n, m, k, alpha, vt, k, u, m, beta, pinva, size(pinva,1,kind=ilp)) + + end subroutine stdlib_linalg_pseudoinvert_${ri}$ + + ! Function interface + module function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Relative tolerance for singular value cutoff + real(${rk}$), optional, intent(in) :: rtol + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + !> Matrix pseudo-inverse + ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) + + ! Use pointer to circumvent svd intent(inout) restriction + ${rt}$, pointer :: ap(:,:) + ap => a + + call stdlib_linalg_pseudoinvert_${ri}$(ap,pinva,rtol,err) + + end function stdlib_linalg_pseudoinverse_${ri}$ + + ! Inverse matrix operator + module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Result pseudo-inverse matrix + ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) + + type(linalg_state_type) :: err + + ! Use pointer to circumvent svd intent(inout) restriction + ${rt}$, pointer :: ap(:,:) + ap => a + + call stdlib_linalg_pseudoinvert_${ri}$(ap,pinva,err=err) + + if (err%error()) then + #:if rt.startswith('real') + pinva = ieee_value(1.0_${rk}$,ieee_quiet_nan) + #:else + pinva = 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_pinv_${ri}$_operator + + #:endfor + +end submodule stdlib_linalg_pseudoinverse diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 345b64344..49f493924 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -6,6 +6,7 @@ set( "test_linalg_eigenvalues.fypp" "test_linalg_solve.fypp" "test_linalg_inverse.fypp" + "test_linalg_pseudoinverse.fypp" "test_linalg_lstsq.fypp" "test_linalg_norm.fypp" "test_linalg_mnorm.fypp" @@ -23,6 +24,7 @@ ADDTEST(linalg_determinant) ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_inverse) +ADDTEST(linalg_pseudoinverse) ADDTEST(linalg_norm) ADDTEST(linalg_mnorm) ADDTEST(linalg_solve) @@ -30,4 +32,4 @@ ADDTEST(linalg_lstsq) ADDTEST(linalg_qr) ADDTEST(linalg_svd) ADDTEST(blas_lapack) -ADDTEST(linalg_sparse) \ No newline at end of file +ADDTEST(linalg_sparse) diff --git a/test/linalg/test_linalg_pseudoinverse.fypp b/test/linalg/test_linalg_pseudoinverse.fypp new file mode 100644 index 000000000..263ccf454 --- /dev/null +++ b/test/linalg/test_linalg_pseudoinverse.fypp @@ -0,0 +1,248 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test Moore-Penrose pseudo matrix inverse +module test_linalg_pseudoinverse + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg + use stdlib_linalg_constants + + implicit none (type,external) + private + + public :: test_pseudoinverse_matrix + + contains + + !> Matrix pseudo-inversion tests + subroutine test_pseudoinverse_matrix(tests) + !> Collertion of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in REAL_KINDS_TYPES + tests = [tests,new_unittest("${ri}$_eye_pseudoinverse",test_${ri}$_eye_pseudoinverse)] + #:endfor + #:for rk,rt,ri in RC_KINDS_TYPES + tests = [tests,new_unittest("${ri}$_square_pseudoinverse",test_${ri}$_square_pseudoinverse), & + new_unittest("${ri}$_tall_pseudoinverse",test_${ri}$_tall_pseudoinverse), & + new_unittest("${ri}$_wide_pseudoinverse",test_${ri}$_wide_pseudoinverse), & + new_unittest("${ri}$_singular_pseudoinverse",test_${ri}$_singular_pseudoinverse)] + #:endfor + + end subroutine test_pseudoinverse_matrix + + !> Invert identity matrix + #:for rk,rt,ri in REAL_KINDS_TYPES + subroutine test_${ri}$_eye_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp) :: i,j + integer(ilp), parameter :: n = 15_ilp + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + + ${rt}$ :: a(n,n),inva(n,n) + + do concurrent (i=1:n,j=1:n) + a(i,j) = merge(1.0_${rk}$,0.0_${rk}$,i==j) + end do + + !> Invert funrtion + inva = pinv(a,err=state) + + call check(error,state%ok(),'${ri}$ pseudoinverse (eye, function): '//state%print()) + if (allocated(error)) return + call check(error,all(abs(a-inva) Inverse subroutine + call pseudoinvert(a,inva,err=state) + + call check(error,state%ok(),'${ri}$ pseudoinverse (eye, subroutine): '//state%print()) + if (allocated(error)) return + call check(error,all(abs(a-inva) Operator + inva = .pinv.a + + call check(error,all(abs(a-inva) Test edge case: square matrix + subroutine test_${ri}$_square_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: n = 10 + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n, n), inva(n, n) + #:if rt.startswith('complex') + real(${rk}$) :: rea(n, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) + #:else + + call random_number(a) + #:endif + + inva = pinv(a, err=state) + call check(error,state%ok(),'${ri}$ pseudoinverse (square): '//state%print()) + if (allocated(error)) return + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (square, convergence): '//state%print()) + if (allocated(error)) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (square, convergence): '//state%print()) + if (allocated(error)) return + + end subroutine test_${ri}$_square_pseudoinverse + + !> Test edge case: tall matrix + subroutine test_${ri}$_tall_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: m = 20, n = 10 + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(m, n), inva(n, m) + #:if rt.startswith('complex') + real(${rk}$) :: rea(m, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) + #:else + + call random_number(a) + #:endif + + inva = pinv(a, err=state) + call check(error,state%ok(),'${ri}$ pseudoinverse (tall): '//state%print()) + if (allocated(error)) return + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (tall, convergence): '//state%print()) + if (allocated(error)) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (tall, convergence): '//state%print()) + if (allocated(error)) return + + end subroutine test_${ri}$_tall_pseudoinverse + + !> Test edge case: wide matrix + subroutine test_${ri}$_wide_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: m = 10, n = 20 + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(m, n), inva(n, m) + #:if rt.startswith('complex') + real(${rk}$) :: rea(m, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) + #:else + + call random_number(a) + #:endif + + inva = pinv(a, err=state) + call check(error,state%ok(),'${ri}$ pseudoinverse (wide): '//state%print()) + if (allocated(error)) return + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (wide, convergence): '//state%print()) + if (allocated(error)) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (wide, convergence): '//state%print()) + if (allocated(error)) return + + end subroutine test_${ri}$_wide_pseudoinverse + + !> Test edge case: singular matrix + subroutine test_${ri}$_singular_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: n = 10 + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) + ${rt}$ :: a(n, n), inva(n, n) + #:if rt.startswith('complex') + real(${rk}$) :: rea(n, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) + #:else + + call random_number(a) + #:endif + + ! Make the matrix singular + a(:, 1) = a(:, 2) + + inva = pinv(a, err=state) + call check(error,state%ok(),'${ri}$ pseudoinverse (singular): '//state%print()) + if (allocated(error)) return + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (singular, convergence): '//state%print()) + if (allocated(error)) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + call check(error,failed==0,'${ri}$ pseudoinverse (singular, convergence): '//state%print()) + if (allocated(error)) return + + end subroutine test_${ri}$_singular_pseudoinverse + + #:endfor + +end module test_linalg_pseudoinverse + +program test_inv + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_pseudoinverse, only : test_pseudoinverse_matrix + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_pseudoinverse", test_pseudoinverse_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 +