From f39bf52951332eef05d11c5183e99be1c6a6a925 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 15:22:13 +0200 Subject: [PATCH 01/29] add source file --- src/CMakeLists.txt | 1 + src/stdlib_linalg_inverse.fypp | 146 +++++++++++++++++++++++++++++++++ 2 files changed, 147 insertions(+) create mode 100644 src/stdlib_linalg_inverse.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a031ab823..09b410a15 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,6 +27,7 @@ set(fppFiles stdlib_linalg_cross_product.fypp stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp + stdlib_linalg_inverse.fypp stdlib_linalg_state.fypp stdlib_optval.fypp stdlib_selection.fypp diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp new file mode 100644 index 000000000..30d9ac9ec --- /dev/null +++ b/src/stdlib_linalg_inverse.fypp @@ -0,0 +1,146 @@ +#:include "common.fypp" +! Compute matrix inverse +module stdlib_linalg_inverse + use stdlib_linalg_constants + use stdlib_linalg_blas + use stdlib_linalg_lapack + use stdlib_linalg_state + use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + implicit none(type,external) + private + + !> Function interface return the matrix inverse + public :: inv + !> Subroutine interface: invert matrix inplace + public :: invert + !> Operator interface: .inv.A returns the matrix inverse of A + public :: operator(.inv.) + + ! Numpy: inv(a) + ! Scipy: inv(a, overwrite_a=False, check_finite=True) + ! IMSL: .i.a + + ! Function interface + interface inv + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_inverse_${ri}$ + #:endfor + end interface inv + + ! Subroutine interface: in-place factorization + interface invert + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_invert_${ri}$ + #:endfor + end interface invert + + ! Operator interface + interface operator(.inv.) + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_inverse_${ri}$_operator + #:endfor + end interface operator(.inv.) + + contains + + #:for rk,rt,ri in ALL_KINDS_TYPES + + ! Compute the in-place square matrix inverse of a + subroutine stdlib_linalg_invert_${ri}$(a,err) + !> Input matrix a[n,n] + ${rt}$, intent(inout) :: a(:,:) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), optional, intent(out) :: err + + !> Local variables + type(linalg_state) :: err0 + integer(ilp) :: lda,n,info,nb,lwork + integer(ilp), allocatable :: ipiv(:) + ${rt}$, allocatable :: work(:) + character(*), parameter :: this = 'invert' + + !> Problem sizes + lda = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + + if (lda<1 .or. n<1 .or. lda/=n) then + err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',lda,',',n,']') + goto 1 + end if + + ! Pivot indices + allocate(ipiv(n)) + + ! Factorize matrix (overwrite result) + call getrf(lda,n,a,lda,ipiv,info) + + ! Return codes from getrf and getri are identical + if (info==0) then + + ! Get optimal worksize (returned in work(1)) (apply 2% safety parameter) + nb = stdlib_ilaenv(1,'${ri}$getri',' ',n,-1,-1,-1) + lwork = nint(1.02*n*nb,kind=ilp) + + allocate(work(lwork)) + + ! Invert matrix + call getri(n,a,lda,ipiv,work,lwork,info) + + endif + + select case (info) + case (0) + ! Success + case (:-1) + err0 = linalg_state(this,LINALG_ERROR,'invalid matrix size a=[',lda,',',n,']') + case (1:) + ! Matrix is singular + err0 = linalg_state(this,LINALG_ERROR,'singular matrix') + case default + err0 = linalg_state(this,LINALG_INTERNAL_ERROR,'catastrophic error') + end select + + ! Process output and return + 1 call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_invert_${ri}$ + + ! Invert matrix in place + 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), optional, intent(out) :: err + + !> Allocate with copy + allocate(inva,source=a) + + !> Compute matrix inverse + call stdlib_linalg_invert_${ri}$(inva,err) + + end function stdlib_linalg_inverse_${ri}$ + + ! Inverse matrix operator + 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) :: err + + inva = stdlib_linalg_inverse_${ri}$(a,err) + + ! On error, return an empty matrix + if (err%error()) then + if (allocated(inva)) deallocate(inva) + allocate(inva(0,0)) + endif + + end function stdlib_linalg_inverse_${ri}$_operator + + #:endfor + +end module stdlib_linalg_inverse From c0da712ec2c9e0d2d6443a0a7585facba328585d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 15:33:01 +0200 Subject: [PATCH 02/29] working implementation --- src/stdlib_linalg_inverse.fypp | 59 ++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 30d9ac9ec..760d7f528 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -1,14 +1,16 @@ #:include "common.fypp" -! Compute matrix inverse +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_linalg_inverse +!! Compute inverse of a square matrix use stdlib_linalg_constants - use stdlib_linalg_blas - use stdlib_linalg_lapack - use stdlib_linalg_state - use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + 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 implicit none(type,external) private + character(*), parameter :: this = 'inverse' + !> Function interface return the matrix inverse public :: inv !> Subroutine interface: invert matrix inplace @@ -22,49 +24,54 @@ module stdlib_linalg_inverse ! Function interface interface inv - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" module procedure stdlib_linalg_inverse_${ri}$ + #:endif #:endfor end interface inv ! Subroutine interface: in-place factorization interface invert - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" module procedure stdlib_linalg_invert_${ri}$ + #:endif #:endfor end interface invert ! Operator interface interface operator(.inv.) - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" module procedure stdlib_linalg_inverse_${ri}$_operator + #:endif #:endfor end interface operator(.inv.) contains - #:for rk,rt,ri in ALL_KINDS_TYPES - + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" ! Compute the in-place square matrix inverse of a subroutine stdlib_linalg_invert_${ri}$(a,err) !> Input matrix a[n,n] - ${rt}$, intent(inout) :: a(:,:) + ${rt}$, intent(inout) :: a(:,:) !> [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 - type(linalg_state) :: err0 + type(linalg_state_type) :: err0 integer(ilp) :: lda,n,info,nb,lwork integer(ilp), allocatable :: ipiv(:) ${rt}$, allocatable :: work(:) - character(*), parameter :: this = 'invert' !> Problem sizes lda = size(a,1,kind=ilp) n = size(a,2,kind=ilp) if (lda<1 .or. n<1 .or. lda/=n) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',lda,',',n,']') + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',lda,',',n,']') goto 1 end if @@ -77,9 +84,9 @@ module stdlib_linalg_inverse ! Return codes from getrf and getri are identical if (info==0) then - ! Get optimal worksize (returned in work(1)) (apply 2% safety parameter) + ! Get optimal worksize (returned in work(1)) (inflate by a 2% safety margin) nb = stdlib_ilaenv(1,'${ri}$getri',' ',n,-1,-1,-1) - lwork = nint(1.02*n*nb,kind=ilp) + lwork = ceiling(1.02*n*nb,kind=ilp) allocate(work(lwork)) @@ -92,12 +99,12 @@ module stdlib_linalg_inverse case (0) ! Success case (:-1) - err0 = linalg_state(this,LINALG_ERROR,'invalid matrix size a=[',lda,',',n,']') + err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',lda,',',n,']') case (1:) ! Matrix is singular - err0 = linalg_state(this,LINALG_ERROR,'singular matrix') + err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix') case default - err0 = linalg_state(this,LINALG_INTERNAL_ERROR,'catastrophic error') + err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') end select ! Process output and return @@ -112,7 +119,7 @@ module stdlib_linalg_inverse !> Output matrix inverse ${rt}$, allocatable :: inva(:,:) !> [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 !> Allocate with copy allocate(inva,source=a) @@ -129,18 +136,16 @@ module stdlib_linalg_inverse !> Result matrix ${rt}$, allocatable :: inva(:,:) - type(linalg_state) :: err + type(linalg_state_type) :: err0 - inva = stdlib_linalg_inverse_${ri}$(a,err) + inva = stdlib_linalg_inverse_${ri}$(a,err0) ! On error, return an empty matrix - if (err%error()) then - if (allocated(inva)) deallocate(inva) - allocate(inva(0,0)) - endif + if (err0%error()) inva = 0.0_${rk}$ end function stdlib_linalg_inverse_${ri}$_operator + #:endif #:endfor end module stdlib_linalg_inverse From 92b1ac5987ac2ca4bdf0d7c4e03f175e593e8a43 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 15:43:57 +0200 Subject: [PATCH 03/29] reorganize as `submodule` --- src/stdlib_linalg.fypp | 47 +++++++++++++++++ src/stdlib_linalg_inverse.fypp | 92 ++++++++++++---------------------- 2 files changed, 79 insertions(+), 60 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 134bf7af5..d69d97868 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -20,6 +20,9 @@ module stdlib_linalg public :: operator(.det.) public :: diag public :: eye + public :: inv + public :: invert + public :: operator(.inv.) public :: lstsq public :: lstsq_space public :: solve @@ -552,6 +555,50 @@ module stdlib_linalg #:endfor end interface + ! Function interface + interface inv + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module function stdlib_linalg_inverse_${ri}$(a,err) result(inva) + !> Input matrix a[n,n] + ${rt}$, intent(in) :: a(:,:) + !> Output matrix inverse + ${rt}$, allocatable :: inva(:,:) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end function stdlib_linalg_inverse_${ri}$ + #:endif + #:endfor + end interface inv + + ! Subroutine interface: in-place factorization + interface invert + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module subroutine stdlib_linalg_invert_${ri}$(a,err) + !> Input matrix a[n,n] + ${rt}$, intent(inout) :: a(:,:) + !> [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_${ri}$ + #:endif + #:endfor + end interface invert + + ! Operator interface + interface operator(.inv.) + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva) + !> Input matrix a[n,n] + ${rt}$, intent(in) :: a(:,:) + !> Result matrix + ${rt}$, allocatable :: inva(:,:) + end function stdlib_linalg_inverse_${ri}$_operator + #:endif + #:endfor + end interface operator(.inv.) + contains diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 760d7f528..5bae6da8a 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -1,60 +1,40 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -module stdlib_linalg_inverse +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 implicit none(type,external) - private character(*), parameter :: this = 'inverse' - !> Function interface return the matrix inverse - public :: inv - !> Subroutine interface: invert matrix inplace - public :: invert - !> Operator interface: .inv.A returns the matrix inverse of A - public :: operator(.inv.) - - ! Numpy: inv(a) - ! Scipy: inv(a, overwrite_a=False, check_finite=True) - ! IMSL: .i.a - - ! Function interface - interface inv - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_inverse_${ri}$ - #:endif - #:endfor - end interface inv - - ! Subroutine interface: in-place factorization - interface invert - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_invert_${ri}$ - #:endif - #:endfor - end interface invert - - ! Operator interface - interface operator(.inv.) - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_inverse_${ri}$_operator - #:endif - #:endfor - end interface operator(.inv.) - 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 #:if rk!="xdp" ! Compute the in-place square matrix inverse of a - subroutine stdlib_linalg_invert_${ri}$(a,err) + module subroutine stdlib_linalg_invert_${ri}$(a,err) !> Input matrix a[n,n] ${rt}$, intent(inout) :: a(:,:) !> [optional] state return flag. On error if not requested, the code will stop @@ -72,7 +52,8 @@ module stdlib_linalg_inverse if (lda<1 .or. n<1 .or. lda/=n) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',lda,',',n,']') - goto 1 + call linalg_error_handling(err0,err) + return end if ! Pivot indices @@ -84,9 +65,9 @@ module stdlib_linalg_inverse ! Return codes from getrf and getri are identical if (info==0) then - ! Get optimal worksize (returned in work(1)) (inflate by a 2% safety margin) + ! Get optimal worksize (returned in work(1)) (inflate by a 5% safety margin) nb = stdlib_ilaenv(1,'${ri}$getri',' ',n,-1,-1,-1) - lwork = ceiling(1.02*n*nb,kind=ilp) + lwork = ceiling(1.05*n*nb,kind=ilp) allocate(work(lwork)) @@ -95,25 +76,16 @@ module stdlib_linalg_inverse endif - select case (info) - case (0) - ! Success - case (:-1) - err0 = linalg_state_type(this,LINALG_ERROR,'invalid matrix size a=[',lda,',',n,']') - case (1:) - ! Matrix is singular - err0 = linalg_state_type(this,LINALG_ERROR,'singular matrix') - case default - err0 = linalg_state_type(this,LINALG_INTERNAL_ERROR,'catastrophic error') - end select + ! Process output + call handle_getri_info(info,lda,n,err0) ! Process output and return - 1 call linalg_error_handling(err0,err) + call linalg_error_handling(err0,err) end subroutine stdlib_linalg_invert_${ri}$ ! Invert matrix in place - function stdlib_linalg_inverse_${ri}$(a,err) result(inva) + module function stdlib_linalg_inverse_${ri}$(a,err) result(inva) !> Input matrix a[n,n] ${rt}$, intent(in) :: a(:,:) !> Output matrix inverse @@ -130,7 +102,7 @@ module stdlib_linalg_inverse end function stdlib_linalg_inverse_${ri}$ ! Inverse matrix operator - function stdlib_linalg_inverse_${ri}$_operator(a) result(inva) + module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva) !> Input matrix a[n,n] ${rt}$, intent(in) :: a(:,:) !> Result matrix @@ -140,7 +112,7 @@ module stdlib_linalg_inverse inva = stdlib_linalg_inverse_${ri}$(a,err0) - ! On error, return an empty matrix + ! On error, return zeros if (err0%error()) inva = 0.0_${rk}$ end function stdlib_linalg_inverse_${ri}$_operator @@ -148,4 +120,4 @@ module stdlib_linalg_inverse #:endif #:endfor -end module stdlib_linalg_inverse +end submodule stdlib_linalg_inverse From e9a3f5bc384e5ddac41f01f3833a002815a65257 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 15:55:10 +0200 Subject: [PATCH 04/29] add tests --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_inverse.fypp | 164 +++++++++++++++++++++++++++ 2 files changed, 166 insertions(+) create mode 100644 test/linalg/test_linalg_inverse.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 9ec242213..934a3052c 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -3,6 +3,7 @@ set( "test_linalg.fypp" "test_blas_lapack.fypp" "test_linalg_solve.fypp" + "test_linalg_inverse.fypp" "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" "test_linalg_matrix_property_checks.fypp" @@ -12,6 +13,7 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_determinant) ADDTEST(linalg_matrix_property_checks) +ADDTEST(linalg_inverse) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(blas_lapack) diff --git a/test/linalg/test_linalg_inverse.fypp b/test/linalg/test_linalg_inverse.fypp new file mode 100644 index 000000000..5cbbd5d86 --- /dev/null +++ b/test/linalg/test_linalg_inverse.fypp @@ -0,0 +1,164 @@ +#: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.) + use stdlib_linalg_state, only: linalg_state_type + + 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 REAL_KINDS_TYPES + #:if rk!="xdp" + tests = [tests,new_unittest("inverse_${ri}$_eye_inverse",test_${ri}$_eye_inverse)] + #:endif + #:endfor + + end subroutine test_inverse_matrix + + !> Invert real identity matrix + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + subroutine test_${ri}$_eye_inverse(error) + type(error_type), allocatable, intent(out) :: error + + type(linalg_state_type) :: state + + integer(ilp), parameter :: n = 25_ilp + integer(ilp) :: i,j + ${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 + + !> 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 + call invert(a,err=state) + + call check(error,state%ok(),'inverse_${ri}$_eye (subroutine): '//state%print()) + if (allocated(error)) return + + call check(error,all(abs(a-inva) Invert complex identity matrix + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + 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)) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_inv From b233156ecd7d207d53b0bd2cd88b0d6c9bda4de9 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 16:06:37 +0200 Subject: [PATCH 05/29] option for pre-allocated pivot indices --- src/stdlib_linalg.fypp | 4 +++- src/stdlib_linalg_inverse.fypp | 30 ++++++++++++++++++++---------- 2 files changed, 23 insertions(+), 11 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index d69d97868..23dd0bd17 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -575,9 +575,11 @@ module stdlib_linalg interface invert #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" - module subroutine stdlib_linalg_invert_${ri}$(a,err) + module subroutine stdlib_linalg_invert_${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_${ri}$ diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 5bae6da8a..7e75469e2 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -34,31 +34,40 @@ submodule (stdlib_linalg) stdlib_linalg_inverse #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" ! Compute the in-place square matrix inverse of a - module subroutine stdlib_linalg_invert_${ri}$(a,err) - !> Input matrix a[n,n] + module subroutine stdlib_linalg_invert_${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 - integer(ilp), allocatable :: ipiv(:) + 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) - if (lda<1 .or. n<1 .or. lda/=n) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',lda,',',n,']') + ! 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 Date: Thu, 23 May 2024 16:45:56 +0200 Subject: [PATCH 07/29] document `.inv.A` --- doc/specs/stdlib_linalg.md | 36 +++++++++++++++++++++++++++++++++++- src/stdlib_linalg.fypp | 16 ++++++++++++++++ 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ffc4f7c42..ccfbd70b6 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -887,7 +887,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. @@ -898,3 +898,37 @@ Exceptions trigger an `error stop`. ```fortran {!example/linalg/example_determinant2.f90!} ``` + +## `.inv.` - Inverse operator of a square matrix + +### Status + +Experimental + +### Description + +This operator returns the inverse of a `real` or `complex` square matrix \( A \). +The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). + +This interface is equivalent to the function version of [[stdlib_linalg(module):inv(interface)]]. + +### Syntax + +`b = ` [[stdlib_linalg(module):operator(.inv.)(interface)]] `(a)` + +### Arguments + +`a`: Shall be a rank-2 square array of any `real` or `complex` kinds. It is an `intent(in)` argument. + +### Return value + +Returns a rank-2 square array with the same type, kind and rank as `a`, that contains the inverse of `a`. + +If the matrix is singular or non-square, the operator returns zeros. + +### Example + +```fortran +{!example/linalg/example_inverse1.f90!} +``` + diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 23dd0bd17..763a11cc2 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -589,6 +589,22 @@ module stdlib_linalg ! Operator interface interface operator(.inv.) + !! version: experimental + !! + !! Inverse operator of a square matrix + !! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-operator-of-a-square-matrix)) + !! + !!### Summary + !! Operator interface for computing the inverse of a square `real` or `complex` matrix. + !! + !!### Description + !! + !! This operator interface provides a convenient way to compute the inverse of a matrix. + !! Supported data types include `real` and `complex`. + !! + !!@note The provided functions are intended for square matrices. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva) From bf0dfd33730454ec194e18f1a555b54387ee524f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 16:56:29 +0200 Subject: [PATCH 08/29] document `invert` --- doc/specs/stdlib_linalg.md | 39 ++++++++++++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 20 +++++++++++++++++++ 2 files changed, 59 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ccfbd70b6..cdbbb002a 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -932,3 +932,42 @@ If the matrix is singular or non-square, the operator returns zeros. {!example/linalg/example_inverse1.f90!} ``` +## `invert` - Inversion of a square matrix. + +### Status + +Experimental + +### Description + +This subroutine inverts a square `real` or `complex` matrix in-place. +The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). + +On return, the input matrix `a` is replaced by its inverse. +The solver is based on LAPACK's `*GETRF` and `*GETRI` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):invert(interface)]] `(a, [, pivot] [, err])` + +### Arguments + +`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. +On output, it is replaced by the inverse of `a`. + +`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Replaces matrix \( A \) with its inverse, \(A^{-1}\). + +Raises `LINALG_ERROR` if the matrix is singular or has invalid size. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_inverse3.f90!} +``` diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 763a11cc2..e923a4ae7 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -573,6 +573,26 @@ module stdlib_linalg ! Subroutine interface: in-place factorization interface invert + !! version: experimental + !! + !! Inversion of a square matrix + !! ([Specification](../page/specs/stdlib_linalg.html#invert-inversion-of-a-square-matrix)) + !! + !!### Summary + !! This interface provides methods for inverting a square `real` or `complex` matrix in-place. + !! The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). + !! + !!### Description + !! + !! This operator interface provides a convenient way to compute the inverse of a matrix. + !! Supported data types include `real` and `complex`. + !! On completion, matrix `a` is replaced by the inverse matrix. Pre-allocated storage may be provided + !! for the array of pivot indices, `pivot`. If all pre-allocated work spaces are provided, no internal + !! memory allocations take place when using this interface. + !! + !!@note The provided subroutines are intended for square matrices. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module subroutine stdlib_linalg_invert_${ri}$(a,pivot,err) From db90b377748ef06aabd9c7e769441f604178e8c0 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 17:07:19 +0200 Subject: [PATCH 09/29] document `inv` --- doc/specs/stdlib_linalg.md | 36 ++++++++++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 22 +++++++++++++++++++++- 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index cdbbb002a..070f1c2b0 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -971,3 +971,39 @@ If `err` is not present, exceptions trigger an `error stop`. ```fortran {!example/linalg/example_inverse3.f90!} ``` + +## `inv` - Inverse of a square matrix. + +### Status + +Experimental + +### Description + +This function returns the inverse of a square `real` or `complex` matrix in-place. +The inverse, \( A^{-1} \), is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). + +The solver is based on LAPACK's `*GETRF` and `*GETRI` backends. + +### Syntax + +`b ` [[stdlib_linalg(module):inv(interface)]] `(a, [, err])` + +### Arguments + +`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns an array value of the same type, kind and rank as `a`, that contains the inverse matrix \(A^{-1}\). + +Raises `LINALG_ERROR` if the matrix is singular or has invalid size. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_inverse2.f90!} +``` diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index e923a4ae7..bdc4a82ea 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -557,6 +557,26 @@ module stdlib_linalg ! Function interface interface inv + !! version: experimental + !! + !! Inverse of a square matrix + !! ([Specification](../page/specs/stdlib_linalg.html#inv-inverse-of-a-square-matrix)) + !! + !!### Summary + !! This interface provides methods for computing the inverse of a square `real` or `complex` matrix. + !! The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). + !! + !!### Description + !! + !! This function interface provides methods that return the inverse of a square matrix. + !! Supported data types include `real` and `complex`. + !! The inverse matrix \( A^{-1} \) is returned as a function result. + !! Exceptions are raised in case of singular matrix or invalid size, and trigger an `error stop` if + !! the state flag `err` is not provided. + !! + !!@note The provided functions are intended for square matrices. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module function stdlib_linalg_inverse_${ri}$(a,err) result(inva) @@ -584,7 +604,7 @@ module stdlib_linalg !! !!### Description !! - !! This operator interface provides a convenient way to compute the inverse of a matrix. + !! This subroutine interface provides a way to compute the inverse of a matrix in-place. !! Supported data types include `real` and `complex`. !! On completion, matrix `a` is replaced by the inverse matrix. Pre-allocated storage may be provided !! for the array of pivot indices, `pivot`. If all pre-allocated work spaces are provided, no internal From 6a1c397f7b52c2f6cd51f6232da7e56efb0801fb Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 23 May 2024 17:49:44 +0200 Subject: [PATCH 10/29] typo --- doc/specs/stdlib_linalg.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 070f1c2b0..08a2a68fe 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -879,7 +879,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 @@ -914,7 +914,7 @@ This interface is equivalent to the function version of [[stdlib_linalg(module) ### Syntax -`b = ` [[stdlib_linalg(module):operator(.inv.)(interface)]] `(a)` +`b = ` [[stdlib_linalg(module):operator(.inv.)(interface)]] `a` ### Arguments From e10e4c439bb7230adb371a244d1db95248bd348c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Mon, 27 May 2024 12:38:46 +0200 Subject: [PATCH 11/29] test for singular matrix; activate `complex` matrix tests --- test/linalg/test_linalg_inverse.fypp | 55 +++++++++++++++++++++++----- 1 file changed, 46 insertions(+), 9 deletions(-) diff --git a/test/linalg/test_linalg_inverse.fypp b/test/linalg/test_linalg_inverse.fypp index 5cbbd5d86..93752240d 100644 --- a/test/linalg/test_linalg_inverse.fypp +++ b/test/linalg/test_linalg_inverse.fypp @@ -4,8 +4,8 @@ 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.) - use stdlib_linalg_state, only: linalg_state_type + use stdlib_linalg, only: inv,invert,operator(.inv.),eye + use stdlib_linalg_state, only: linalg_state_type,LINALG_ERROR implicit none (type,external) private @@ -21,29 +21,27 @@ module test_linalg_inverse allocate(tests(0)) - #:for rk,rt,ri in REAL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" - tests = [tests,new_unittest("inverse_${ri}$_eye_inverse",test_${ri}$_eye_inverse)] + tests = [tests,new_unittest("${ri}$_eye_inverse",test_${ri}$_eye_inverse), & + new_unittest("${ri}$_singular_inverse",test_${ri}$_singular_inverse)] #:endif #:endfor end subroutine test_inverse_matrix - !> Invert real identity matrix #:for rk,rt,ri in REAL_KINDS_TYPES #:if rk!="xdp" + !> 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 - integer(ilp) :: i,j ${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 + a = eye(n) !> Inverse function inva = inv(a,err=state) @@ -64,6 +62,27 @@ module test_linalg_inverse end subroutine test_${ri}$_eye_inverse + !> 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 + #:endif #:endfor @@ -132,6 +151,24 @@ module test_linalg_inverse end subroutine test_${ci}$_eye_inverse + !> 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 + #:endif #:endfor From 5f320f9f3729ac171865dfb1b9a9790387392ee3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 09:45:36 +0200 Subject: [PATCH 12/29] subroutine version, non in-place --- doc/specs/stdlib_linalg.md | 11 ++++++--- src/stdlib_linalg.fypp | 25 +++++++++++++++----- src/stdlib_linalg_inverse.fypp | 42 +++++++++++++++++++++++++++++++--- 3 files changed, 66 insertions(+), 12 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 08a2a68fe..d8ce700f8 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -948,12 +948,16 @@ The solver is based on LAPACK's `*GETRF` and `*GETRI` backends. ### Syntax -`call ` [[stdlib_linalg(module):invert(interface)]] `(a, [, pivot] [, err])` +`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. It is an `intent(inout)` argument. -On output, it is replaced by the inverse of `a`. +`a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. +If `inva` is provided, it is an `intent(in)` argument. +If `inva` is not provided, it is an `intent(inout)` argument: on output, it is replaced by the inverse of `a`. + +`inva` (optional): Shall be a rank-2, square, `real` or `complex` array with the same size, and kind as `a`. +On output, it contains the inverse of `a`. `pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices. @@ -964,6 +968,7 @@ On output, it is replaced by the inverse of `a`. Replaces matrix \( A \) with its inverse, \(A^{-1}\). Raises `LINALG_ERROR` if the matrix is singular or has invalid size. +Raises `LINALG_VALUE_ERROR` if `inva` and `a` do not have the same size. If `err` is not present, exceptions trigger an `error stop`. ### Example diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index bdc4a82ea..e010e0d6d 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -604,25 +604,38 @@ module stdlib_linalg !! !!### Description !! - !! This subroutine interface provides a way to compute the inverse of a matrix in-place. + !! This subroutine interface provides a way to compute the inverse of a matrix. !! Supported data types include `real` and `complex`. - !! On completion, matrix `a` is replaced by the inverse matrix. 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. + !! The user may provide a unique matrix argument `a`. In this case, `a` is replaced by the inverse matrix. + !! on output. Otherwise, one may provide two separate arguments: an input matrix `a` and an output matrix + !! `inva`. In this case, `a` will not be modified, and the inverse is returned in `inva`. + !! Pre-allocated storage may be provided for the array of pivot indices, `pivot`. If all pre-allocated + !! work spaces are provided, no internal memory allocations take place when using this interface. !! !!@note The provided subroutines are intended for square matrices. !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" - module subroutine stdlib_linalg_invert_${ri}$(a,pivot,err) + 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_${ri}$ + end subroutine stdlib_linalg_invert_inplace_${ri}$ + ! Compute the square matrix inverse of a + module subroutine stdlib_linalg_invert_split_${ri}$(a,inva,pivot,err) + !> Input matrix a[n,n]. + ${rt}$, intent(in) :: a(:,:) + !> Inverse matrix a[n,n]. + ${rt}$, intent(out) :: inva(:,:) + !> [optional] Storage array for the diagonal pivot indices + integer(ilp), optional, intent(inout), target :: pivot(:) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_invert_split_${ri}$ #:endif #:endfor end interface invert diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 7e75469e2..2f72f8b3a 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -34,7 +34,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" ! Compute the in-place square matrix inverse of a - module subroutine stdlib_linalg_invert_${ri}$(a,pivot,err) + 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 @@ -92,7 +92,43 @@ submodule (stdlib_linalg) stdlib_linalg_inverse if (.not.present(pivot)) deallocate(ipiv) call linalg_error_handling(err0,err) - end subroutine stdlib_linalg_invert_${ri}$ + 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 + + 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) @@ -107,7 +143,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse allocate(inva,source=a) !> Compute matrix inverse - call stdlib_linalg_invert_${ri}$(inva,err=err) + call stdlib_linalg_invert_inplace_${ri}$(inva,err=err) end function stdlib_linalg_inverse_${ri}$ From 47ee61a0a54ee88fb127ee68a9a245cc85dbac40 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 09:45:44 +0200 Subject: [PATCH 13/29] test subroutine version, non in-place --- test/linalg/test_linalg_inverse.fypp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/test/linalg/test_linalg_inverse.fypp b/test/linalg/test_linalg_inverse.fypp index 93752240d..aa5e955a0 100644 --- a/test/linalg/test_linalg_inverse.fypp +++ b/test/linalg/test_linalg_inverse.fypp @@ -49,15 +49,24 @@ module test_linalg_inverse if (allocated(error)) return call check(error,all(abs(a-inva) Inverse subroutine - call invert(a,err=state) + if (allocated(error)) return + + !> 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) Date: Wed, 5 Jun 2024 09:47:58 +0200 Subject: [PATCH 14/29] `.inv.` operator: `error stop` on issues --- doc/specs/stdlib_linalg.md | 2 +- src/stdlib_linalg_inverse.fypp | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index d8ce700f8..43acfa513 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -924,7 +924,7 @@ This interface is equivalent to the function version of [[stdlib_linalg(module) Returns a rank-2 square array with the same type, kind and rank as `a`, that contains the inverse of `a`. -If the matrix is singular or non-square, the operator returns zeros. +Exceptions always trigger an `error stop`. ### Example diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 2f72f8b3a..846c872ce 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -154,12 +154,8 @@ submodule (stdlib_linalg) stdlib_linalg_inverse !> Result matrix ${rt}$, allocatable :: inva(:,:) - type(linalg_state_type) :: err0 - - inva = stdlib_linalg_inverse_${ri}$(a,err0) - - ! On error, return zeros - if (err0%error()) inva = 0.0_${rk}$ + ! Do not provide an error handler (error stop on issues) + inva = stdlib_linalg_inverse_${ri}$(a) end function stdlib_linalg_inverse_${ri}$_operator From 8b138f714761455dd7c787b5867b579dbb60e58f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:35:25 +0200 Subject: [PATCH 15/29] add non-in-place subroutine example --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_inverse4.f90 | 22 ++++++++++++++++++++++ 2 files changed, 23 insertions(+) create mode 100644 example/linalg/example_inverse4.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 29e91117b..50c58dd5f 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -15,6 +15,7 @@ ADD_EXAMPLE(is_triangular) ADD_EXAMPLE(inverse1) ADD_EXAMPLE(inverse2) ADD_EXAMPLE(inverse3) +ADD_EXAMPLE(inverse4) ADD_EXAMPLE(outer_product) ADD_EXAMPLE(trace) ADD_EXAMPLE(state1) diff --git a/example/linalg/example_inverse4.f90 b/example/linalg/example_inverse4.f90 new file mode 100644 index 000000000..5713589db --- /dev/null +++ b/example/linalg/example_inverse4.f90 @@ -0,0 +1,22 @@ +! Matrix inversion example: subroutine interface +program example_inverse4 + 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_inverse4 From 513d77bc351d22bd7e676f6963c9a3302a26f9e3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:38:31 +0200 Subject: [PATCH 16/29] rename example programs --- doc/specs/stdlib_linalg.md | 10 +++++++--- example/linalg/CMakeLists.txt | 8 ++++---- ...ample_inverse2.f90 => example_inverse_function.f90} | 4 ++-- ...xample_inverse3.f90 => example_inverse_inplace.f90} | 4 ++-- ...ample_inverse1.f90 => example_inverse_operator.f90} | 4 ++-- ...ple_inverse4.f90 => example_inverse_subroutine.f90} | 4 ++-- 6 files changed, 19 insertions(+), 15 deletions(-) rename example/linalg/{example_inverse2.f90 => example_inverse_function.f90} (88%) rename example/linalg/{example_inverse3.f90 => example_inverse_inplace.f90} (88%) rename example/linalg/{example_inverse1.f90 => example_inverse_operator.f90} (88%) rename example/linalg/{example_inverse4.f90 => example_inverse_subroutine.f90} (87%) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 2d59f3861..ee77dd38d 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1020,7 +1020,7 @@ Exceptions always trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_inverse1.f90!} +{!example/linalg/example_inverse_operator.f90!} ``` ## `invert` - Inversion of a square matrix. @@ -1065,7 +1065,11 @@ If `err` is not present, exceptions trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_inverse3.f90!} +{!example/linalg/example_inverse_inplace.f90!} +``` + +```fortran +{!example/linalg/example_inverse_subroutine.f90!} ``` ## `inv` - Inverse of a square matrix. @@ -1101,5 +1105,5 @@ If `err` is not present, exceptions trigger an `error stop`. ### Example ```fortran -{!example/linalg/example_inverse2.f90!} +{!example/linalg/example_inverse_function.f90!} ``` diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 50c58dd5f..f54e06183 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -12,10 +12,10 @@ ADD_EXAMPLE(is_skew_symmetric) ADD_EXAMPLE(is_square) ADD_EXAMPLE(is_symmetric) ADD_EXAMPLE(is_triangular) -ADD_EXAMPLE(inverse1) -ADD_EXAMPLE(inverse2) -ADD_EXAMPLE(inverse3) -ADD_EXAMPLE(inverse4) +ADD_EXAMPLE(inverse_operator) +ADD_EXAMPLE(inverse_function) +ADD_EXAMPLE(inverse_inplace) +ADD_EXAMPLE(inverse_subroutine) ADD_EXAMPLE(outer_product) ADD_EXAMPLE(trace) ADD_EXAMPLE(state1) diff --git a/example/linalg/example_inverse2.f90 b/example/linalg/example_inverse_function.f90 similarity index 88% rename from example/linalg/example_inverse2.f90 rename to example/linalg/example_inverse_function.f90 index fa5aba8e7..a59b0ebdb 100644 --- a/example/linalg/example_inverse2.f90 +++ b/example/linalg/example_inverse_function.f90 @@ -1,5 +1,5 @@ ! Matrix inversion example: function interface -program example_inverse2 +program example_inverse_function use stdlib_linalg_constants, only: dp use stdlib_linalg, only: inv,eye implicit none @@ -19,4 +19,4 @@ program example_inverse2 ! Final check print *, 'CHECK passed? ',matmul(A,Am1)==eye(2) -end program example_inverse2 +end program example_inverse_function diff --git a/example/linalg/example_inverse3.f90 b/example/linalg/example_inverse_inplace.f90 similarity index 88% rename from example/linalg/example_inverse3.f90 rename to example/linalg/example_inverse_inplace.f90 index e06760738..9a3816a34 100644 --- a/example/linalg/example_inverse3.f90 +++ b/example/linalg/example_inverse_inplace.f90 @@ -1,5 +1,5 @@ ! Matrix inversion example: in-place inversion -program example_inverse3 +program example_inverse_inplace use stdlib_linalg_constants, only: dp use stdlib_linalg, only: invert,eye implicit none @@ -20,4 +20,4 @@ program example_inverse3 ! Final check print *, 'CHECK passed? ',matmul(A,Am1)==eye(2) -end program example_inverse3 +end program example_inverse_inplace diff --git a/example/linalg/example_inverse1.f90 b/example/linalg/example_inverse_operator.f90 similarity index 88% rename from example/linalg/example_inverse1.f90 rename to example/linalg/example_inverse_operator.f90 index a27009718..a4f71c75d 100644 --- a/example/linalg/example_inverse1.f90 +++ b/example/linalg/example_inverse_operator.f90 @@ -1,5 +1,5 @@ ! Matrix inversion example: operator interface -program example_inverse1 +program example_inverse_operator use stdlib_linalg_constants, only: dp use stdlib_linalg, only: operator(.inv.),eye implicit none @@ -19,4 +19,4 @@ program example_inverse1 ! Final check print *, 'CHECK passed? ',matmul(A,Am1)==eye(2) -end program example_inverse1 +end program example_inverse_operator diff --git a/example/linalg/example_inverse4.f90 b/example/linalg/example_inverse_subroutine.f90 similarity index 87% rename from example/linalg/example_inverse4.f90 rename to example/linalg/example_inverse_subroutine.f90 index 5713589db..1c75d04f0 100644 --- a/example/linalg/example_inverse4.f90 +++ b/example/linalg/example_inverse_subroutine.f90 @@ -1,5 +1,5 @@ ! Matrix inversion example: subroutine interface -program example_inverse4 +program example_inverse_subroutine use stdlib_linalg_constants, only: dp use stdlib_linalg, only: invert,eye implicit none @@ -19,4 +19,4 @@ program example_inverse4 ! Final check print *, 'CHECK passed? ',matmul(A,Am1)==eye(2) -end program example_inverse4 +end program example_inverse_subroutine From 5d7d7ba2ae8562f27d009f7660bf0b374e9cb719 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:40:50 +0200 Subject: [PATCH 17/29] 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 2d59f3861..a67e16cad 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1001,7 +1001,7 @@ Experimental This operator returns the inverse of a `real` or `complex` square matrix \( A \). The inverse \( A^{-1} \) is defined such that \( A \cdot A^{-1} = A^{-1} \cdot A = I_n \). -This interface is equivalent to the function version of [[stdlib_linalg(module):inv(interface)]]. +This interface is equivalent to the function [[stdlib_linalg(module):inv(interface)]]. ### Syntax From 63006d2d540a93390a47c65146386c629d528f3d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:41:03 +0200 Subject: [PATCH 18/29] 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 a67e16cad..72c111183 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1023,7 +1023,7 @@ Exceptions always trigger an `error stop`. {!example/linalg/example_inverse1.f90!} ``` -## `invert` - Inversion of a square matrix. +## `invert` - Inversion of a square matrix ### Status From 41f502e65866e7ac2a44f7568bc6c80677b0a7cd Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:41:16 +0200 Subject: [PATCH 19/29] 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 72c111183..3e7e1a389 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1050,7 +1050,7 @@ If `inva` is not provided, it is an `intent(inout)` argument: on output, it is r `inva` (optional): Shall be a rank-2, square, `real` or `complex` array with the same size, and kind as `a`. On output, it contains the inverse of `a`. -`pivot` (optional): Shall be a rank-1 array of the same kind and matrix dimension as `a`, providing storage for the diagonal pivot indices. It is an `intent(inout)` arguments, and returns the diagonal pivot indices. +`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. From 860cbb02f079ea29edf1c5de511c828a3398f9b2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:41:52 +0200 Subject: [PATCH 20/29] 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 3e7e1a389..f31197a6e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1056,7 +1056,7 @@ On output, it contains the inverse of `a`. ### Return value -Replaces matrix \( A \) with its inverse, \(A^{-1}\). +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. From 80a7742bafe842b66766ab77b7c2c6d6b1199661 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 15:22:24 +0200 Subject: [PATCH 21/29] new test: random SPD matrix (real, complex) --- test/linalg/test_linalg_inverse.fypp | 119 ++++++++++++++++++++++++++- 1 file changed, 118 insertions(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_inverse.fypp b/test/linalg/test_linalg_inverse.fypp index aa5e955a0..885b4c552 100644 --- a/test/linalg/test_linalg_inverse.fypp +++ b/test/linalg/test_linalg_inverse.fypp @@ -24,7 +24,8 @@ module test_linalg_inverse #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" tests = [tests,new_unittest("${ri}$_eye_inverse",test_${ri}$_eye_inverse), & - new_unittest("${ri}$_singular_inverse",test_${ri}$_singular_inverse)] + new_unittest("${ri}$_singular_inverse",test_${ri}$_singular_inverse), & + new_unittest("${ri}$_random_spd_inverse",test_${ri}$_random_spd_inverse)] #:endif #:endfor @@ -91,7 +92,53 @@ module test_linalg_inverse 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)) 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 From 3c0bcdbffba80e9249bbc485e7e22756eaa5e8d7 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 15:22:49 +0200 Subject: [PATCH 22/29] 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 f31197a6e..dd4aad5cb 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1089,7 +1089,7 @@ The solver is based on LAPACK's `*GETRF` and `*GETRI` backends. `a`: Shall be a rank-2, square, `real` or `complex` array containing the coefficient matrix. It is an `intent(inout)` argument. -`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. +`err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. ### Return value From 505f30a2a38d571429439d36e51670be8df4baf2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 15:24:45 +0200 Subject: [PATCH 23/29] lwork: overflow-safer evaluation --- src/stdlib_linalg_inverse.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 846c872ce..01a8925b1 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -76,7 +76,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse ! Get optimal worksize (returned in work(1)) (inflate by a 5% safety margin) nb = stdlib_ilaenv(1,'${ri}$getri',' ',n,-1,-1,-1) - lwork = ceiling(1.05*n*nb,kind=ilp) + lwork = max(1,min(huge(0_ilp),ceiling(1.05_${rk}$*real(n,${rk}$)*nb,kind=ilp))) allocate(work(lwork)) From d0af9bee5d4ebbf8a8e5fefb7e0ce0f261de23ed Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 15:31:10 +0200 Subject: [PATCH 24/29] `operator(.inv.)`: return `NaN`s on exceptions --- doc/specs/stdlib_linalg.md | 3 ++- src/stdlib_linalg.fypp | 3 ++- src/stdlib_linalg_inverse.fypp | 13 +++++++++++-- 3 files changed, 15 insertions(+), 4 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ee77dd38d..e8c0ed233 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1015,7 +1015,8 @@ This interface is equivalent to the function version of [[stdlib_linalg(module) Returns a rank-2 square array with the same type, kind and rank as `a`, that contains the inverse of `a`. -Exceptions always trigger an `error stop`. +If an exception occurred on input errors, or singular matrix, NaNs will be returned. + ### Example diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 6d2f1356c..a4e813978 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -655,7 +655,8 @@ module stdlib_linalg !!### Description !! !! This operator interface provides a convenient way to compute the inverse of a matrix. - !! Supported data types include `real` and `complex`. + !! Supported data types include `real` and `complex`. On input errors or singular matrix, + !! NaNs will be returned. !! !!@note The provided functions are intended for square matrices. !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index 01a8925b1..f018cd035 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -6,6 +6,7 @@ submodule (stdlib_linalg) stdlib_linalg_inverse 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' @@ -153,9 +154,17 @@ submodule (stdlib_linalg) stdlib_linalg_inverse ${rt}$, intent(in) :: a(:,:) !> Result matrix ${rt}$, allocatable :: inva(:,:) + + type(linalg_state_type) :: err - ! Do not provide an error handler (error stop on issues) - inva = stdlib_linalg_inverse_${ri}$(a) + ! 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)),source=ieee_value(1.0_${rk}$,ieee_quiet_nan)) + endif end function stdlib_linalg_inverse_${ri}$_operator From 406e17034513c4b95683af6844bcc5e1c5bab36d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 15:51:46 +0200 Subject: [PATCH 25/29] typo: set complex NaN --- src/stdlib_linalg_inverse.fypp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index f018cd035..b200700e6 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -163,7 +163,14 @@ submodule (stdlib_linalg) stdlib_linalg_inverse ! 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)),source=ieee_value(1.0_${rk}$,ieee_quiet_nan)) + 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 From 3e7c82eee4924b9a51cb06df4ecb8e599bb551fe Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 16:03:13 +0200 Subject: [PATCH 26/29] lstsq: explain singular values better --- doc/specs/stdlib_linalg.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ec6575292..4461b8d03 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -777,8 +777,7 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \( `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument. -`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `m -al(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` argument. +`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. From a292e9df496d5b048e249b898c25b7d3513d1cd4 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 27 Jun 2024 08:47:18 +0200 Subject: [PATCH 27/29] Update stdlib_linalg.md --- doc/specs/stdlib_linalg.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 4461b8d03..13432619c 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1014,7 +1014,9 @@ This interface is equivalent to the function [[stdlib_linalg(module):inv(interf Returns a rank-2 square array with the same type, kind and rank as `a`, that contains the inverse of `a`. -If an exception occurred on input errors, or singular matrix, NaNs will be returned. +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 From 7386d2d8a59ca76dcd1c92bc573082f17cacd050 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 2 Jul 2024 15:08:38 +0200 Subject: [PATCH 28/29] Update src/stdlib_linalg.fypp Co-authored-by: Jeremie Vandenplas --- 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 97cb8d4a2..3d9741971 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -597,7 +597,7 @@ module stdlib_linalg #:endfor end interface inv - ! Matrix Inverse: Subroutine interface - in-place factorization + ! Matrix Inverse: Subroutine interface - in-place inversion interface invert !! version: experimental !! From bb3f7e4c5614ae58baa9e3b29b5029eb6023cb63 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 5 Jul 2024 12:22:51 +0200 Subject: [PATCH 29/29] remove `xdp` restriction --- doc/specs/stdlib_linalg.md | 3 ++- src/stdlib_linalg.fypp | 9 --------- src/stdlib_linalg_inverse.fypp | 2 -- test/linalg/test_linalg_inverse.fypp | 6 ------ 4 files changed, 2 insertions(+), 18 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 02e4f302a..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 diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 28b2cf7c0..9cb7a13a7 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -567,10 +567,8 @@ module stdlib_linalg !! the state flag `err` is not provided. !! !!@note The provided functions are intended for square matrices. - !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). !! #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" module function stdlib_linalg_inverse_${ri}$(a,err) result(inva) !> Input matrix a[n,n] ${rt}$, intent(in) :: a(:,:) @@ -579,7 +577,6 @@ module stdlib_linalg !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err end function stdlib_linalg_inverse_${ri}$ - #:endif #:endfor end interface inv @@ -605,10 +602,8 @@ module stdlib_linalg !! work spaces are provided, no internal memory allocations take place when using this interface. !! !!@note The provided subroutines are intended for square matrices. - !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). !! #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" module subroutine stdlib_linalg_invert_inplace_${ri}$(a,pivot,err) !> Input matrix a[n,n] ${rt}$, intent(inout) :: a(:,:) @@ -628,7 +623,6 @@ module stdlib_linalg !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err end subroutine stdlib_linalg_invert_split_${ri}$ - #:endif #:endfor end interface invert @@ -649,17 +643,14 @@ module stdlib_linalg !! NaNs will be returned. !! !!@note The provided functions are intended for square matrices. - !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). !! #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" module function stdlib_linalg_inverse_${ri}$_operator(a) result(inva) !> Input matrix a[n,n] ${rt}$, intent(in) :: a(:,:) !> Result matrix ${rt}$, allocatable :: inva(:,:) end function stdlib_linalg_inverse_${ri}$_operator - #:endif #:endfor end interface operator(.inv.) diff --git a/src/stdlib_linalg_inverse.fypp b/src/stdlib_linalg_inverse.fypp index b200700e6..38ad3dda7 100644 --- a/src/stdlib_linalg_inverse.fypp +++ b/src/stdlib_linalg_inverse.fypp @@ -33,7 +33,6 @@ submodule (stdlib_linalg) stdlib_linalg_inverse end subroutine handle_getri_info #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" ! 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 @@ -175,7 +174,6 @@ submodule (stdlib_linalg) stdlib_linalg_inverse end function stdlib_linalg_inverse_${ri}$_operator - #:endif #:endfor end submodule stdlib_linalg_inverse diff --git a/test/linalg/test_linalg_inverse.fypp b/test/linalg/test_linalg_inverse.fypp index 885b4c552..679a586c5 100644 --- a/test/linalg/test_linalg_inverse.fypp +++ b/test/linalg/test_linalg_inverse.fypp @@ -22,17 +22,14 @@ module test_linalg_inverse allocate(tests(0)) #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" 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)] - #:endif #:endfor end subroutine test_inverse_matrix #:for rk,rt,ri in REAL_KINDS_TYPES - #:if rk!="xdp" !> Invert real identity matrix subroutine test_${ri}$_eye_inverse(error) type(error_type), allocatable, intent(out) :: error @@ -139,12 +136,10 @@ module test_linalg_inverse end subroutine test_${ri}$_random_spd_inverse - #:endif #:endfor !> Invert complex identity matrix #:for ck,ct,ci in CMPLX_KINDS_TYPES - #:if ck!="xdp" subroutine test_${ci}$_eye_inverse(error) type(error_type), allocatable, intent(out) :: error @@ -295,7 +290,6 @@ module test_linalg_inverse end subroutine test_${ci}$_singular_inverse - #:endif #:endfor end module test_linalg_inverse