diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ffc4f7c42..a8fd633c7 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -898,3 +898,95 @@ Exceptions trigger an `error stop`. ```fortran {!example/linalg/example_determinant2.f90!} ``` + +## `svd` - Compute the singular value decomposition of a rank-2 array (matrix). + +### Status + +Experimental + +### Description + +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 \). +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])` + +### 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.`. + +`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(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(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. + +`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 vectors of `a` along its columns. +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. +Exceptions trigger an `error stop`, unless argument `err` is present. + +### Example + +```fortran +{!example/linalg/example_svd.f90!} +``` + +## `svdvals` - Compute the singular values of a rank-2 array (matrix). + +### Status + +Experimental + +### Description + +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 \). + +### 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_svdvals.f90!} +``` + diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 2d80f067c..e9ad22a03 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -23,5 +23,7 @@ ADD_EXAMPLE(lstsq2) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) +ADD_EXAMPLE(svd) +ADD_EXAMPLE(svdvals) ADD_EXAMPLE(determinant) ADD_EXAMPLE(determinant2) diff --git a/example/linalg/example_svd.f90 b/example/linalg/example_svd.f90 new file mode 100644 index 000000000..ba8a14dbb --- /dev/null +++ b/example/linalg/example_svd.f90 @@ -0,0 +1,50 @@ +! Singular Value Decomposition +program example_svd + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: 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: + ! + ! 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 fmt, ' ' + print fmt, 'S = ',s + print fmt, ' ' + + ! Left vectors (may be flipped): + ! [Ã2/2 Ã2/2] + ! U = [Ã2/2 -Ã2/2] + ! + print fmt, ' ' + print fmt, 'U = ',u(1,:) + print fmt, ' ',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 fmt, ' ' + print fmt, ' ',vt(1,:) + print fmt, 'VT= ',vt(2,:) + print fmt, ' ',vt(3,:) + print fmt, ' ' + + +end program example_svd diff --git a/example/linalg/example_svdvals.f90 b/example/linalg/example_svdvals.f90 new file mode 100644 index 000000000..fe228c16e --- /dev/null +++ b/example/linalg/example_svdvals.f90 @@ -0,0 +1,26 @@ +! Singular Values +program example_svdvals + use stdlib_linalg_constants, only: dp + use stdlib_linalg, only: svdvals + implicit none + + real(dp), allocatable :: A(:,:),s(:) + character(*), parameter :: fmt="(a,*(1x,f12.8))" + + ! 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 fmt, ' ' + print fmt, 'S = ',s + print fmt, ' ' + +end program example_svdvals diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ced6e7b43..5e27927b1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -30,6 +30,7 @@ set(fppFiles stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp stdlib_linalg_state.fypp + stdlib_linalg_svd.fypp stdlib_optval.fypp stdlib_selection.fypp stdlib_sorting.fypp diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 134bf7af5..10cd16650 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -26,6 +26,8 @@ module stdlib_linalg public :: solve_lu public :: solve_lstsq public :: trace + public :: svd + public :: svdvals public :: outer_product public :: kronecker_product public :: cross_product @@ -552,6 +554,138 @@ module stdlib_linalg #:endfor end interface + ! Singular value decomposition + interface svd + !! 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-rank-2-array-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 left singular vectors + ${rt}$, optional, intent(out), target :: u(:,:) + !> 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 + !> [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 + !! + !! Computes the singular values of a `real` or `complex` 2d matrix. + !! ([Specification](../page/specs/stdlib_linalg.html#svdvals-compute-the-singular-values-of-a-rank-2-array-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 + 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_lapack.fypp b/src/stdlib_linalg_lapack.fypp index ac6f43f77..eb4c53c7d 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 new file mode 100644 index 000000000..b93343c15 --- /dev/null +++ b/src/stdlib_linalg_svd.fypp @@ -0,0 +1,293 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +submodule(stdlib_linalg) stdlib_linalg_svd +!! Singular-Value Decomposition + use stdlib_linalg_constants + 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) + + character(*), parameter :: this = 'svd' + + !> List of internal GESDD tasks: + + !> 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' + + contains + + !> Process GESDD output flag + elemental subroutine handle_gesdd_info(err,info,m,n) + !> Error handler + type(linalg_state_type), 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_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]) + case (-8) + 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]) + case (-4) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'A contains invalid/NaN values.') + case (1:) + err = linalg_state_type(this,LINALG_ERROR,'SVD computation did not converge.') + case default + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by GESDD.') + end select + + end subroutine handle_gesdd_info + + + #:for rk,rt,ri in RC_KINDS_TYPES + #:if rk!="xdp" + + !> 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 + type(linalg_state_type), 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 + 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 left singular vectors + ${rt}$, optional, intent(out), target :: u(:,:) + !> 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 + !> [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 + + !> Local variables + type(linalg_state_type) :: err0 + integer(ilp) :: m,n,lda,ldu,ldvt,info,k,lwork,liwork,lrwork + integer(ilp), allocatable :: iwork(:) + 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(:) + 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_type(this,LINALG_VALUE_ERROR,'invalid or matrix size: a=',[m,n]) + 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) + call linalg_error_handling(err0,err) + return + endif + + ! Integer storage + liwork = 8*k + allocate(iwork(liwork)) + + ! Can A be overwritten? By default, do not overwrite + overwrite_a_ = .false. + if (present(overwrite_a)) overwrite_a_ = overwrite_a + + ! Initialize a matrix temporary? + if (overwrite_a_) then + amat => a + else + allocate(amat(m,n),source=a) + endif + + ! 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 + 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 + ! User input + 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 + temp_u = .false. + elseif (.not.full_storage) then + ! Allocate with minimum size + allocate(umat(m,k)) + temp_u = .true. + else + ! Allocate with regular size + allocate(umat(m,m)) + temp_u = .true. + end if + + if (present(vt)) then + ! User input + vtmat => vt + temp_vt = .false. + elseif ((m vt_dummy + temp_vt = .false. + elseif (.not.full_storage) then + ! Allocate with minimum size + allocate(vtmat(k,n)) + temp_vt = .true. + else + ! Allocate with regular size + allocate(vtmat(n,n)) + temp_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_amat) 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)) + #: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) + + ! 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('complex')}#rwork,#{endif}#iwork,info) + call handle_gesdd_info(err0,info,m,n) + + endif + + ! Finalize storage and process output flag + 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}$ + #:endif + #:endfor + +end submodule stdlib_linalg_svd diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 9ec242213..eacfff53c 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -5,6 +5,7 @@ set( "test_linalg_solve.fypp" "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" + "test_linalg_svd.fypp" "test_linalg_matrix_property_checks.fypp" ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -14,4 +15,5 @@ ADDTEST(linalg_determinant) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_solve) 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..9fe8a889a --- /dev/null +++ b/test/linalg/test_linalg_svd.fypp @@ -0,0 +1,269 @@ +#: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,svd,svdvals + use stdlib_linalg_state, only: linalg_state_type + implicit none (type,external) + + public :: test_svd + + contains + + !> 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" + tests = [tests,new_unittest("test_svd_${ri}$",test_svd_${ri}$)] + #:endif + #:endfor + + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + tests = [tests,new_unittest("test_complex_svd_${ci}$",test_complex_svd_${ci}$)] + #:endif + #:endfor + + end subroutine test_svd + + !> Real matrix svd + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + subroutine test_svd_${ri}$(error) + type(error_type), allocatable, 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 + character(:), allocatable :: test + 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) + + 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) + + 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) + + 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) + + 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) + + 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) + + 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) + + 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) + + 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) + + 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}$ + + #:endif + #:endfor + + !> Test complex svd + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + subroutine test_complex_svd_${ci}$(error) + 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 :: 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([csqrt2,isqrt2,isqrt2,csqrt2],[2,2]) + ${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) + + !> Initialize matrix + A = A_mat + + !> Simple subroutine version + call svd(A,s,err=state) + + 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) + + 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) + + 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