From f5de00f5b5337e6af7080c167ae43c6e78e26596 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 16:18:04 +0100 Subject: [PATCH 01/18] copy in code --- src/CMakeLists.txt | 1 + src/stdlib_linalg_pinv.fypp | 152 ++++++++++++++++++++++++++++++++++++ 2 files changed, 153 insertions(+) create mode 100644 src/stdlib_linalg_pinv.fypp 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_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp new file mode 100644 index 000000000..61bd72b7f --- /dev/null +++ b/src/stdlib_linalg_pinv.fypp @@ -0,0 +1,152 @@ +#:include "common.fypp" +! Compute the (Moore-Penrose) pseudo-inverse of a matrix. +module stdlib_linalg_pseudoinverse + use stdlib_linalg_constants + use stdlib_linalg_blas + use stdlib_linalg_lapack + use stdlib_linalg_state + use stdlib_linalg_svd, only: svd + use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + implicit none(type,external) + private + + !> Pseudo-inverse: Function interface + public :: pinv + !> Pseudo-inverse: Subroutine interface (pre-allocated) + public :: pseudoinvert + !> Operator interface: .pinv.A returns the pseudo-inverse of A + public :: operator(.pinv.) + + ! Function interface + interface pinv + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_pseudoinverse_${ri}$ + #:endfor + end interface pinv + + ! Subroutine interface + interface pseudoinvert + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_pseudoinvert_${ri}$ + #:endfor + end interface pseudoinvert + + ! Operator interface + interface operator(.pinv.) + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_pinv_${ri}$_operator + #:endfor + end interface operator(.pinv.) + + character(*), parameter :: this = 'pseudo-inverse' + + contains + + #:for rk,rt,ri in ALL_KINDS_TYPES + + ! Compute the in-place pseudo-inverse of matrix a + subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err) + !> Input matrix a[m,n] + ${rt}$, intent(inout) :: a(:,:) + !> Output pseudo-inverse matrix + ${rt}$, intent(inout) :: pinva(:,:) + !> [optional] .... + real(${rk}$), optional, intent(in) :: rtol + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), optional, intent(out) :: err + + ! Local variables + real(${rk}$) :: tolerance,cutoff + real(${rk}$), allocatable :: s(:) + ${rt}$, allocatable :: u(:,:),vt(:,:) + type(linalg_state) :: err0 + integer(ilp) :: m,n,k,i,j + + ! 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(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(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(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/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 + forall (i=1:m,j=1:k) u(i,j) = s(j)*u(i,j) + + ! 2) commutate matmul: A_pinv = V^H * (U * diag(1/s))^H = ((U * diag(1/s)) * V^H)^H. + ! This avoids one matrix transpose + #:if rt.startswith('complex') + pinva = conjg(transpose(matmul(u,vt))) + #:else + pinva = transpose(matmul(u,vt)) + #:endif + + end subroutine stdlib_linalg_pseudoinvert_${ri}$ + + ! Function interface + function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] .... + real(${rk}$), optional, intent(in) :: rtol + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), 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 + function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) + !> Input matrix a[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Result matrix + ${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) + + end function stdlib_linalg_pinv_${ri}$_operator + + #:endfor + +end module stdlib_linalg_pseudoinverse From 9d350849baf723297700739aa183ad356c637f05 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 16:19:08 +0100 Subject: [PATCH 02/18] copy in tests --- test/linalg/CMakeLists.txt | 4 +- test/linalg/test_linalg_pseudoinverse.fypp | 208 +++++++++++++++++++++ 2 files changed, 211 insertions(+), 1 deletion(-) create mode 100644 test/linalg/test_linalg_pseudoinverse.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index ad41e5bb0..9d1fa5c6c 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_determinant.fypp" @@ -22,10 +23,11 @@ ADDTEST(linalg_determinant) ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_inverse) +ADDTEST(linalg_pseudoinverse) ADDTEST(linalg_norm) ADDTEST(linalg_solve) 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..6056ef7e2 --- /dev/null +++ b/test/linalg/test_linalg_pseudoinverse.fypp @@ -0,0 +1,208 @@ +#:include "common.fypp" +! Test inverse matrix +module test_linalg_pseudoinverse + use stdlib_linalg_interface + + implicit none (type,external) + + contains + + !> Matrix inversion tests + subroutine test_pseudoinverse_matrix(error) + logical, intent(out) :: error + + real :: t0,t1 + + call cpu_time(t0) + + #:for rk,rt,ri in REAL_KINDS_TYPES + call test_${ri}$_eye_pseudoinverse(error) + if (error) return + #:endfor + #:for ck,ct,ci in ALL_KINDS_TYPES + call test_${ci}$_square_pseudoinverse(error) + if (error) return + call test_${ci}$_tall_pseudoinverse(error) + if (error) return + call test_${ci}$_wide_pseudoinverse(error) + if (error) return + call test_${ci}$_singular_pseudoinverse(error) + if (error) return + + #:endfor + + call cpu_time(t1) + + print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error) + + 1 format('Pseudo-Inverse matrix tests completed in ',f9.4,' milliseconds, result=',a) + + end subroutine test_pseudoinverse_matrix + + !> Invert identity matrix + #:for rk,rt,ri in REAL_KINDS_TYPES + subroutine test_${ri}$_eye_pseudoinverse(error) + logical, intent(out) :: error + + type(linalg_state) :: state + + integer(ilp) :: i,j + integer(ilp), parameter :: n = 15_ilp + real(${rk}$), parameter :: tol = 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 function + inva = pinv(a,err=state) + error = state%error() .or. .not.all(abs(a-inva) Inverse subroutine + call pseudoinvert(a,inva,err=state) + error = state%error() .or. .not.all(abs(a-inva) Operator + inva = .pinv.a + error = .not.all(abs(a-inva) Test edge case: square matrix + subroutine test_${ci}$_square_pseudoinverse(error) + logical, intent(out) :: error + + type(linalg_state) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: n = 10 + real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) + ${ct}$ :: a(n, n), inva(n, n) + #:if ct.startswith('complex') + real(${ck}$) :: rea(n, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$) + #:else + + call random_number(a) + #:endif + + inva = pinv(a, err=state) + error = state%error(); if (error) return + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + error = failed > 0; if (error) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + error = failed > 0; if (error) return + + end subroutine test_${ci}$_square_pseudoinverse + + !> Test edge case: tall matrix + subroutine test_${ci}$_tall_pseudoinverse(error) + logical, intent(out) :: error + + type(linalg_state) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: m = 20, n = 10 + real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) + ${ct}$ :: a(m, n), inva(n, m) + #:if ct.startswith('complex') + real(${ck}$) :: rea(m, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$) + #:else + + call random_number(a) + #:endif + + inva = pinv(a, err=state) + error = state%error(); if (error) return + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + error = failed > 0; if (error) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + error = failed > 0; if (error) return + + end subroutine test_${ci}$_tall_pseudoinverse + + !> Test edge case: wide matrix + subroutine test_${ci}$_wide_pseudoinverse(error) + logical, intent(out) :: error + + type(linalg_state) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: m = 10, n = 20 + real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) + ${ct}$ :: a(m, n), inva(n, m) + #:if ct.startswith('complex') + real(${ck}$) :: rea(m, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$) + #:else + + call random_number(a) + #:endif + + inva = pinv(a, err=state) + error = state%error(); if (error) return + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + error = failed > 0; if (error) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + error = failed > 0; if (error) return + + end subroutine test_${ci}$_wide_pseudoinverse + + !> Test edge case: singular matrix + subroutine test_${ci}$_singular_pseudoinverse(error) + logical, intent(out) :: error + + type(linalg_state) :: state + + integer(ilp) :: failed + integer(ilp), parameter :: n = 10 + real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) + ${ct}$ :: a(n, n), inva(n, n) + #:if ct.startswith('complex') + real(${ck}$) :: rea(n, n, 2) + + call random_number(rea) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${ck}$) + #:else + + call random_number(a) + #:endif + + ! Make the matrix singular + a(:, 1) = a(:, 2) + + inva = pinv(a, err=state) + + failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) + error = failed > 0; if (error) return + + failed = count(abs(inva - matmul(inva, matmul(a, inva))) > tol) + error = failed > 0; if (error) return + + end subroutine test_${ci}$_singular_pseudoinverse + + #:endfor + +end module test_linalg_pseudoinverse + + From c54b95919cc56602f77ae7799a051d87ea9289f3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 16:47:10 +0100 Subject: [PATCH 03/18] adjust templates --- src/stdlib_linalg_pinv.fypp | 25 +-- test/linalg/test_linalg_pseudoinverse.fypp | 209 ++++++++++++--------- 2 files changed, 138 insertions(+), 96 deletions(-) diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index 61bd72b7f..5cc7d4cfa 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -1,11 +1,12 @@ #:include "common.fypp" -! Compute the (Moore-Penrose) pseudo-inverse of a matrix. +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module 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 stdlib_linalg_svd, only: svd + use stdlib_linalg, only: svd use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit implicit none(type,external) private @@ -19,21 +20,21 @@ module stdlib_linalg_pseudoinverse ! Function interface interface pinv - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES module procedure stdlib_linalg_pseudoinverse_${ri}$ #:endfor end interface pinv ! Subroutine interface interface pseudoinvert - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES module procedure stdlib_linalg_pseudoinvert_${ri}$ #:endfor end interface pseudoinvert ! Operator interface interface operator(.pinv.) - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES module procedure stdlib_linalg_pinv_${ri}$_operator #:endfor end interface operator(.pinv.) @@ -42,7 +43,7 @@ module stdlib_linalg_pseudoinverse contains - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES ! Compute the in-place pseudo-inverse of matrix a subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err) @@ -53,13 +54,13 @@ module stdlib_linalg_pseudoinverse !> [optional] .... real(${rk}$), optional, intent(in) :: rtol !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state), optional, intent(out) :: err + 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) :: err0 + type(linalg_state_type) :: err0 integer(ilp) :: m,n,k,i,j ! Problem size @@ -67,13 +68,13 @@ module stdlib_linalg_pseudoinverse n = size(a,2,kind=ilp) k = min(m,n) if (m<1 .or. n<1) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + 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(this,LINALG_VALUE_ERROR,'invalid pinv size:',shape(pinva),'should be',[n,m]) + 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 @@ -89,7 +90,7 @@ module stdlib_linalg_pseudoinverse 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(this,LINALG_ERROR,'svd failure -',err0%message) + err0 = linalg_state_type(this,LINALG_ERROR,'svd failure -',err0%message) call linalg_error_handling(err0,err) return endif @@ -120,7 +121,7 @@ module stdlib_linalg_pseudoinverse !> [optional] .... real(${rk}$), optional, intent(in) :: rtol !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state), optional, intent(out) :: err + type(linalg_state_type), optional, intent(out) :: err !> Matrix pseudo-inverse ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) diff --git a/test/linalg/test_linalg_pseudoinverse.fypp b/test/linalg/test_linalg_pseudoinverse.fypp index 6056ef7e2..e9524e39e 100644 --- a/test/linalg/test_linalg_pseudoinverse.fypp +++ b/test/linalg/test_linalg_pseudoinverse.fypp @@ -1,50 +1,44 @@ #:include "common.fypp" -! Test inverse matrix +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test Moore-Penrose pseudo matrix inverse module test_linalg_pseudoinverse - use stdlib_linalg_interface + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg + use stdlib_linalg_constants + use stdlib_linalg_pseudoinverse implicit none (type,external) + private + + public :: test_pseudoinverse_matrix contains - !> Matrix inversion tests - subroutine test_pseudoinverse_matrix(error) - logical, intent(out) :: error - - real :: t0,t1 - - call cpu_time(t0) + !> 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 - call test_${ri}$_eye_pseudoinverse(error) - if (error) return + tests = [tests,new_unittest("${ri}$_eye_pseudoinverse",test_${ri}$_eye_pseudoinverse)] #:endfor - #:for ck,ct,ci in ALL_KINDS_TYPES - call test_${ci}$_square_pseudoinverse(error) - if (error) return - call test_${ci}$_tall_pseudoinverse(error) - if (error) return - call test_${ci}$_wide_pseudoinverse(error) - if (error) return - call test_${ci}$_singular_pseudoinverse(error) - if (error) return - + #: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 - call cpu_time(t1) - - print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error) - - 1 format('Pseudo-Inverse matrix tests completed in ',f9.4,' milliseconds, result=',a) - end subroutine test_pseudoinverse_matrix !> Invert identity matrix #:for rk,rt,ri in REAL_KINDS_TYPES subroutine test_${ri}$_eye_pseudoinverse(error) - logical, intent(out) :: error + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp) :: i,j integer(ilp), parameter :: n = 15_ilp @@ -56,133 +50,151 @@ module test_linalg_pseudoinverse a(i,j) = merge(1.0_${rk}$,0.0_${rk}$,i==j) end do - !> Invert function + !> Invert funrtion inva = pinv(a,err=state) - error = state%error() .or. .not.all(abs(a-inva) Inverse subroutine call pseudoinvert(a,inva,err=state) - error = state%error() .or. .not.all(abs(a-inva) Operator inva = .pinv.a - error = .not.all(abs(a-inva) Test edge case: square matrix - subroutine test_${ci}$_square_pseudoinverse(error) - logical, intent(out) :: error + subroutine test_${ri}$_square_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp) :: failed integer(ilp), parameter :: n = 10 - real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) - ${ct}$ :: a(n, n), inva(n, n) - #:if ct.startswith('complex') - real(${ck}$) :: rea(n, n, 2) + real(${rk}$), parameter :: tol = 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=${ck}$) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) #:else call random_number(a) #:endif inva = pinv(a, err=state) - error = state%error(); if (error) return + call check(error,state%ok(),'${ri}$ pseudoinverse (square): '//state%print()) + if (allocated(error)) return failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) - error = failed > 0; if (error) return + 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) - error = failed > 0; if (error) return + call check(error,failed==0,'${ri}$ pseudoinverse (square, convergence): '//state%print()) + if (allocated(error)) return - end subroutine test_${ci}$_square_pseudoinverse + end subroutine test_${ri}$_square_pseudoinverse !> Test edge case: tall matrix - subroutine test_${ci}$_tall_pseudoinverse(error) - logical, intent(out) :: error + subroutine test_${ri}$_tall_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp) :: failed integer(ilp), parameter :: m = 20, n = 10 - real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) - ${ct}$ :: a(m, n), inva(n, m) - #:if ct.startswith('complex') - real(${ck}$) :: rea(m, n, 2) + real(${rk}$), parameter :: tol = 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=${ck}$) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) #:else call random_number(a) #:endif inva = pinv(a, err=state) - error = state%error(); if (error) return + call check(error,state%ok(),'${ri}$ pseudoinverse (tall): '//state%print()) + if (allocated(error)) return failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) - error = failed > 0; if (error) return + 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) - error = failed > 0; if (error) return + call check(error,failed==0,'${ri}$ pseudoinverse (tall, convergence): '//state%print()) + if (allocated(error)) return - end subroutine test_${ci}$_tall_pseudoinverse + end subroutine test_${ri}$_tall_pseudoinverse !> Test edge case: wide matrix - subroutine test_${ci}$_wide_pseudoinverse(error) - logical, intent(out) :: error + subroutine test_${ri}$_wide_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp) :: failed integer(ilp), parameter :: m = 10, n = 20 - real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) - ${ct}$ :: a(m, n), inva(n, m) - #:if ct.startswith('complex') - real(${ck}$) :: rea(m, n, 2) + real(${rk}$), parameter :: tol = 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=${ck}$) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) #:else call random_number(a) #:endif inva = pinv(a, err=state) - error = state%error(); if (error) return + call check(error,state%ok(),'${ri}$ pseudoinverse (wide): '//state%print()) + if (allocated(error)) return failed = count(abs(a - matmul(a, matmul(inva, a))) > tol) - error = failed > 0; if (error) return + 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) - error = failed > 0; if (error) return + call check(error,failed==0,'${ri}$ pseudoinverse (wide, convergence): '//state%print()) + if (allocated(error)) return - end subroutine test_${ci}$_wide_pseudoinverse + end subroutine test_${ri}$_wide_pseudoinverse !> Test edge case: singular matrix - subroutine test_${ci}$_singular_pseudoinverse(error) - logical, intent(out) :: error + subroutine test_${ri}$_singular_pseudoinverse(error) + type(error_type), allocatable, intent(out) :: error - type(linalg_state) :: state + type(linalg_state_type) :: state integer(ilp) :: failed integer(ilp), parameter :: n = 10 - real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) - ${ct}$ :: a(n, n), inva(n, n) - #:if ct.startswith('complex') - real(${ck}$) :: rea(n, n, 2) + real(${rk}$), parameter :: tol = 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=${ck}$) + a = cmplx(rea(:, :, 1), rea(:, :, 2), kind=${rk}$) #:else call random_number(a) @@ -192,17 +204,46 @@ module test_linalg_pseudoinverse 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) - error = failed > 0; if (error) return + 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) - error = failed > 0; if (error) return - - end subroutine test_${ci}$_singular_pseudoinverse + 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 From 37612f6d329980104851e8ac88cbbb911b0b5539 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:03:07 +0100 Subject: [PATCH 04/18] submodule --- src/stdlib_linalg.fypp | 47 ++++++++++++++++++++++ src/stdlib_linalg_pinv.fypp | 39 +++--------------- test/linalg/test_linalg_pseudoinverse.fypp | 1 - 3 files changed, 52 insertions(+), 35 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ad36ad677..e6d415ae5 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 @@ -809,6 +812,50 @@ module stdlib_linalg end interface operator(.inv.) + ! Moose-Penrose Pseudo-Inverse: Function interface + interface pinv + #: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] .... + 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)) + end function stdlib_linalg_pseudoinverse_${ri}$ + #:endfor + end interface pinv + + ! Matrix Inverse: Subroutine interface - in-place inversion + interface pseudoinvert + #: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 + ${rt}$, intent(inout) :: pinva(:,:) + !> [optional] .... + 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 + #:endfor + end interface pseudoinvert + + ! Operator interface + interface operator(.pinv.) + #: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 matrix + ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) + end function + #: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 index 5cc7d4cfa..16bdcb68d 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -1,6 +1,6 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -module stdlib_linalg_pseudoinverse +submodule(stdlib_linalg) stdlib_linalg_pseudoinverse !! Compute the (Moore-Penrose) pseudo-inverse of a matrix. use stdlib_linalg_constants use stdlib_linalg_blas @@ -9,35 +9,6 @@ module stdlib_linalg_pseudoinverse use stdlib_linalg, only: svd use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit implicit none(type,external) - private - - !> Pseudo-inverse: Function interface - public :: pinv - !> Pseudo-inverse: Subroutine interface (pre-allocated) - public :: pseudoinvert - !> Operator interface: .pinv.A returns the pseudo-inverse of A - public :: operator(.pinv.) - - ! Function interface - interface pinv - #:for rk,rt,ri in RC_KINDS_TYPES - module procedure stdlib_linalg_pseudoinverse_${ri}$ - #:endfor - end interface pinv - - ! Subroutine interface - interface pseudoinvert - #:for rk,rt,ri in RC_KINDS_TYPES - module procedure stdlib_linalg_pseudoinvert_${ri}$ - #:endfor - end interface pseudoinvert - - ! Operator interface - interface operator(.pinv.) - #:for rk,rt,ri in RC_KINDS_TYPES - module procedure stdlib_linalg_pinv_${ri}$_operator - #:endfor - end interface operator(.pinv.) character(*), parameter :: this = 'pseudo-inverse' @@ -46,7 +17,7 @@ module stdlib_linalg_pseudoinverse #:for rk,rt,ri in RC_KINDS_TYPES ! Compute the in-place pseudo-inverse of matrix a - subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err) + module subroutine stdlib_linalg_pseudoinvert_${ri}$(a,pinva,rtol,err) !> Input matrix a[m,n] ${rt}$, intent(inout) :: a(:,:) !> Output pseudo-inverse matrix @@ -115,7 +86,7 @@ module stdlib_linalg_pseudoinverse end subroutine stdlib_linalg_pseudoinvert_${ri}$ ! Function interface - function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva) + module function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva) !> Input matrix a[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] .... @@ -134,7 +105,7 @@ module stdlib_linalg_pseudoinverse end function stdlib_linalg_pseudoinverse_${ri}$ ! Inverse matrix operator - function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) + module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) !> Input matrix a[m,n] ${rt}$, intent(in), target :: a(:,:) !> Result matrix @@ -150,4 +121,4 @@ module stdlib_linalg_pseudoinverse #:endfor -end module stdlib_linalg_pseudoinverse +end submodule stdlib_linalg_pseudoinverse diff --git a/test/linalg/test_linalg_pseudoinverse.fypp b/test/linalg/test_linalg_pseudoinverse.fypp index e9524e39e..c0bdf883a 100644 --- a/test/linalg/test_linalg_pseudoinverse.fypp +++ b/test/linalg/test_linalg_pseudoinverse.fypp @@ -5,7 +5,6 @@ module test_linalg_pseudoinverse use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg use stdlib_linalg_constants - use stdlib_linalg_pseudoinverse implicit none (type,external) private From 73f9742d298cd2d859cb0c863077df2cb1f6f1a1 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:12:23 +0100 Subject: [PATCH 05/18] document interfaces --- src/stdlib_linalg.fypp | 115 +++++++++++++++++++++++++++++++++++------ 1 file changed, 98 insertions(+), 17 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index e6d415ae5..33b973b09 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -812,50 +812,131 @@ module stdlib_linalg end interface operator(.inv.) - ! Moose-Penrose Pseudo-Inverse: Function interface + ! Moore-Penrose Pseudo-Inverse: Function interface interface pinv + !! version: experimental + !! + !! Pseudo-inverse of a matrix + !! ([Specification](../page/specs/stdlib_linalg.html#pinv-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] .... + !> [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 + !> [optional] State return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err - !> Matrix pseudo-inverse + !> 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 - ! Matrix Inverse: Subroutine interface - in-place inversion + ! Moore-Penrose Pseudo-Inverse: Subroutine interface interface pseudoinvert + !! version: experimental + !! + !! Computation of the Moore-Penrose pseudo-inverse + !! ([Specification](../page/specs/stdlib_linalg.html#pseudoinvert-computation-of-a-matrix-pseudo-inverse)) + !! + !!### 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 - ${rt}$, intent(inout) :: pinva(:,:) - !> [optional] .... + ${rt}$, intent(in) :: 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 + !> [optional] State return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err - end subroutine + end subroutine stdlib_linalg_pseudoinvert_${ri}$ #:endfor end interface pseudoinvert - ! Operator interface + ! Moore-Penrose Pseudo-Inverse: Operator interface interface operator(.pinv.) - #:for rk,rt,ri in RC_KINDS_TYPES - module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) + !! version: experimental + !! + !! Pseudo-inverse operator of a matrix + !! ([Specification](../page/specs/stdlib_linalg.html#pinv-pseudo-inverse-operator-of-a-matrix)) + !! + !!### 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 matrix + !> Result pseudo-inverse matrix ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) - end function - #:endfor + 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 From 2e4f1a083f3c8619873247a98a4b0a800fc14a97 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:12:37 +0100 Subject: [PATCH 06/18] sync operator behavior with `.inv.` (return NaNs on error) --- src/stdlib_linalg_pinv.fypp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index 16bdcb68d..ce201b762 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -110,12 +110,23 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse ${rt}$, intent(in), target :: a(:,:) !> Result 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) + 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 From 61bf91c330951d0ee87658d93cc3fa5e448f8874 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:29:56 +0100 Subject: [PATCH 07/18] write specs --- doc/specs/stdlib_linalg.md | 136 ++++++++++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 6 +- src/stdlib_linalg_pinv.fypp | 10 +-- 3 files changed, 144 insertions(+), 8 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 2f03f0fad..d507c006e 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 `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!} +``` + +## `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/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 33b973b09..bb3be85b1 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -817,7 +817,7 @@ module stdlib_linalg !! version: experimental !! !! Pseudo-inverse of a matrix - !! ([Specification](../page/specs/stdlib_linalg.html#pinv-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. @@ -862,7 +862,7 @@ module stdlib_linalg !! version: experimental !! !! Computation of the Moore-Penrose pseudo-inverse - !! ([Specification](../page/specs/stdlib_linalg.html#pseudoinvert-computation-of-a-matrix-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 @@ -894,7 +894,7 @@ module stdlib_linalg #: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(in) :: a(:,:) + ${rt}$, intent(inout) :: a(:,:) !> Output pseudo-inverse matrix [n,m] ${rt}$, intent(out) :: pinva(:,:) !> [optional] Relative tolerance for singular value cutoff diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index ce201b762..570ce3f72 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -7,7 +7,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse use stdlib_linalg_lapack use stdlib_linalg_state use stdlib_linalg, only: svd - use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + use ieee_arithmetic, only: ieee_value, ieee_quiet_nan implicit none(type,external) character(*), parameter :: this = 'pseudo-inverse' @@ -21,8 +21,8 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse !> Input matrix a[m,n] ${rt}$, intent(inout) :: a(:,:) !> Output pseudo-inverse matrix - ${rt}$, intent(inout) :: pinva(:,:) - !> [optional] .... + ${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 @@ -89,7 +89,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse module function stdlib_linalg_pseudoinverse_${ri}$(a,rtol,err) result(pinva) !> Input matrix a[m,n] ${rt}$, intent(in), target :: a(:,:) - !> [optional] .... + !> [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 @@ -108,7 +108,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse module function stdlib_linalg_pinv_${ri}$_operator(a) result(pinva) !> Input matrix a[m,n] ${rt}$, intent(in), target :: a(:,:) - !> Result matrix + !> Result pseudo-inverse matrix ${rt}$ :: pinva(size(a,2,kind=ilp),size(a,1,kind=ilp)) type(linalg_state_type) :: err From 59a881d1218908cf4ac95e3683f61f893e58d674 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:36:29 +0100 Subject: [PATCH 08/18] add example --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_pseudoinverse.f90 | 32 ++++++++++++++++++++++++ 2 files changed, 33 insertions(+) create mode 100644 example/linalg/example_pseudoinverse.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index fc3f823d9..de9eb4f99 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 From 523bf4ccc5bab625be689a28ad6be54ab40adf35 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:37:18 +0100 Subject: [PATCH 09/18] fix spec link --- src/stdlib_linalg.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index bb3be85b1..4a01aa723 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -910,7 +910,7 @@ module stdlib_linalg !! version: experimental !! !! Pseudo-inverse operator of a matrix - !! ([Specification](../page/specs/stdlib_linalg.html#pinv-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. From 098d9b8b7425e28fcc85da87cca0a4c9cabde377 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:53:34 +0100 Subject: [PATCH 10/18] relax test tolerance (gcc-12 failure) --- test/linalg/test_linalg_pseudoinverse.fypp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/linalg/test_linalg_pseudoinverse.fypp b/test/linalg/test_linalg_pseudoinverse.fypp index c0bdf883a..263ccf454 100644 --- a/test/linalg/test_linalg_pseudoinverse.fypp +++ b/test/linalg/test_linalg_pseudoinverse.fypp @@ -41,7 +41,7 @@ module test_linalg_pseudoinverse integer(ilp) :: i,j integer(ilp), parameter :: n = 15_ilp - real(${rk}$), parameter :: tol = sqrt(epsilon(0.0_${rk}$)) + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) ${rt}$ :: a(n,n),inva(n,n) @@ -85,7 +85,7 @@ module test_linalg_pseudoinverse integer(ilp) :: failed integer(ilp), parameter :: n = 10 - real(${rk}$), parameter :: tol = sqrt(epsilon(0.0_${rk}$)) + 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) @@ -119,7 +119,7 @@ module test_linalg_pseudoinverse integer(ilp) :: failed integer(ilp), parameter :: m = 20, n = 10 - real(${rk}$), parameter :: tol = sqrt(epsilon(0.0_${rk}$)) + 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) @@ -153,7 +153,7 @@ module test_linalg_pseudoinverse integer(ilp) :: failed integer(ilp), parameter :: m = 10, n = 20 - real(${rk}$), parameter :: tol = sqrt(epsilon(0.0_${rk}$)) + 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) @@ -187,7 +187,7 @@ module test_linalg_pseudoinverse integer(ilp) :: failed integer(ilp), parameter :: n = 10 - real(${rk}$), parameter :: tol = sqrt(epsilon(0.0_${rk}$)) + 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) From 76368a0f4d2d12e69f2b429be6980de6d995aa48 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 27 Nov 2024 17:59:32 +0100 Subject: [PATCH 11/18] do not `use stdlib_linalg` --- src/stdlib_linalg_pinv.fypp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index 570ce3f72..91002ff2f 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -6,7 +6,6 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse use stdlib_linalg_blas use stdlib_linalg_lapack use stdlib_linalg_state - use stdlib_linalg, only: svd use ieee_arithmetic, only: ieee_value, ieee_quiet_nan implicit none(type,external) From 3d10d399c8ed8945438c0847f3a978cfea36ad12 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 18 Dec 2024 10:57:50 +0100 Subject: [PATCH 12/18] Update src/stdlib_linalg_pinv.fypp --- src/stdlib_linalg_pinv.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index 91002ff2f..76e956b6e 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -67,7 +67,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse !> Discard singular values cutoff = tolerance*maxval(s) - s = merge(1/s,0.0_${rk}$,s>cutoff) + 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 From 38bae4da2ca0b434d64bf614341780a0477e132d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 18 Dec 2024 10:58:01 +0100 Subject: [PATCH 13/18] Update src/stdlib_linalg_pinv.fypp --- src/stdlib_linalg_pinv.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index 76e956b6e..8bf30fd82 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -72,7 +72,7 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse ! 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 - forall (i=1:m,j=1:k) u(i,j) = s(j)*u(i,j) + do concurrent (i=1:m,j=1:k); u(i,j) = s(j)*u(i,j); end do ! 2) commutate matmul: A_pinv = V^H * (U * diag(1/s))^H = ((U * diag(1/s)) * V^H)^H. ! This avoids one matrix transpose From bf02cd775d4f1e88c6a56cafe3884f2ffd14b945 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 20 Dec 2024 16:38:56 +0100 Subject: [PATCH 14/18] use `gemm` instead of `matmul` --- src/stdlib_linalg_pinv.fypp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/stdlib_linalg_pinv.fypp b/src/stdlib_linalg_pinv.fypp index 8bf30fd82..e1d41f0a1 100644 --- a/src/stdlib_linalg_pinv.fypp +++ b/src/stdlib_linalg_pinv.fypp @@ -32,6 +32,8 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse ${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) @@ -74,13 +76,9 @@ submodule(stdlib_linalg) stdlib_linalg_pseudoinverse ! 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^H * (U * diag(1/s))^H = ((U * diag(1/s)) * V^H)^H. + ! 2) commutate matmul: A_pinv = V * (U * diag(1/s))^H = ((U * diag(1/s)) * V^H)^H. ! This avoids one matrix transpose - #:if rt.startswith('complex') - pinva = conjg(transpose(matmul(u,vt))) - #:else - pinva = transpose(matmul(u,vt)) - #:endif + 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}$ From 16a394aaa3d6a8b9fcf8b16244f3cc7872aca12a Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 21 Dec 2024 12:37:54 +0100 Subject: [PATCH 15/18] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 11e47235e..3b8e040a6 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1480,7 +1480,7 @@ singular value, are treated as zero. ### Syntax -`b ` [[stdlib_linalg(module):pinv(interface)]] `(a, [, rtol, err])` +`b =` [[stdlib_linalg(module):pinv(interface)]] `(a, [, rtol, err])` ### Arguments From 8912cf991b3c21d4ea626ea052d9507901d9f6fa Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sat, 21 Dec 2024 12:38:13 +0100 Subject: [PATCH 16/18] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 3b8e040a6..e23bbe2a3 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1459,7 +1459,7 @@ If `err` is not present, exceptions trigger an `error stop`. {!example/linalg/example_inverse_function.f90!} ``` -## `pinv` - Moore-Penrose pseudo-inverse of a matrix. +## `pinv` - Moore-Penrose pseudo-inverse of a matrix ### Status From 9dc5cc960293e7e20ac76c6f63c6f9c97f99fe30 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 26 Dec 2024 09:47:35 +0100 Subject: [PATCH 17/18] fix comment --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index e23bbe2a3..a7c8c0ed3 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1498,7 +1498,7 @@ where \( \epsilon \) is the machine precision for the element type of `a`. It is 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 `pinva` and `a` have degenerate or incompatible sizes. +Raises `LINALG_VALUE_ERROR` if `a` has invalid size. If `err` is not present, exceptions trigger an `error stop`. ### Example From 82977634ccab08f34700dab5db93b073e59e6109 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 26 Dec 2024 09:58:34 +0100 Subject: [PATCH 18/18] `loadtxt/savetxt`: export array size in error message --- src/stdlib_io.fypp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) 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