linalg: Singular Value Decomposition (#808)

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!}
+```
+

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)

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

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

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

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
 
 

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