From 42aa3e5b0a69ceee1eb5cfaff09a25614b6623de Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 08:55:27 +0200 Subject: [PATCH 01/26] add `svd` --- src/CMakeLists.txt | 3 +- src/stdlib_linalg_svd.fypp | 276 +++++++++++++++++++++++++++++++++++++ 2 files changed, 278 insertions(+), 1 deletion(-) create mode 100644 src/stdlib_linalg_svd.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b58ba0b9c..9a23eae99 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -26,7 +26,8 @@ set(fppFiles stdlib_linalg_kronecker.fypp stdlib_linalg_cross_product.fypp stdlib_linalg_determinant.fypp - stdlib_linalg_state.fypp + stdlib_linalg_state.fypp + stdlib_linalg_svd.fypp stdlib_optval.fypp stdlib_selection.fypp stdlib_sorting.fypp diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp new file mode 100644 index 000000000..2d77e70f8 --- /dev/null +++ b/src/stdlib_linalg_svd.fypp @@ -0,0 +1,276 @@ +#:include "common.fypp" +module stdlib_linalg_svd + 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 + + !> Singular value decomposition + public :: svd + !> Singular values + public :: svdvals + + ! Numpy: svd(a, full_matrices=True, compute_uv=True, hermitian=False) + ! Scipy: svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd') + + interface svd + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_svd_${ri}$ + #:endfor + end interface svd + + interface svdvals + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_svdvals_${ri}$ + #:endfor + end interface svdvals + + !> Return full matrices U, V^T to separate storage + character, parameter :: GESDD_FULL_MATRICES = 'A' + + !> Return shrunk matrices U, V^T to k = min(m,n) + character, parameter :: GESDD_SHRINK_MATRICES = 'S' + + !> Overwrite A storage with U (if M>=N) or VT (if M Do not return either U or VT (singular values array only) + character, parameter :: GESDD_SINGVAL_ONLY = 'N' + + character(*), parameter :: this = 'svd' + + + contains + + !> Process GESDD output flag + elemental subroutine gesdd_info(err,info,m,n) + !> Error handler + type(linalg_state), intent(inout) :: err + !> GESDD return flag + integer(ilp), intent(in) :: info + !> Input matrix size + integer(ilp), intent(in) :: m,n + + select case (info) + case (0) + ! Success! + err%state = LINALG_SUCCESS + case (-1) + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid task ID on input to GESDD.') + case (-5,-3:-2) + err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',m,',',n,']') + case (-8) + err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix U size, with a=[',m,',',n,']') + case (-10) + err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix V size, with a=[',m,',',n,']') + case (-4) + err = linalg_state(this,LINALG_VALUE_ERROR,'A contains invalid/NaN values.') + case (1:) + err = linalg_state(this,LINALG_ERROR,'SVD computation did not converge.') + case default + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Unknown error returned by GESDD.') + end select + + end subroutine gesdd_info + + + #:for rk,rt,ri in ALL_KINDS_TYPES + + !> Singular values of matrix A + function stdlib_linalg_svdvals_${ri}$(a,err) result(s) + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), optional, intent(out) :: err + !> Array of singular values + real(${rk}$), allocatable :: s(:) + + !> Create + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(s(k)) + + !> Compute singular values + call stdlib_linalg_svd_${ri}$(amat,s,overwrite_a=.false.,err=err) + + end function stdlib_linalg_svdvals_${ri}$ + + !> SVD of matrix A = U S V^T, returning S and optionally U and V^T + subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err) + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of singular values + real(${rk}$), intent(out) :: s(:) + !> The columns of U contain the eigenvectors of A A^T + ${rt}$, optional, intent(out), target :: u(:,:) + !> The rows of V^T contain the eigenvectors of A^T A + ${rt}$, optional, intent(out), target :: vt(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] full matrices have shape(u)==[m,m], shape(vh)==[n,n] (default); otherwise + !> they are shape(u)==[m,k] and shape(vh)==[k,n] with k=min(m,n) + logical(lk), optional, intent(in) :: full_matrices + !> [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) :: m,n,lda,ldu,ldvt,info,k,lwork,liwork,lrwork + integer(ilp), allocatable :: iwork(:) + logical(lk) :: copy_a,full_storage,compute_uv,alloc_u,alloc_vt,can_overwrite_a + character :: task + ${rt}$, target :: work_dummy(1),u_dummy(1,1),vt_dummy(1,1) + ${rt}$, allocatable :: work(:) + real(${rk}$), allocatable :: rwork(:) + ${rt}$, pointer :: amat(:,:),umat(:,:),vtmat(:,:) + + !> Matrix determinant size + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + lda = m + + if (.not.k>0) then + err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or matrix size: a=[',m,',',n,']') + goto 1 + end if + + if (.not.size(s,kind=ilp)>=k) then + err0 = linalg_state(this,LINALG_VALUE_ERROR,'singular value array has insufficient size:',& + ' s=[',size(s,kind=ilp),'], k=',k) + goto 1 + endif + + ! Integer storage + liwork = 8*k + allocate(iwork(liwork)) + + ! Can A be overwritten? By default, do not overwrite + if (present(overwrite_a)) then + copy_a = .not.overwrite_a + else + copy_a = .true._lk + endif + + ! Initialize a matrix temporary + if (copy_a) then + allocate(amat(m,n),source=a) + + ! Check if we can overwrite A with data that will be lost + can_overwrite_a = merge(.not.present(u),.not.present(vt),m>=n) + + else + amat => a + + can_overwrite_a = .false. + + endif + + ! Full-size matrices + if (present(full_matrices)) then + full_storage = full_matrices + else + full_storage = .true. + endif + + ! Decide if U, VT matrices should be computed + compute_uv = present(u) .or. present(vt) + + ! U, VT storage + if (present(u)) then + umat => u + alloc_u = .false. + elseif ((copy_a .and. m>=n) .or. .not.compute_uv) then + ! U not wanted, and A can be overwritten: do not allocate + umat => u_dummy + alloc_u = .false. + elseif (.not.full_storage) then + allocate(umat(m,k)) + alloc_u = .true. + else + allocate(umat(m,m)) + alloc_u = .true. + end if + + if (present(vt)) then + vtmat => vt + alloc_vt = .false. + elseif ((copy_a .and. m vt_dummy + alloc_vt = .false. + elseif (.not.full_storage) then + allocate(vtmat(k,n)) + alloc_vt = .true. + else + allocate(vtmat(n,n)) + alloc_vt = .true. + end if + + ldu = size(umat ,1,kind=ilp) + ldvt = size(vtmat,1,kind=ilp) + + ! Decide SVD task + if (.not.compute_uv) then + task = GESDD_SINGVAL_ONLY + elseif (can_overwrite_a) then + ! A is a copy: we can overwrite its storage + task = GESDD_OVERWRITE_A + elseif (.not.full_storage) then + task = GESDD_SHRINK_MATRICES + else + task = GESDD_FULL_MATRICES + end if + + ! Compute workspace + #:if rt.startswith('complex') + if (task==GESDD_SINGVAL_ONLY) then + lrwork = max(1,7*k) + else + lrwork = max(1,5*k*(k+1),2*k*(k+max(m,n))+k) + endif + allocate(rwork(lrwork)) + #:endif + + lwork = -1_ilp + + call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,& + work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#iwork,info) + call gesdd_info(err0,info,m,n) + + ! Compute SVD + if (info==0) then + + !> Prepare working storage + lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) + allocate(work(lwork)) + + !> Compute SVD + call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,& + work,lwork,#{if rt.startswith('comp')}#rwork,#{endif}#iwork,info) + call gesdd_info(err0,info,m,n) + + endif + + ! Finalize storage and process output flag + if (copy_a) deallocate(amat) + if (alloc_u) deallocate(umat) + if (alloc_vt) deallocate(vtmat) + 1 call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_svd_${ri}$ + #:endfor + +end module stdlib_linalg_svd From e1d49d85a3cc61a55288a1928b1f3ed520b185f5 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 09:27:23 +0200 Subject: [PATCH 02/26] base implementation --- src/stdlib_linalg_svd.fypp | 72 ++++++++++++++++++++------------------ 1 file changed, 38 insertions(+), 34 deletions(-) diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index 2d77e70f8..c4a9d5056 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -1,12 +1,14 @@ #:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_linalg_svd +!! Singular-Value Decomposition 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: gesdd + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS implicit none(type,external) - private + + character(*), parameter :: this = 'svd' !> Singular value decomposition public :: svd @@ -17,14 +19,18 @@ module stdlib_linalg_svd ! Scipy: svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd') interface svd - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" module procedure stdlib_linalg_svd_${ri}$ + #:endif #:endfor end interface svd interface svdvals - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" module procedure stdlib_linalg_svdvals_${ri}$ + #:endif #:endfor end interface svdvals @@ -40,15 +46,15 @@ module stdlib_linalg_svd !> Do not return either U or VT (singular values array only) character, parameter :: GESDD_SINGVAL_ONLY = 'N' - character(*), parameter :: this = 'svd' + contains !> Process GESDD output flag - elemental subroutine gesdd_info(err,info,m,n) + elemental subroutine handle_gesdd_info(err,info,m,n) !> Error handler - type(linalg_state), intent(inout) :: err + type(linalg_state_type), intent(inout) :: err !> GESDD return flag integer(ilp), intent(in) :: info !> Input matrix size @@ -59,32 +65,33 @@ module stdlib_linalg_svd ! Success! err%state = LINALG_SUCCESS case (-1) - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid task ID on input to GESDD.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID on input to GESDD.') case (-5,-3:-2) - err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',m,',',n,']') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',m,',',n,']') case (-8) - err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix U size, with a=[',m,',',n,']') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix U size, with a=[',m,',',n,']') case (-10) - err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix V size, with a=[',m,',',n,']') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix V size, with a=[',m,',',n,']') case (-4) - err = linalg_state(this,LINALG_VALUE_ERROR,'A contains invalid/NaN values.') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'A contains invalid/NaN values.') case (1:) - err = linalg_state(this,LINALG_ERROR,'SVD computation did not converge.') + err = linalg_state_type(this,LINALG_ERROR,'SVD computation did not converge.') case default - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Unknown error returned by GESDD.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by GESDD.') end select - end subroutine gesdd_info + end subroutine handle_gesdd_info - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" !> Singular values of matrix A function stdlib_linalg_svdvals_${ri}$(a,err) result(s) !> Input matrix A[m,n] ${rt}$, intent(in), target :: 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 !> Array of singular values real(${rk}$), allocatable :: s(:) @@ -123,10 +130,10 @@ module stdlib_linalg_svd !> they are shape(u)==[m,k] and shape(vh)==[k,n] with k=min(m,n) logical(lk), optional, intent(in) :: full_matrices !> [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) :: m,n,lda,ldu,ldvt,info,k,lwork,liwork,lrwork integer(ilp), allocatable :: iwork(:) logical(lk) :: copy_a,full_storage,compute_uv,alloc_u,alloc_vt,can_overwrite_a @@ -143,12 +150,12 @@ module stdlib_linalg_svd lda = m if (.not.k>0) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or matrix size: a=[',m,',',n,']') + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size: a=[',m,',',n,']') goto 1 end if if (.not.size(s,kind=ilp)>=k) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'singular value array has insufficient size:',& + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'singular value array has insufficient size:',& ' s=[',size(s,kind=ilp),'], k=',k) goto 1 endif @@ -167,17 +174,13 @@ module stdlib_linalg_svd ! Initialize a matrix temporary if (copy_a) then allocate(amat(m,n),source=a) - - ! Check if we can overwrite A with data that will be lost - can_overwrite_a = merge(.not.present(u),.not.present(vt),m>=n) - else amat => a - - can_overwrite_a = .false. - endif - + + ! Check if we can overwrite A with data that will be lost + can_overwrite_a = copy_a .and. merge(.not.present(u),.not.present(vt),m>=n) + ! Full-size matrices if (present(full_matrices)) then full_storage = full_matrices @@ -248,7 +251,7 @@ module stdlib_linalg_svd call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,& work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#iwork,info) - call gesdd_info(err0,info,m,n) + call handle_gesdd_info(err0,info,m,n) ! Compute SVD if (info==0) then @@ -260,7 +263,7 @@ module stdlib_linalg_svd !> Compute SVD call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,& work,lwork,#{if rt.startswith('comp')}#rwork,#{endif}#iwork,info) - call gesdd_info(err0,info,m,n) + call handle_gesdd_info(err0,info,m,n) endif @@ -271,6 +274,7 @@ module stdlib_linalg_svd 1 call linalg_error_handling(err0,err) end subroutine stdlib_linalg_svd_${ri}$ + #:endif #:endfor end module stdlib_linalg_svd From f5adf8f95c8d57f182f888d930edbfc56ce4030d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 09:28:18 +0200 Subject: [PATCH 03/26] 1d matrix shape --- src/stdlib_linalg_svd.fypp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index c4a9d5056..41e8e297c 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -67,11 +67,11 @@ module stdlib_linalg_svd case (-1) err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID on input to GESDD.') case (-5,-3:-2) - err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=[',m,',',n,']') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) case (-8) - err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix U size, with a=[',m,',',n,']') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix U size, with a=',[m,n]) case (-10) - err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix V size, with a=[',m,',',n,']') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix V size, with a=',[m,n]) case (-4) err = linalg_state_type(this,LINALG_VALUE_ERROR,'A contains invalid/NaN values.') case (1:) @@ -150,7 +150,7 @@ module stdlib_linalg_svd lda = m if (.not.k>0) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size: a=[',m,',',n,']') + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size: a=',[m,n]) goto 1 end if From 64adda4471cc7b16d4fc2925a59fd726bbced914 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 09:42:16 +0200 Subject: [PATCH 04/26] cleanup --- src/stdlib_linalg_lapack.fypp | 18 +++++++++--------- src/stdlib_linalg_svd.fypp | 14 +++++++++++--- 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/src/stdlib_linalg_lapack.fypp b/src/stdlib_linalg_lapack.fypp index 8d0b4ce2c..a089d6b58 100644 --- a/src/stdlib_linalg_lapack.fypp +++ b/src/stdlib_linalg_lapack.fypp @@ -19,16 +19,16 @@ module stdlib_linalg_lapack interface bbcsd !! BBCSD computes the CS decomposition of a unitary matrix in !! bidiagonal-block form, - !! [ B11 | B12 0 0 ] - !! [ 0 | 0 -I 0 ] + !! [ B11 | B12 0 0 ] + !! [ 0 | 0 -I 0 ] !! X = [----------------] - !! [ B21 | B22 0 0 ] - !! [ 0 | 0 0 I ] - !! [ C | -S 0 0 ] - !! [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H - !! = [---------] [---------------] [---------] . - !! [ | U2 ] [ S | C 0 0 ] [ | V2 ] - !! [ 0 | 0 0 I ] + !! [ B21 | B22 0 0 ] + !! [ 0 | 0 0 I ] + !! [ C | -S 0 0 ] + !! [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H + !! = [---------] [---------------] [---------] . + !! [ | U2 ] [ S | C 0 0 ] [ | V2 ] + !! [ 0 | 0 0 I ] !! X is M-by-M, its top-left block is P-by-Q, and Q must be no larger !! than P, M-P, or M-Q. (If Q is not the smallest index, then X must be !! transposed and/or permuted. This can be done in constant time using diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index 41e8e297c..f082f418b 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -156,7 +156,7 @@ module stdlib_linalg_svd if (.not.size(s,kind=ilp)>=k) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'singular value array has insufficient size:',& - ' s=[',size(s,kind=ilp),'], k=',k) + ' s=',shape(s,kind=ilp),', k=',k) goto 1 endif @@ -193,6 +193,7 @@ module stdlib_linalg_svd ! U, VT storage if (present(u)) then + ! User input umat => u alloc_u = .false. elseif ((copy_a .and. m>=n) .or. .not.compute_uv) then @@ -200,14 +201,17 @@ module stdlib_linalg_svd umat => u_dummy alloc_u = .false. elseif (.not.full_storage) then + ! Allocate with minimum size allocate(umat(m,k)) alloc_u = .true. else + ! Allocate with regular size allocate(umat(m,m)) alloc_u = .true. end if if (present(vt)) then + ! User input vtmat => vt alloc_vt = .false. elseif ((copy_a .and. m vt_dummy alloc_vt = .false. elseif (.not.full_storage) then + ! Allocate with minimum size allocate(vtmat(k,n)) alloc_vt = .true. else + ! Allocate with regular size allocate(vtmat(n,n)) alloc_vt = .true. end if @@ -245,10 +251,12 @@ module stdlib_linalg_svd lrwork = max(1,5*k*(k+1),2*k*(k+max(m,n))+k) endif allocate(rwork(lrwork)) + #:else + lrwork = -1_ilp ! not needed #:endif + ! First call: request working storage space lwork = -1_ilp - call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,& work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#iwork,info) call handle_gesdd_info(err0,info,m,n) @@ -262,7 +270,7 @@ module stdlib_linalg_svd !> Compute SVD call gesdd(task,m,n,amat,lda,s,umat,ldu,vtmat,ldvt,& - work,lwork,#{if rt.startswith('comp')}#rwork,#{endif}#iwork,info) + work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#iwork,info) call handle_gesdd_info(err0,info,m,n) endif From 036c574cb24ffc66abe3433debc69e198c45d0cc Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 09:46:31 +0200 Subject: [PATCH 05/26] add test programs --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_svd.fypp | 201 +++++++++++++++++++++++++++++++ 2 files changed, 203 insertions(+) create mode 100644 test/linalg/test_linalg_svd.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 5a07b5bc4..e524886c3 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -4,6 +4,7 @@ set( "test_blas_lapack.fypp" "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" + "test_linalg_svd.fypp" "test_linalg_matrix_property_checks.fypp" ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -12,4 +13,5 @@ ADDTEST(linalg) ADDTEST(linalg_determinant) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_lstsq) +ADDTEST(linalg_svd) ADDTEST(blas_lapack) diff --git a/test/linalg/test_linalg_svd.fypp b/test/linalg/test_linalg_svd.fypp new file mode 100644 index 000000000..9c0d72009 --- /dev/null +++ b/test/linalg/test_linalg_svd.fypp @@ -0,0 +1,201 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test singular value decomposition +module test_linalg_svd + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg, only: diag + use stdlib_linalg_svd, only: svd,svdvals + use stdlib_linalg_state, only: linalg_state_type + + implicit none (type,external) + + contains + + !> SVD tests + subroutine test_svd(error) + logical, intent(out) :: error + + real :: t0,t1 + + call cpu_time(t0) + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + call test_svd_${ri}$(error) + if (error) return + #:endif + #:endfor + + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + call test_complex_svd_${ci}$(error) + if (error) return + #:endif + #:endfor + + call cpu_time(t1) + + print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error) + + 1 format('SVD tests completed in ',f9.4,' milliseconds, result=',a) + + end subroutine test_svd + + !> Real matrix svd + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + subroutine test_svd_${ri}$(error) + logical,intent(out) :: error + + !> Reference solution + ${rt}$, parameter :: tol = sqrt(epsilon(0.0_${rk}$)) + ${rt}$, parameter :: third = 1.0_${rk}$/3.0_${rk}$ + ${rt}$, parameter :: twothd = 2*third + ${rt}$, parameter :: rsqrt2 = 1.0_${rk}$/sqrt(2.0_${rk}$) + ${rt}$, parameter :: rsqrt18 = 1.0_${rk}$/sqrt(18.0_${rk}$) + + ${rt}$, parameter :: A_mat(2,3) = reshape([${rt}$ :: 3,2, 2,3, 2,-2],[2,3]) + ${rt}$, parameter :: s_sol(2) = [${rt}$ :: 5, 3] + ${rt}$, parameter :: u_sol(2,2) = reshape(rsqrt2*[1,1,1,-1],[2,2]) + ${rt}$, parameter :: vt_sol(3,3) = reshape([rsqrt2,rsqrt18,twothd, & + rsqrt2,-rsqrt18,-twothd,& + 0.0_${rk}$,4*rsqrt18,-third],[3,3]) + + !> Local variables + type(linalg_state_type) :: state + ${rt}$ :: A(2,3),s(2),u(2,2),vt(3,3) + + !> Initialize matrix + A = A_mat + + !> Simple subroutine version + call svd(A,s,err=state) + error = state%error() .or. .not. all(abs(s-s_sol)<=tol) + if (error) return + + !> Function interface + s = svdvals(A,err=state) + error = state%error() .or. .not. all(abs(s-s_sol)<=tol) + if (error) return + + !> [S, U]. Singular vectors could be all flipped + call svd(A,s,u,err=state) + error = state%error() .or. & + .not. all(abs(s-s_sol)<=tol) .or. & + .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) + if (error) return + + !> [S, U]. Overwrite A matrix + call svd(A,s,u,overwrite_a=.true.,err=state) + error = state%error() .or. & + .not. all(abs(s-s_sol)<=tol) .or. & + .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) + if (error) return + + !> [S, U, V^T] + A = A_mat + call svd(A,s,u,vt,overwrite_a=.true.,err=state) + error = state%error() .or. & + .not. all(abs(s-s_sol)<=tol) .or. & + .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) .or. & + .not.(all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) + if (error) return + + !> [S, V^T]. Do not overwrite A matrix + A = A_mat + call svd(A,s,vt=vt,err=state) + error = state%error() .or. & + .not. all(abs(s-s_sol)<=tol) .or. & + .not.(all(abs(vt+vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) + if (error) return + + !> [S, V^T]. Overwrite A matrix + call svd(A,s,vt=vt,overwrite_a=.true.,err=state) + error = state%error() .or. & + .not. all(abs(s-s_sol)<=tol) .or. & + .not.(all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) + if (error) return + + !> [U, S, V^T]. + A = A_mat + call svd(A,s,u,vt,err=state) + error = state%error() .or. & + .not. all(abs(s-s_sol)<=tol) .or. & + .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) .or. & + .not.(all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) + if (error) return + + !> [U, S, V^T]. Partial storage -> compare until k=2 columns of U rows of V^T + A = A_mat + u = 0 + vt = 0 + call svd(A,s,u,vt,full_matrices=.false.,err=state) + error = state%error() & + .or. .not. all(abs(s-s_sol)<=tol) & + .or. .not.(all(abs( u(:,:2)- u_sol(:,:2))<=tol) .or. all(abs( u(:,:2)+ u_sol(:,:2))<=tol)) & + .or. .not.(all(abs(vt(:2,:)-vt_sol(:2,:))<=tol) .or. all(abs(vt(:2,:)+vt_sol(:2,:))<=tol)) + + if (error) return + + end subroutine test_svd_${ri}$ + + #:endif + #:endfor + + !> Test complex svd + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + subroutine test_complex_svd_${ci}$(error) + logical,intent(out) :: error + + !> Reference solution + real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) + real(${ck}$), parameter :: one = 1.0_${ck}$ + real(${ck}$), parameter :: zero = 0.0_${ck}$ + real(${ck}$), parameter :: sqrt2 = sqrt(2.0_${ck}$) + real(${ck}$), parameter :: rsqrt2 = one/sqrt2 + ${ct}$, parameter :: cone = (1.0_${ck}$,0.0_${ck}$) + ${ct}$, parameter :: cimg = (0.0_${ck}$,1.0_${ck}$) + ${ct}$, parameter :: czero = (0.0_${ck}$,0.0_${ck}$) + + real(${ck}$), parameter :: s_sol(2) = [sqrt2,sqrt2] + ${ct}$, parameter :: A_mat(2,2) = reshape([cone,cimg,cimg,cone],[2,2]) + ${ct}$, parameter :: u_sol(2,2) = reshape(rsqrt2*[cone,cimg,cimg,cone],[2,2]) + ${ct}$, parameter :: vt_sol(2,2) = reshape([cone,czero,czero,cone],[2,2]) + + !> Local variables + type(linalg_state_type) :: state + ${ct}$ :: A(2,2),u(2,2),vt(2,2) + real(${ck}$) :: s(2) + + !> Initialize matrix + A = A_mat + + !> Simple subroutine version + call svd(A,s,err=state) + error = state%error() .or. .not. all(abs(s-s_sol)<=tol) + if (error) return + + !> Function interface + s = svdvals(A,err=state) + error = state%error() .or. .not. all(abs(s-s_sol)<=tol) + if (error) return + + !> [S, U, V^T] + A = A_mat + call svd(A,s,u,vt,overwrite_a=.true.,err=state) + error = state%error() .or. & + .not. all(abs(s-s_sol)<=tol) .or. & + .not. all(abs(matmul(u,matmul(diag(s),vt)) - A_mat)<=tol) + if (error) return + + end subroutine test_complex_svd_${ci}$ + + #:endif + #:endfor + + +end module test_linalg_svd + + From 4a90c0953b332ab8d4bd915308a1366c6837659c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 10:01:56 +0200 Subject: [PATCH 06/26] tests: replace with `testdrive` --- test/linalg/test_linalg_svd.fypp | 206 ++++++++++++++++++++----------- 1 file changed, 137 insertions(+), 69 deletions(-) diff --git a/test/linalg/test_linalg_svd.fypp b/test/linalg/test_linalg_svd.fypp index 9c0d72009..eb3ec7e29 100644 --- a/test/linalg/test_linalg_svd.fypp +++ b/test/linalg/test_linalg_svd.fypp @@ -9,44 +9,37 @@ module test_linalg_svd use stdlib_linalg_state, only: linalg_state_type implicit none (type,external) + + public :: test_svd contains - !> SVD tests - subroutine test_svd(error) - logical, intent(out) :: error - - real :: t0,t1 - - call cpu_time(t0) + !> Solve several SVD problems + subroutine test_svd(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" - call test_svd_${ri}$(error) - if (error) return + tests = [tests,new_unittest("test_svd_${ri}$",test_svd_${ri}$)] #:endif #:endfor #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if ck!="xdp" - call test_complex_svd_${ci}$(error) - if (error) return + tests = [tests,new_unittest("test_complex_svd_${ci}$",test_complex_svd_${ci}$)] #:endif #:endfor - call cpu_time(t1) - - print 1, 1000*(t1-t0), merge('SUCCESS','ERROR ',.not.error) - - 1 format('SVD tests completed in ',f9.4,' milliseconds, result=',a) - end subroutine test_svd !> Real matrix svd #:for rk,rt,ri in REAL_KINDS_TYPES #:if rk!="xdp" subroutine test_svd_${ri}$(error) - logical,intent(out) :: error + type(error_type), allocatable, intent(out) :: error !> Reference solution ${rt}$, parameter :: tol = sqrt(epsilon(0.0_${rk}$)) @@ -63,6 +56,7 @@ module test_linalg_svd 0.0_${rk}$,4*rsqrt18,-third],[3,3]) !> Local variables + character(:), allocatable :: test type(linalg_state_type) :: state ${rt}$ :: A(2,3),s(2),u(2,2),vt(3,3) @@ -71,72 +65,110 @@ module test_linalg_svd !> Simple subroutine version call svd(A,s,err=state) - error = state%error() .or. .not. all(abs(s-s_sol)<=tol) - if (error) return - + + test = 'subroutine version' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + !> Function interface s = svdvals(A,err=state) - error = state%error() .or. .not. all(abs(s-s_sol)<=tol) - if (error) return + + test = 'function interface' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return !> [S, U]. Singular vectors could be all flipped call svd(A,s,u,err=state) - error = state%error() .or. & - .not. all(abs(s-s_sol)<=tol) .or. & - .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) - if (error) return + + test = 'subroutine with singular vectors' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol), test//': U') + if (allocated(error)) return !> [S, U]. Overwrite A matrix call svd(A,s,u,overwrite_a=.true.,err=state) - error = state%error() .or. & - .not. all(abs(s-s_sol)<=tol) .or. & - .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) - if (error) return + + test = 'subroutine, overwrite_a' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol), test//': U') + if (allocated(error)) return !> [S, U, V^T] A = A_mat call svd(A,s,u,vt,overwrite_a=.true.,err=state) - error = state%error() .or. & - .not. all(abs(s-s_sol)<=tol) .or. & - .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) .or. & - .not.(all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) - if (error) return - + + test = '[S, U, V^T]' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol), test//': U') + if (allocated(error)) return + call check(error, all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol), test//': V^T') + if (allocated(error)) return + !> [S, V^T]. Do not overwrite A matrix A = A_mat call svd(A,s,vt=vt,err=state) - error = state%error() .or. & - .not. all(abs(s-s_sol)<=tol) .or. & - .not.(all(abs(vt+vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) - if (error) return + + test = '[S, V^T], overwrite_a=.false.' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol), test//': V^T') + if (allocated(error)) return !> [S, V^T]. Overwrite A matrix call svd(A,s,vt=vt,overwrite_a=.true.,err=state) - error = state%error() .or. & - .not. all(abs(s-s_sol)<=tol) .or. & - .not.(all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) - if (error) return - + + test = '[S, V^T], overwrite_a=.true.' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol), test//': V^T') + if (allocated(error)) return + !> [U, S, V^T]. A = A_mat call svd(A,s,u,vt,err=state) - error = state%error() .or. & - .not. all(abs(s-s_sol)<=tol) .or. & - .not.(all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol)) .or. & - .not.(all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol)) - if (error) return + + test = '[U, S, V^T]' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(u-u_sol)<=tol) .or. all(abs(u+u_sol)<=tol), test//': U') + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(vt-vt_sol)<=tol) .or. all(abs(vt+vt_sol)<=tol), test//': V^T') + if (allocated(error)) return !> [U, S, V^T]. Partial storage -> compare until k=2 columns of U rows of V^T A = A_mat u = 0 vt = 0 call svd(A,s,u,vt,full_matrices=.false.,err=state) - error = state%error() & - .or. .not. all(abs(s-s_sol)<=tol) & - .or. .not.(all(abs( u(:,:2)- u_sol(:,:2))<=tol) .or. all(abs( u(:,:2)+ u_sol(:,:2))<=tol)) & - .or. .not.(all(abs(vt(:2,:)-vt_sol(:2,:))<=tol) .or. all(abs(vt(:2,:)+vt_sol(:2,:))<=tol)) - - if (error) return + + test = '[U, S, V^T], partial storage' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(u(:,:2)-u_sol(:,:2))<=tol) .or. all(abs(u(:,:2)+u_sol(:,:2))<=tol), test//': U(:,:2)') + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(vt(:2,:)-vt_sol(:2,:))<=tol) .or. all(abs(vt(:2,:)+vt_sol(:2,:))<=tol), test//': V^T(:2,:)') + if (allocated(error)) return end subroutine test_svd_${ri}$ @@ -147,7 +179,7 @@ module test_linalg_svd #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if ck!="xdp" subroutine test_complex_svd_${ci}$(error) - logical,intent(out) :: error + type(error_type), allocatable, intent(out) :: error !> Reference solution real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) @@ -165,6 +197,7 @@ module test_linalg_svd ${ct}$, parameter :: vt_sol(2,2) = reshape([cone,czero,czero,cone],[2,2]) !> Local variables + character(:), allocatable :: test type(linalg_state_type) :: state ${ct}$ :: A(2,2),u(2,2),vt(2,2) real(${ck}$) :: s(2) @@ -174,28 +207,63 @@ module test_linalg_svd !> Simple subroutine version call svd(A,s,err=state) - error = state%error() .or. .not. all(abs(s-s_sol)<=tol) - if (error) return - + + test = '[S], complex subroutine' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + !> Function interface s = svdvals(A,err=state) - error = state%error() .or. .not. all(abs(s-s_sol)<=tol) - if (error) return + + test = 'svdvals, complex function' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return !> [S, U, V^T] A = A_mat call svd(A,s,u,vt,overwrite_a=.true.,err=state) - error = state%error() .or. & - .not. all(abs(s-s_sol)<=tol) .or. & - .not. all(abs(matmul(u,matmul(diag(s),vt)) - A_mat)<=tol) - if (error) return + + test = '[S, U, V^T], complex' + call check(error,state%ok(),test//': '//state%print()) + if (allocated(error)) return + call check(error, all(abs(s-s_sol)<=tol), test//': S') + if (allocated(error)) return + call check(error, all(abs(matmul(u,matmul(diag(s),vt))-A_mat)<=tol), test//': U*S*V^T') + if (allocated(error)) return end subroutine test_complex_svd_${ci}$ #:endif #:endfor - end module test_linalg_svd - +program test_lstsq + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_svd, only : test_svd + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_svd", test_svd) & + ] + + 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_lstsq From 063b4216774e54260680ee60307f115ab642b933 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 10:14:38 +0200 Subject: [PATCH 07/26] create `submodule` --- src/stdlib_linalg.fypp | 54 ++++++++++++++++++++++++++++++++ src/stdlib_linalg_svd.fypp | 37 ++++------------------ test/linalg/test_linalg_svd.fypp | 4 +-- 3 files changed, 61 insertions(+), 34 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 460aa0593..bc9ce5950 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -24,6 +24,8 @@ module stdlib_linalg public :: lstsq_space public :: solve_lstsq public :: trace + public :: svd + public :: svdvals public :: outer_product public :: kronecker_product public :: cross_product @@ -456,6 +458,58 @@ module stdlib_linalg #:endfor end interface + ! Singular value decomposition + interface svd + !! version: experimental + !! + !!### Description + !! + !! + + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err) + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of singular values + real(${rk}$), intent(out) :: s(:) + !> The columns of U contain the eigenvectors of A A^T + ${rt}$, optional, intent(out), target :: u(:,:) + !> The rows of V^T contain the eigenvectors of A^T A + ${rt}$, optional, intent(out), target :: vt(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] full matrices have shape(u)==[m,m], shape(vh)==[n,n] (default); otherwise + !> they are shape(u)==[m,k] and shape(vh)==[k,n] with k=min(m,n) + logical(lk), optional, intent(in) :: full_matrices + !> [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_svd_${ri}$ + #:endif + #:endfor + end interface svd + + ! Singular values + interface svdvals + !! version: experimental + !! + !!### Description + !! + !! + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + module function stdlib_linalg_svdvals_${ri}$(a,err) result(s) + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + !> Array of singular values + real(${rk}$), allocatable :: s(:) + end function stdlib_linalg_svdvals_${ri}$ + #:endif + #:endfor + end interface svdvals + contains diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index f082f418b..8bc72ffdc 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -1,6 +1,6 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -module stdlib_linalg_svd +submodule(stdlib_linalg) stdlib_linalg_svd !! Singular-Value Decomposition use stdlib_linalg_constants use stdlib_linalg_lapack, only: gesdd @@ -9,30 +9,8 @@ module stdlib_linalg_svd implicit none(type,external) character(*), parameter :: this = 'svd' - - !> Singular value decomposition - public :: svd - !> Singular values - public :: svdvals - - ! Numpy: svd(a, full_matrices=True, compute_uv=True, hermitian=False) - ! Scipy: svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd') - - interface svd - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_svd_${ri}$ - #:endif - #:endfor - end interface svd - - interface svdvals - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_svdvals_${ri}$ - #:endif - #:endfor - end interface svdvals + + !> List of internal GESDD tasks: !> Return full matrices U, V^T to separate storage character, parameter :: GESDD_FULL_MATRICES = 'A' @@ -46,9 +24,6 @@ module stdlib_linalg_svd !> Do not return either U or VT (singular values array only) character, parameter :: GESDD_SINGVAL_ONLY = 'N' - - - contains !> Process GESDD output flag @@ -87,7 +62,7 @@ module stdlib_linalg_svd #:if rk!="xdp" !> Singular values of matrix A - function stdlib_linalg_svdvals_${ri}$(a,err) result(s) + module function stdlib_linalg_svdvals_${ri}$(a,err) result(s) !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] state return flag. On error if not requested, the code will stop @@ -115,7 +90,7 @@ module stdlib_linalg_svd end function stdlib_linalg_svdvals_${ri}$ !> SVD of matrix A = U S V^T, returning S and optionally U and V^T - subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err) + module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err) !> Input matrix A[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Array of singular values @@ -285,4 +260,4 @@ module stdlib_linalg_svd #:endif #:endfor -end module stdlib_linalg_svd +end submodule stdlib_linalg_svd diff --git a/test/linalg/test_linalg_svd.fypp b/test/linalg/test_linalg_svd.fypp index eb3ec7e29..088c6db14 100644 --- a/test/linalg/test_linalg_svd.fypp +++ b/test/linalg/test_linalg_svd.fypp @@ -4,10 +4,8 @@ module test_linalg_svd use testdrive, only: error_type, check, new_unittest, unittest_type use stdlib_linalg_constants - use stdlib_linalg, only: diag - use stdlib_linalg_svd, only: svd,svdvals + use stdlib_linalg, only: diag,svd,svdvals use stdlib_linalg_state, only: linalg_state_type - implicit none (type,external) public :: test_svd From 6abe7dff5a77c870fb7289b0a1f81c4304a4b9b3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 10:48:41 +0200 Subject: [PATCH 08/26] document `svd` and interface --- src/stdlib_linalg.fypp | 50 ++++++++++++++++++++++++++++++++++---- src/stdlib_linalg_svd.fypp | 22 +++++++++++++++-- 2 files changed, 65 insertions(+), 7 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index bc9ce5950..57f3f24ff 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -458,24 +458,64 @@ module stdlib_linalg #:endfor end interface - ! Singular value decomposition - interface svd + ! Singular value decomposition + interface svd !! version: experimental !! + !! Computes the singular value decomposition of a `real` or `complex` 2d matrix. + !! + !!### Summary + !! Interface for computing the singular value decomposition of a `real` or `complex` 2d matrix. + !! !!### Description !! + !! This interface provides methods for computing the singular value decomposition of a matrix. + !! Supported data types include `real` and `complex`. The subroutine returns a `real` array of + !! singular values, and optionally, left- and right- singular vector matrices, `U` and `V`. + !! For a matrix `A` with size [m,n], full matrix storage for `U` and `V` should be [m,m] and [n,n]. + !! It is possible to use partial storage [m,k] and [k,n], `k=min(m,n)`, choosing `full_matrices=.false.`. !! - + !!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + !!### Example + !! + !!```fortran + !! real(sp) :: a(2,3), s(2), u(2,2), vt(3,3) + !! a = reshape([3,2, 2,3, 2,-2],[2,3]) + !! + !! call svd(A,s,u,v) + !! print *, 'singular values = ',s + !!``` + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err) + !!### Summary + !! Compute singular value decomposition of a matrix \( A = U \cdot S \cdot \V^T \) + !! + !!### Description + !! + !! This function computes the singular value decomposition of a `real` or `complex` matrix \( A \), + !! and returns the array of singular values, and optionally the left matrix \( U \) containing the + !! left unitary singular vectors, and the right matrix \( V^T \), containing the right unitary + !! singular vectors. + !! + !! param: a Input matrix of size [m,n]. + !! param: s Output `real` array of size [min(m,n)] returning a list of singular values. + !! param: u [optional] Output left singular matrix of size [m,m] or [m,min(m,n)] (.not.full_matrices). Contains singular vectors as columns. + !! param: vt [optional] Output right singular matrix of size [n,n] or [min(m,n),n] (.not.full_matrices). Contains singular vectors as rows. + !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten. + !! param: full_matrices [optional] If `.true.` (default), matrices \( U \) and \( V^T \) have size [m,m], [n,n]. Otherwise, they are [m,k], [k,n] with `k=min(m,n)`. + !! param: err [optional] State return flag. + !! !> Input matrix A[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Array of singular values real(${rk}$), intent(out) :: s(:) - !> The columns of U contain the eigenvectors of A A^T + !> The columns of U contain the left singular vectors ${rt}$, optional, intent(out), target :: u(:,:) - !> The rows of V^T contain the eigenvectors of A^T A + !> The rows of V^T contain the right singular vectors ${rt}$, optional, intent(out), target :: vt(:,:) !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index 8bc72ffdc..f8538e134 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -91,13 +91,31 @@ submodule(stdlib_linalg) stdlib_linalg_svd !> SVD of matrix A = U S V^T, returning S and optionally U and V^T module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err) + !!### Summary + !! Compute singular value decomposition of a matrix \( A = U \cdot S \cdot \V^T \) + !! + !!### Description + !! + !! This function computes the singular value decomposition of a `real` or `complex` matrix \( A \), + !! and returns the array of singular values, and optionally the left matrix \( U \) containing the + !! left unitary singular vectors, and the right matrix \( V^T \), containing the right unitary + !! singular vectors. + !! + !! param: a Input matrix of size [m,n]. + !! param: s Output `real` array of size [min(m,n)] returning a list of singular values. + !! param: u [optional] Output left singular matrix of size [m,m] or [m,min(m,n)] (.not.full_matrices). Contains singular vectors as columns. + !! param: vt [optional] Output right singular matrix of size [n,n] or [min(m,n),n] (.not.full_matrices). Contains singular vectors as rows. + !! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten. + !! param: full_matrices [optional] If `.true.` (default), matrices \( U \) and \( V^T \) have size [m,m], [n,n]. Otherwise, they are [m,k], [k,n] with `k=min(m,n)`. + !! param: err [optional] State return flag. + !! !> Input matrix A[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Array of singular values real(${rk}$), intent(out) :: s(:) - !> The columns of U contain the eigenvectors of A A^T + !> The columns of U contain the left singular vectors ${rt}$, optional, intent(out), target :: u(:,:) - !> The rows of V^T contain the eigenvectors of A^T A + !> The rows of V^T contain the right singular vectors ${rt}$, optional, intent(out), target :: vt(:,:) !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a From 194f4a1c86052154785997b58eae90f5d4024089 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 10:56:02 +0200 Subject: [PATCH 09/26] document `svdvals` and interface --- src/stdlib_linalg.fypp | 38 ++++++++++++++++++++++++++++++++++++++ src/stdlib_linalg_svd.fypp | 15 +++++++++++++++ 2 files changed, 53 insertions(+) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 57f3f24ff..9f261f4de 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -533,12 +533,50 @@ module stdlib_linalg interface svdvals !! version: experimental !! + !! Computes the singular values of a `real` or `complex` 2d matrix. + !! + !!### Summary + !! + !! Function interface for computing the array of singular values from the singular value decomposition + !! of a `real` or `complex` 2d matrix. + !! !!### Description !! + !! This interface provides methods for computing the singular values a 2d matrix. + !! Supported data types include `real` and `complex`. The function returns a `real` array of + !! singular values, with size [min(m,n)]. !! + !!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! + !!### Example + !! + !!```fortran + !! real(sp) :: a(2,3), s(2) + !! a = reshape([3,2, 2,3, 2,-2],[2,3]) + !! + !! s = svdvals(A) + !! print *, 'singular values = ',s + !!``` + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module function stdlib_linalg_svdvals_${ri}$(a,err) result(s) + !!### Summary + !! Compute singular values \(S \) from the singular-value decomposition of a matrix \( A = U \cdot S \cdot \V^T \). + !! + !!### Description + !! + !! This function returns the array of singular values from the singular value decomposition of a `real` + !! or `complex` matrix \( A = U \cdot S \cdot V^T \). + !! + !! param: a Input matrix of size [m,n]. + !! param: err [optional] State return flag. + !! + !!### Return value + !! + !! param: s `real` array of size [min(m,n)] returning a list of singular values. + !! !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] state return flag. On error if not requested, the code will stop diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index f8538e134..74135a3e2 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -63,6 +63,21 @@ submodule(stdlib_linalg) stdlib_linalg_svd !> Singular values of matrix A module function stdlib_linalg_svdvals_${ri}$(a,err) result(s) + !!### Summary + !! Compute singular values \(S \) from the singular-value decomposition of a matrix \( A = U \cdot S \cdot \V^T \). + !! + !!### Description + !! + !! This function returns the array of singular values from the singular value decomposition of a `real` + !! or `complex` matrix \( A = U \cdot S \cdot V^T \). + !! + !! param: a Input matrix of size [m,n]. + !! param: err [optional] State return flag. + !! + !!### Return value + !! + !! param: s `real` array of size [min(m,n)] returning a list of singular values. + !! !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] state return flag. On error if not requested, the code will stop From 9de160b85815205fec0e008ca19d49a3f6a004ee Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 10:58:16 +0200 Subject: [PATCH 10/26] remove `goto` --- src/stdlib_linalg_svd.fypp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index 74135a3e2..1f6bb38e1 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -159,13 +159,13 @@ submodule(stdlib_linalg) stdlib_linalg_svd if (.not.k>0) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size: a=',[m,n]) - goto 1 - end if - - if (.not.size(s,kind=ilp)>=k) then + call linalg_error_handling(err0,err) + return + elseif (.not.size(s,kind=ilp)>=k) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'singular value array has insufficient size:',& ' s=',shape(s,kind=ilp),', k=',k) - goto 1 + call linalg_error_handling(err0,err) + return endif ! Integer storage @@ -287,7 +287,7 @@ submodule(stdlib_linalg) stdlib_linalg_svd if (copy_a) deallocate(amat) if (alloc_u) deallocate(umat) if (alloc_vt) deallocate(vtmat) - 1 call linalg_error_handling(err0,err) + call linalg_error_handling(err0,err) end subroutine stdlib_linalg_svd_${ri}$ #:endif From cf40f75392a41a96c311d1846dc2b9e9f99d6cca Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 11:16:08 +0200 Subject: [PATCH 11/26] add example --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_svd1.f90 | 50 +++++++++++++++++++++++++++++++++ 2 files changed, 51 insertions(+) create mode 100644 example/linalg/example_svd1.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 729f117d3..3261ff6d9 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -20,5 +20,6 @@ ADD_EXAMPLE(blas_gemv) ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) +ADD_EXAMPLE(svd1) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) diff --git a/example/linalg/example_svd1.f90 b/example/linalg/example_svd1.f90 new file mode 100644 index 000000000..261d5ee40 --- /dev/null +++ b/example/linalg/example_svd1.f90 @@ -0,0 +1,50 @@ +! Singular Value Decomposition +program example_svd1 + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: svd + implicit none + + real(dp), allocatable :: A(:,:),s(:),u(:,:),vt(:,:) + + ! We want to find the singular value decomposition of matrix: + ! + ! A = [ 3 2 2] + ! [ 2 3 -2] + ! + A = transpose(reshape([ 3, 2, 2, & + 2, 3,-2], [3,2])) + + ! Prepare arrays + allocate(s(2),u(2,2),vt(3,3)) + + ! Get singular value decomposition + call svd(A,s,u,vt) + + ! Singular values: [5, 3] + print 1, ' ' + print 1, 'S = ',s + print 1, ' ' + + ! Left vectors (may be flipped): + ! [Ã2/2 Ã2/2] + ! U = [Ã2/2 -Ã2/2] + ! + print 1, ' ' + print 1, 'U = ',u(1,:) + print 1, ' ',u(2,:) + + + ! Right vectors (may be flipped): + ! [Ã2/2 Ã2/2 0] + ! V = [1/Ã18 -1/Ã18 4/Ã18] + ! [ 2/3 -2/3 -1/3] + ! + print 1, ' ' + print 1, ' ',vt(1,:) + print 1, 'VT= ',vt(2,:) + print 1, ' ',vt(3,:) + print 1, ' ' + + 1 format(a,*(1x,f12.8)) + +end program example_svd1 From 32784cd93a25414e0417f354375c1a6da6b35d35 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 11:45:17 +0200 Subject: [PATCH 12/26] add `svdvals` example --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_svd2.f90 | 27 +++++++++++++++++++++++++++ 2 files changed, 28 insertions(+) create mode 100644 example/linalg/example_svd2.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 3261ff6d9..c8d204faa 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -21,5 +21,6 @@ ADD_EXAMPLE(lapack_getrf) ADD_EXAMPLE(lstsq1) ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(svd1) +ADD_EXAMPLE(svd2) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) diff --git a/example/linalg/example_svd2.f90 b/example/linalg/example_svd2.f90 new file mode 100644 index 000000000..1b22d5f9e --- /dev/null +++ b/example/linalg/example_svd2.f90 @@ -0,0 +1,27 @@ +! Singular Values +program example_svd2 + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: svd + implicit none + + real(dp), allocatable :: A(:,:),s(:) + + ! We want to find the singular values of matrix: + ! + ! A = [ 3 2 2] + ! [ 2 3 -2] + ! + A = transpose(reshape([ 3, 2, 2, & + 2, 3,-2], [3,2])) + + ! Get singular values + s = svdvals(A) + + ! Singular values: [5, 3] + print 1, ' ' + print 1, 'S = ',s + print 1, ' ' + + 1 format(a,*(1x,f12.8)) + +end program example_svd2 From bb47ba53e49da7a6026cebe935344414f1d05d8c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 11:46:43 +0200 Subject: [PATCH 13/26] specs: `svd`, `svdvals` --- doc/specs/stdlib_linalg.md | 89 ++++++++++++++++++++++++++++++++++++++ src/stdlib_linalg.fypp | 2 + 2 files changed, 91 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index a081291d0..8c4eadfdf 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -797,3 +797,92 @@ Exceptions trigger an `error stop`. ```fortran {!example/linalg/example_determinant2.f90!} ``` + +## `svd` - Compute the singular value decomposition of a 2d matrix. + +### Status + +Experimental + +### Description + +This subroutine computes the singular value decomposition of a `real` or `complex` 2d matrix \( A = U \cdot S \cdot \V^T \). +The solver is based on LAPACK's `*GESDD` backends. + +Result vector `s` returns the array of singular values on the diagonal of \( S \). +If requested, `u` contains the left singular vectors, as columns of \( U \). +If requested, `vt` contains the right singular vectors, as rows of \( V^T \). + +### Syntax + +`call ` [[stdlib_linalg(module):svd(interface)]] `(a, s, [, u, vt, overwrite_a, full_matrices, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(inout)` argument, but returns unchanged unless `overwrite_a=.true.`. + +`s`: Shall be a rank-1 `real` array, returning the list of `k = min(m,n)` singular values. It is an `intent(inout)` argument. + +`u` (optional): Shall be a rank-2 array of same kind as `a`, returning the left singular vectors of `a` as columns. Its size should be `[m,m]` unless `full_matrices=.false.`, in which case, it can be `[m,min(m,n)]`. It is an `intent(inout)` argument. + +`vt` (optional): Shall be a rank-2 array of same kind as `a`, returning the right singular vectors of `a` as rows. Its size should be `[n,n]` unless `full_matrices=.false.`, in which case, it can be `[min(m,n),n]`. It is an `intent(inout)` 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. By default, `overwrite_a=.false.`. This is an `intent(in)` argument. + +`full_matrices` (optional): Shall be an input `logical` flag. If `.true.` (default), matrices `u` and `vt` shall be full-sized. Otherwise, their secondary dimension can be resized to `min(m,n)`. See `u`, `v` for details. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return values + +Returns an array `s` that contains the list of singular values of matrix `a`. +If requested, returns a rank-2 array `u` that contains the left singular values of `a` along its columns. +If requested, returns a rank-2 array `vt` that contains the right singular values of `a` along its rows. + +Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. +Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes. +Exceptions trigger an `error stop`, unless argument `err` is present. + +### Example + +```fortran +{!example/linalg/example_svd1.f90!} +``` + +## `svdvals` - Compute the singular values of a 2d matrix. + +### Status + +Experimental + +### Description + +This subroutine computes the singular values of a `real` or `complex` 2d matrix from its singular +value decomposition \( A = U \cdot S \cdot \V^T \). The solver is based on LAPACK's `*GESDD` backends. + +Result vector `s` returns the array of singular values on the diagonal of \( S \). + +### Syntax + +`s = ` [[stdlib_linalg(module):svdvals(interface)]] `(a [, err])` + +### Arguments + +`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return values + +Returns an array `s` that contains the list of singular values of matrix `a`. + +Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. +Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes. +Exceptions trigger an `error stop`, unless argument `err` is present. + +### Example + +```fortran +{!example/linalg/example_svd2.f90!} +``` + diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 9f261f4de..593f1f57b 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -463,6 +463,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the singular value decomposition of a `real` or `complex` 2d matrix. + !! ([Specification](../page/specs/stdlib_linalg.html#svd-compute-the-singular-value-decomposition-of-a-2d-matrix)) !! !!### Summary !! Interface for computing the singular value decomposition of a `real` or `complex` 2d matrix. @@ -534,6 +535,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the singular values of a `real` or `complex` 2d matrix. + !! ([Specification](../page/specs/stdlib_linalg.html#svdvals-compute-the-singular-values-of-a-2d-matrix)) !! !!### Summary !! From 3b365aac1d00d2070de66980552b572621db5421 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 12:50:39 +0200 Subject: [PATCH 14/26] fix --- example/linalg/example_svd2.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/linalg/example_svd2.f90 b/example/linalg/example_svd2.f90 index 1b22d5f9e..f483c1118 100644 --- a/example/linalg/example_svd2.f90 +++ b/example/linalg/example_svd2.f90 @@ -1,7 +1,7 @@ ! Singular Values program example_svd2 use stdlib_linalg_constants, only: dp - use stdlib_linalg, only: svd + use stdlib_linalg, only: svdvals implicit none real(dp), allocatable :: A(:,:),s(:) From ba5e5e7078a4637f17ba70eabd3abc14ac2864ff Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 20:04:09 +0200 Subject: [PATCH 15/26] remove complex qp test --- test/linalg/test_linalg_svd.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/linalg/test_linalg_svd.fypp b/test/linalg/test_linalg_svd.fypp index 088c6db14..273c884f4 100644 --- a/test/linalg/test_linalg_svd.fypp +++ b/test/linalg/test_linalg_svd.fypp @@ -26,7 +26,7 @@ module test_linalg_svd #:endfor #:for ck,ct,ci in CMPLX_KINDS_TYPES - #:if ck!="xdp" + #:if not ck in ["xdp","qp"] tests = [tests,new_unittest("test_complex_svd_${ci}$",test_complex_svd_${ci}$)] #:endif #:endfor @@ -175,7 +175,7 @@ module test_linalg_svd !> Test complex svd #:for ck,ct,ci in CMPLX_KINDS_TYPES - #:if ck!="xdp" + #:if not ck in ["xdp","qp"] subroutine test_complex_svd_${ci}$(error) type(error_type), allocatable, intent(out) :: error From 6556e8b52d6157f917159284cb296e14c34ffed0 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 20:13:27 +0200 Subject: [PATCH 16/26] restore `qp` subroutine (test still unactive) --- test/linalg/test_linalg_svd.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/linalg/test_linalg_svd.fypp b/test/linalg/test_linalg_svd.fypp index 273c884f4..db88df9c4 100644 --- a/test/linalg/test_linalg_svd.fypp +++ b/test/linalg/test_linalg_svd.fypp @@ -175,7 +175,7 @@ module test_linalg_svd !> Test complex svd #:for ck,ct,ci in CMPLX_KINDS_TYPES - #:if not ck in ["xdp","qp"] + #:if ck!="xdp" subroutine test_complex_svd_${ci}$(error) type(error_type), allocatable, intent(out) :: error From 9c2f1332b942349d54119705776bf020a012d052 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 9 May 2024 21:46:22 +0200 Subject: [PATCH 17/26] restore all, avoid `real128` complex parameter math --- test/linalg/test_linalg_svd.fypp | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/test/linalg/test_linalg_svd.fypp b/test/linalg/test_linalg_svd.fypp index db88df9c4..9fe8a889a 100644 --- a/test/linalg/test_linalg_svd.fypp +++ b/test/linalg/test_linalg_svd.fypp @@ -26,7 +26,7 @@ module test_linalg_svd #:endfor #:for ck,ct,ci in CMPLX_KINDS_TYPES - #:if not ck in ["xdp","qp"] + #:if ck!="xdp" tests = [tests,new_unittest("test_complex_svd_${ci}$",test_complex_svd_${ci}$)] #:endif #:endfor @@ -180,18 +180,20 @@ module test_linalg_svd type(error_type), allocatable, intent(out) :: error !> Reference solution - real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) - real(${ck}$), parameter :: one = 1.0_${ck}$ - real(${ck}$), parameter :: zero = 0.0_${ck}$ - real(${ck}$), parameter :: sqrt2 = sqrt(2.0_${ck}$) - real(${ck}$), parameter :: rsqrt2 = one/sqrt2 - ${ct}$, parameter :: cone = (1.0_${ck}$,0.0_${ck}$) - ${ct}$, parameter :: cimg = (0.0_${ck}$,1.0_${ck}$) - ${ct}$, parameter :: czero = (0.0_${ck}$,0.0_${ck}$) + real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$)) + real(${ck}$), parameter :: one = 1.0_${ck}$ + real(${ck}$), parameter :: zero = 0.0_${ck}$ + real(${ck}$), parameter :: sqrt2 = sqrt(2.0_${ck}$) + real(${ck}$), parameter :: rsqrt2 = one/sqrt2 + ${ct}$, parameter :: csqrt2 = (rsqrt2,zero) + ${ct}$, parameter :: isqrt2 = (zero,rsqrt2) + ${ct}$, parameter :: cone = (1.0_${ck}$,0.0_${ck}$) + ${ct}$, parameter :: cimg = (0.0_${ck}$,1.0_${ck}$) + ${ct}$, parameter :: czero = (0.0_${ck}$,0.0_${ck}$) real(${ck}$), parameter :: s_sol(2) = [sqrt2,sqrt2] ${ct}$, parameter :: A_mat(2,2) = reshape([cone,cimg,cimg,cone],[2,2]) - ${ct}$, parameter :: u_sol(2,2) = reshape(rsqrt2*[cone,cimg,cimg,cone],[2,2]) + ${ct}$, parameter :: u_sol(2,2) = reshape([csqrt2,isqrt2,isqrt2,csqrt2],[2,2]) ${ct}$, parameter :: vt_sol(2,2) = reshape([cone,czero,czero,cone],[2,2]) !> Local variables From cf15eb1636f53a2342301d028ed119ef567f6d76 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Sun, 26 May 2024 12:33:15 +0200 Subject: [PATCH 18/26] clearer `logical` flags --- src/stdlib_linalg_svd.fypp | 57 ++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 30 deletions(-) diff --git a/src/stdlib_linalg_svd.fypp b/src/stdlib_linalg_svd.fypp index 1f6bb38e1..b93343c15 100644 --- a/src/stdlib_linalg_svd.fypp +++ b/src/stdlib_linalg_svd.fypp @@ -144,7 +144,7 @@ submodule(stdlib_linalg) stdlib_linalg_svd type(linalg_state_type) :: err0 integer(ilp) :: m,n,lda,ldu,ldvt,info,k,lwork,liwork,lrwork integer(ilp), allocatable :: iwork(:) - logical(lk) :: copy_a,full_storage,compute_uv,alloc_u,alloc_vt,can_overwrite_a + logical(lk) :: overwrite_a_,full_storage,compute_uv,temp_u,temp_vt,can_overwrite_amat character :: task ${rt}$, target :: work_dummy(1),u_dummy(1,1),vt_dummy(1,1) ${rt}$, allocatable :: work(:) @@ -173,21 +173,18 @@ submodule(stdlib_linalg) stdlib_linalg_svd allocate(iwork(liwork)) ! Can A be overwritten? By default, do not overwrite - if (present(overwrite_a)) then - copy_a = .not.overwrite_a - else - copy_a = .true._lk - endif + overwrite_a_ = .false. + if (present(overwrite_a)) overwrite_a_ = overwrite_a - ! Initialize a matrix temporary - if (copy_a) then - allocate(amat(m,n),source=a) - else + ! Initialize a matrix temporary? + if (overwrite_a_) then amat => a + else + allocate(amat(m,n),source=a) endif - ! Check if we can overwrite A with data that will be lost - can_overwrite_a = copy_a .and. merge(.not.present(u),.not.present(vt),m>=n) + ! Check if we can overwrite amat with data that will be lost + can_overwrite_amat = (.not.overwrite_a_) .and. merge(.not.present(u),.not.present(vt),m>=n) ! Full-size matrices if (present(full_matrices)) then @@ -202,38 +199,38 @@ submodule(stdlib_linalg) stdlib_linalg_svd ! U, VT storage if (present(u)) then ! User input - umat => u - alloc_u = .false. - elseif ((copy_a .and. m>=n) .or. .not.compute_uv) then + umat => u + temp_u = .false. + elseif ((m>=n .and. .not.overwrite_a_) .or. .not.compute_uv) then ! U not wanted, and A can be overwritten: do not allocate - umat => u_dummy - alloc_u = .false. + umat => u_dummy + temp_u = .false. elseif (.not.full_storage) then ! Allocate with minimum size allocate(umat(m,k)) - alloc_u = .true. + temp_u = .true. else ! Allocate with regular size allocate(umat(m,m)) - alloc_u = .true. + temp_u = .true. end if if (present(vt)) then ! User input - vtmat => vt - alloc_vt = .false. - elseif ((copy_a .and. m vt + temp_vt = .false. + elseif ((m vt_dummy - alloc_vt = .false. + vtmat => vt_dummy + temp_vt = .false. elseif (.not.full_storage) then ! Allocate with minimum size allocate(vtmat(k,n)) - alloc_vt = .true. + temp_vt = .true. else ! Allocate with regular size allocate(vtmat(n,n)) - alloc_vt = .true. + temp_vt = .true. end if ldu = size(umat ,1,kind=ilp) @@ -242,7 +239,7 @@ submodule(stdlib_linalg) stdlib_linalg_svd ! Decide SVD task if (.not.compute_uv) then task = GESDD_SINGVAL_ONLY - elseif (can_overwrite_a) then + elseif (can_overwrite_amat) then ! A is a copy: we can overwrite its storage task = GESDD_OVERWRITE_A elseif (.not.full_storage) then @@ -284,9 +281,9 @@ submodule(stdlib_linalg) stdlib_linalg_svd endif ! Finalize storage and process output flag - if (copy_a) deallocate(amat) - if (alloc_u) deallocate(umat) - if (alloc_vt) deallocate(vtmat) + if (.not.overwrite_a_) deallocate(amat) + if (temp_u) deallocate(umat) + if (temp_vt) deallocate(vtmat) call linalg_error_handling(err0,err) end subroutine stdlib_linalg_svd_${ri}$ From 980a555f037820e54c93a5d2eaf3943dad39ae7b Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 21:55:34 +0200 Subject: [PATCH 19/26] 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 866859352..5eb0ab751 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -928,7 +928,7 @@ If requested, `vt` contains the right singular vectors, as rows of \( V^T \). `vt` (optional): Shall be a rank-2 array of same kind as `a`, returning the right singular vectors of `a` as rows. Its size should be `[n,n]` unless `full_matrices=.false.`, in which case, it can be `[min(m,n),n]`. It is an `intent(inout)` 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. By default, `overwrite_a=.false.`. This is an `intent(in)` 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. By default, `overwrite_a=.false.`. It is an `intent(in)` argument. `full_matrices` (optional): Shall be an input `logical` flag. If `.true.` (default), matrices `u` and `vt` shall be full-sized. Otherwise, their secondary dimension can be resized to `min(m,n)`. See `u`, `v` for details. From 28ae2cb661a9147e9c257e30f97830a21bc1292f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 21:56:14 +0200 Subject: [PATCH 20/26] 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 5eb0ab751..eab8495a8 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -937,7 +937,7 @@ If requested, `vt` contains the right singular vectors, as rows of \( V^T \). ### Return values Returns an array `s` that contains the list of singular values of matrix `a`. -If requested, returns a rank-2 array `u` that contains the left singular values of `a` along its columns. +If requested, returns a rank-2 array `u` that contains the left singular vectors of `a` along its columns. If requested, returns a rank-2 array `vt` that contains the right singular values of `a` along its rows. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. From 69458ad88c53bbad4824e29ec870a4f008257339 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 21:56:29 +0200 Subject: [PATCH 21/26] 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 eab8495a8..fa0997c30 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -938,7 +938,7 @@ If requested, `vt` contains the right singular vectors, as rows of \( V^T \). Returns an array `s` that contains the list of singular values of matrix `a`. If requested, returns a rank-2 array `u` that contains the left singular vectors of `a` along its columns. -If requested, returns a rank-2 array `vt` that contains the right singular values of `a` along its rows. +If requested, returns a rank-2 array `vt` that contains the right singular vectors of `a` along its rows. Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge. Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes. From b579d9866cb51fb6b66a968c78b20a26dd88f083 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 21:59:18 +0200 Subject: [PATCH 22/26] rename `svd` examples --- doc/specs/stdlib_linalg.md | 4 ++-- example/linalg/CMakeLists.txt | 4 ++-- example/linalg/{example_svd1.f90 => example_svd.f90} | 4 ++-- example/linalg/{example_svd2.f90 => example_svdvals.f90} | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) rename example/linalg/{example_svd1.f90 => example_svd.f90} (95%) rename example/linalg/{example_svd2.f90 => example_svdvals.f90} (90%) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index fa0997c30..0b0ecaacf 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -947,7 +947,7 @@ Exceptions trigger an `error stop`, unless argument `err` is present. ### Example ```fortran -{!example/linalg/example_svd1.f90!} +{!example/linalg/example_svd.f90!} ``` ## `svdvals` - Compute the singular values of a 2d matrix. @@ -984,6 +984,6 @@ Exceptions trigger an `error stop`, unless argument `err` is present. ### Example ```fortran -{!example/linalg/example_svd2.f90!} +{!example/linalg/example_svdvals.f90!} ``` diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index ce9cf1315..e9ad22a03 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -23,7 +23,7 @@ ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) -ADD_EXAMPLE(svd1) -ADD_EXAMPLE(svd2) +ADD_EXAMPLE(svd) +ADD_EXAMPLE(svdvals) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) diff --git a/example/linalg/example_svd1.f90 b/example/linalg/example_svd.f90 similarity index 95% rename from example/linalg/example_svd1.f90 rename to example/linalg/example_svd.f90 index 261d5ee40..1f03f56e2 100644 --- a/example/linalg/example_svd1.f90 +++ b/example/linalg/example_svd.f90 @@ -1,5 +1,5 @@ ! Singular Value Decomposition -program example_svd1 +program example_svd use stdlib_linalg_constants, only: dp use stdlib_linalg, only: svd implicit none @@ -47,4 +47,4 @@ program example_svd1 1 format(a,*(1x,f12.8)) -end program example_svd1 +end program example_svd diff --git a/example/linalg/example_svd2.f90 b/example/linalg/example_svdvals.f90 similarity index 90% rename from example/linalg/example_svd2.f90 rename to example/linalg/example_svdvals.f90 index f483c1118..ef0d755f5 100644 --- a/example/linalg/example_svd2.f90 +++ b/example/linalg/example_svdvals.f90 @@ -1,5 +1,5 @@ ! Singular Values -program example_svd2 +program example_svdvals use stdlib_linalg_constants, only: dp use stdlib_linalg, only: svdvals implicit none @@ -24,4 +24,4 @@ program example_svd2 1 format(a,*(1x,f12.8)) -end program example_svd2 +end program example_svdvals From 8387797c86192bed56701fdf01572194a86c0d79 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 22:02:04 +0200 Subject: [PATCH 23/26] fix intent in docs --- doc/specs/stdlib_linalg.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 0b0ecaacf..13f443887 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -922,11 +922,11 @@ If requested, `vt` contains the right singular vectors, as rows of \( V^T \). `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(inout)` argument, but returns unchanged unless `overwrite_a=.true.`. -`s`: Shall be a rank-1 `real` array, returning the list of `k = min(m,n)` singular values. It is an `intent(inout)` argument. +`s`: Shall be a rank-1 `real` array, returning the list of `k = min(m,n)` singular values. It is an `intent(out)` argument. -`u` (optional): Shall be a rank-2 array of same kind as `a`, returning the left singular vectors of `a` as columns. Its size should be `[m,m]` unless `full_matrices=.false.`, in which case, it can be `[m,min(m,n)]`. It is an `intent(inout)` argument. +`u` (optional): Shall be a rank-2 array of same kind as `a`, returning the left singular vectors of `a` as columns. Its size should be `[m,m]` unless `full_matrices=.false.`, in which case, it can be `[m,min(m,n)]`. It is an `intent(out)` argument. -`vt` (optional): Shall be a rank-2 array of same kind as `a`, returning the right singular vectors of `a` as rows. Its size should be `[n,n]` unless `full_matrices=.false.`, in which case, it can be `[min(m,n),n]`. It is an `intent(inout)` argument. +`vt` (optional): Shall be a rank-2 array of same kind as `a`, returning the right singular vectors of `a` as rows. Its size should be `[n,n]` unless `full_matrices=.false.`, in which case, it can be `[min(m,n),n]`. 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. By default, `overwrite_a=.false.`. It is an `intent(in)` argument. From e0d800c949cae73cdd30b3a62c61592264354f43 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 22:03:07 +0200 Subject: [PATCH 24/26] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index fa0997c30..7289d4f41 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -918,6 +918,9 @@ If requested, `vt` contains the right singular vectors, as rows of \( V^T \). `call ` [[stdlib_linalg(module):svd(interface)]] `(a, s, [, u, vt, overwrite_a, full_matrices, err])` +### Class +Subroutine + ### Arguments `a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(inout)` argument, but returns unchanged unless `overwrite_a=.true.`. From 8f35e9856aaba14f1e715d647896d89037bddaed Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 22:08:27 +0200 Subject: [PATCH 25/26] format -> character parameter --- example/linalg/example_svd.f90 | 24 ++++++++++++------------ example/linalg/example_svdvals.f90 | 9 ++++----- 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/example/linalg/example_svd.f90 b/example/linalg/example_svd.f90 index 1f03f56e2..ba8a14dbb 100644 --- a/example/linalg/example_svd.f90 +++ b/example/linalg/example_svd.f90 @@ -5,6 +5,7 @@ program example_svd implicit none real(dp), allocatable :: A(:,:),s(:),u(:,:),vt(:,:) + character(*), parameter :: fmt = "(a,*(1x,f12.8))" ! We want to find the singular value decomposition of matrix: ! @@ -21,17 +22,17 @@ program example_svd call svd(A,s,u,vt) ! Singular values: [5, 3] - print 1, ' ' - print 1, 'S = ',s - print 1, ' ' + print fmt, ' ' + print fmt, 'S = ',s + print fmt, ' ' ! Left vectors (may be flipped): ! [Ã2/2 Ã2/2] ! U = [Ã2/2 -Ã2/2] ! - print 1, ' ' - print 1, 'U = ',u(1,:) - print 1, ' ',u(2,:) + print fmt, ' ' + print fmt, 'U = ',u(1,:) + print fmt, ' ',u(2,:) ! Right vectors (may be flipped): @@ -39,12 +40,11 @@ program example_svd ! V = [1/Ã18 -1/Ã18 4/Ã18] ! [ 2/3 -2/3 -1/3] ! - print 1, ' ' - print 1, ' ',vt(1,:) - print 1, 'VT= ',vt(2,:) - print 1, ' ',vt(3,:) - print 1, ' ' + print fmt, ' ' + print fmt, ' ',vt(1,:) + print fmt, 'VT= ',vt(2,:) + print fmt, ' ',vt(3,:) + print fmt, ' ' - 1 format(a,*(1x,f12.8)) end program example_svd diff --git a/example/linalg/example_svdvals.f90 b/example/linalg/example_svdvals.f90 index ef0d755f5..fe228c16e 100644 --- a/example/linalg/example_svdvals.f90 +++ b/example/linalg/example_svdvals.f90 @@ -5,6 +5,7 @@ program example_svdvals implicit none real(dp), allocatable :: A(:,:),s(:) + character(*), parameter :: fmt="(a,*(1x,f12.8))" ! We want to find the singular values of matrix: ! @@ -18,10 +19,8 @@ program example_svdvals s = svdvals(A) ! Singular values: [5, 3] - print 1, ' ' - print 1, 'S = ',s - print 1, ' ' + print fmt, ' ' + print fmt, 'S = ',s + print fmt, ' ' - 1 format(a,*(1x,f12.8)) - end program example_svdvals From bdffb4f52c9dcf07c30e15056e4e1f99ca7912ee Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Tue, 28 May 2024 22:12:04 +0200 Subject: [PATCH 26/26] matrix -> rank-2 array --- doc/specs/stdlib_linalg.md | 8 ++++---- src/stdlib_linalg.fypp | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 5a4da2405..395d4d6b3 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -899,7 +899,7 @@ Exceptions trigger an `error stop`. {!example/linalg/example_determinant2.f90!} ``` -## `svd` - Compute the singular value decomposition of a 2d matrix. +## `svd` - Compute the singular value decomposition of a rank-2 array (matrix). ### Status @@ -907,7 +907,7 @@ Experimental ### Description -This subroutine computes the singular value decomposition of a `real` or `complex` 2d matrix \( A = U \cdot S \cdot \V^T \). +This subroutine computes the singular value decomposition of a `real` or `complex` rank-2 array (matrix) \( A = U \cdot S \cdot \V^T \). The solver is based on LAPACK's `*GESDD` backends. Result vector `s` returns the array of singular values on the diagonal of \( S \). @@ -953,7 +953,7 @@ Exceptions trigger an `error stop`, unless argument `err` is present. {!example/linalg/example_svd.f90!} ``` -## `svdvals` - Compute the singular values of a 2d matrix. +## `svdvals` - Compute the singular values of a rank-2 array (matrix). ### Status @@ -961,7 +961,7 @@ Experimental ### Description -This subroutine computes the singular values of a `real` or `complex` 2d matrix from its singular +This subroutine computes the singular values of a `real` or `complex` rank-2 array (matrix) from its singular value decomposition \( A = U \cdot S \cdot \V^T \). The solver is based on LAPACK's `*GESDD` backends. Result vector `s` returns the array of singular values on the diagonal of \( S \). diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index f27b8fa15..10cd16650 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -559,7 +559,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the singular value decomposition of a `real` or `complex` 2d matrix. - !! ([Specification](../page/specs/stdlib_linalg.html#svd-compute-the-singular-value-decomposition-of-a-2d-matrix)) + !! ([Specification](../page/specs/stdlib_linalg.html#svd-compute-the-singular-value-decomposition-of-a-rank-2-array-matrix)) !! !!### Summary !! Interface for computing the singular value decomposition of a `real` or `complex` 2d matrix. @@ -631,7 +631,7 @@ module stdlib_linalg !! version: experimental !! !! Computes the singular values of a `real` or `complex` 2d matrix. - !! ([Specification](../page/specs/stdlib_linalg.html#svdvals-compute-the-singular-values-of-a-2d-matrix)) + !! ([Specification](../page/specs/stdlib_linalg.html#svdvals-compute-the-singular-values-of-a-rank-2-array-matrix)) !! !!### Summary !!