Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Sign up
Appearance settings

Commit c0c0647

Browse files
authored
linalg: Singular Value Decomposition (#808)
2 parents efb53f5 + 9f51efe commit c0c0647

File tree

10 files changed

+878
-9
lines changed

10 files changed

+878
-9
lines changed

‎doc/specs/stdlib_linalg.md

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -898,3 +898,95 @@ Exceptions trigger an `error stop`.
898898
```fortran
899899
{!example/linalg/example_determinant2.f90!}
900900
```
901+
902+
## `svd` - Compute the singular value decomposition of a rank-2 array (matrix).
903+
904+
### Status
905+
906+
Experimental
907+
908+
### Description
909+
910+
This subroutine computes the singular value decomposition of a `real` or `complex` rank-2 array (matrix) \( A = U \cdot S \cdot \V^T \).
911+
The solver is based on LAPACK's `*GESDD` backends.
912+
913+
Result vector `s` returns the array of singular values on the diagonal of \( S \).
914+
If requested, `u` contains the left singular vectors, as columns of \( U \).
915+
If requested, `vt` contains the right singular vectors, as rows of \( V^T \).
916+
917+
### Syntax
918+
919+
`call ` [[stdlib_linalg(module):svd(interface)]] `(a, s, [, u, vt, overwrite_a, full_matrices, err])`
920+
921+
### Class
922+
Subroutine
923+
924+
### Arguments
925+
926+
`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.`.
927+
928+
`s`: Shall be a rank-1 `real` array, returning the list of `k = min(m,n)` singular values. It is an `intent(out)` argument.
929+
930+
`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.
931+
932+
`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.
933+
934+
`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.
935+
936+
`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.
937+
938+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
939+
940+
### Return values
941+
942+
Returns an array `s` that contains the list of singular values of matrix `a`.
943+
If requested, returns a rank-2 array `u` that contains the left singular vectors of `a` along its columns.
944+
If requested, returns a rank-2 array `vt` that contains the right singular vectors of `a` along its rows.
945+
946+
Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
947+
Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes.
948+
Exceptions trigger an `error stop`, unless argument `err` is present.
949+
950+
### Example
951+
952+
```fortran
953+
{!example/linalg/example_svd.f90!}
954+
```
955+
956+
## `svdvals` - Compute the singular values of a rank-2 array (matrix).
957+
958+
### Status
959+
960+
Experimental
961+
962+
### Description
963+
964+
This subroutine computes the singular values of a `real` or `complex` rank-2 array (matrix) from its singular
965+
value decomposition \( A = U \cdot S \cdot \V^T \). The solver is based on LAPACK's `*GESDD` backends.
966+
967+
Result vector `s` returns the array of singular values on the diagonal of \( S \).
968+
969+
### Syntax
970+
971+
`s = ` [[stdlib_linalg(module):svdvals(interface)]] `(a [, err])`
972+
973+
### Arguments
974+
975+
`a`: Shall be a rank-2 `real` or `complex` array containing the coefficient matrix of size `[m,n]`. It is an `intent(in)` argument.
976+
977+
`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument.
978+
979+
### Return values
980+
981+
Returns an array `s` that contains the list of singular values of matrix `a`.
982+
983+
Raises `LINALG_ERROR` if the underlying Singular Value Decomposition process did not converge.
984+
Raises `LINALG_VALUE_ERROR` if the matrix or any of the output arrays invalid/incompatible sizes.
985+
Exceptions trigger an `error stop`, unless argument `err` is present.
986+
987+
### Example
988+
989+
```fortran
990+
{!example/linalg/example_svdvals.f90!}
991+
```
992+

‎example/linalg/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,5 +23,7 @@ ADD_EXAMPLE(lstsq2)
2323
ADD_EXAMPLE(solve1)
2424
ADD_EXAMPLE(solve2)
2525
ADD_EXAMPLE(solve3)
26+
ADD_EXAMPLE(svd)
27+
ADD_EXAMPLE(svdvals)
2628
ADD_EXAMPLE(determinant)
2729
ADD_EXAMPLE(determinant2)

‎example/linalg/example_svd.f90

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
! Singular Value Decomposition
2+
program example_svd
3+
use stdlib_linalg_constants, only: dp
4+
use stdlib_linalg, only: svd
5+
implicit none
6+
7+
real(dp), allocatable :: A(:,:),s(:),u(:,:),vt(:,:)
8+
character(*), parameter :: fmt = "(a,*(1x,f12.8))"
9+
10+
! We want to find the singular value decomposition of matrix:
11+
!
12+
! A = [ 3 2 2]
13+
! [ 2 3 -2]
14+
!
15+
A = transpose(reshape([ 3, 2, 2, &
16+
2, 3,-2], [3,2]))
17+
18+
! Prepare arrays
19+
allocate(s(2),u(2,2),vt(3,3))
20+
21+
! Get singular value decomposition
22+
call svd(A,s,u,vt)
23+
24+
! Singular values: [5, 3]
25+
print fmt, ' '
26+
print fmt, 'S = ',s
27+
print fmt, ' '
28+
29+
! Left vectors (may be flipped):
30+
! [Ã2/2 Ã2/2]
31+
! U = [Ã2/2 -Ã2/2]
32+
!
33+
print fmt, ' '
34+
print fmt, 'U = ',u(1,:)
35+
print fmt, ' ',u(2,:)
36+
37+
38+
! Right vectors (may be flipped):
39+
! [Ã2/2 Ã2/2 0]
40+
! V = [1/Ã18 -1/Ã18 4/Ã18]
41+
! [ 2/3 -2/3 -1/3]
42+
!
43+
print fmt, ' '
44+
print fmt, ' ',vt(1,:)
45+
print fmt, 'VT= ',vt(2,:)
46+
print fmt, ' ',vt(3,:)
47+
print fmt, ' '
48+
49+
50+
end program example_svd

‎example/linalg/example_svdvals.f90

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
! Singular Values
2+
program example_svdvals
3+
use stdlib_linalg_constants, only: dp
4+
use stdlib_linalg, only: svdvals
5+
implicit none
6+
7+
real(dp), allocatable :: A(:,:),s(:)
8+
character(*), parameter :: fmt="(a,*(1x,f12.8))"
9+
10+
! We want to find the singular values of matrix:
11+
!
12+
! A = [ 3 2 2]
13+
! [ 2 3 -2]
14+
!
15+
A = transpose(reshape([ 3, 2, 2, &
16+
2, 3,-2], [3,2]))
17+
18+
! Get singular values
19+
s = svdvals(A)
20+
21+
! Singular values: [5, 3]
22+
print fmt, ' '
23+
print fmt, 'S = ',s
24+
print fmt, ' '
25+
26+
end program example_svdvals

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ set(fppFiles
3030
stdlib_linalg_solve.fypp
3131
stdlib_linalg_determinant.fypp
3232
stdlib_linalg_state.fypp
33+
stdlib_linalg_svd.fypp
3334
stdlib_optval.fypp
3435
stdlib_selection.fypp
3536
stdlib_sorting.fypp

‎src/stdlib_linalg.fypp

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@ module stdlib_linalg
2626
public :: solve_lu
2727
public :: solve_lstsq
2828
public :: trace
29+
public :: svd
30+
public :: svdvals
2931
public :: outer_product
3032
public :: kronecker_product
3133
public :: cross_product
@@ -552,6 +554,138 @@ module stdlib_linalg
552554
#:endfor
553555
end interface
554556

557+
! Singular value decomposition
558+
interface svd
559+
!! version: experimental
560+
!!
561+
!! Computes the singular value decomposition of a `real` or `complex` 2d matrix.
562+
!! ([Specification](../page/specs/stdlib_linalg.html#svd-compute-the-singular-value-decomposition-of-a-rank-2-array-matrix))
563+
!!
564+
!!### Summary
565+
!! Interface for computing the singular value decomposition of a `real` or `complex` 2d matrix.
566+
!!
567+
!!### Description
568+
!!
569+
!! This interface provides methods for computing the singular value decomposition of a matrix.
570+
!! Supported data types include `real` and `complex`. The subroutine returns a `real` array of
571+
!! singular values, and optionally, left- and right- singular vector matrices, `U` and `V`.
572+
!! For a matrix `A` with size [m,n], full matrix storage for `U` and `V` should be [m,m] and [n,n].
573+
!! It is possible to use partial storage [m,k] and [k,n], `k=min(m,n)`, choosing `full_matrices=.false.`.
574+
!!
575+
!!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
576+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
577+
!!
578+
!!### Example
579+
!!
580+
!!```fortran
581+
!! real(sp) :: a(2,3), s(2), u(2,2), vt(3,3)
582+
!! a = reshape([3,2, 2,3, 2,-2],[2,3])
583+
!!
584+
!! call svd(A,s,u,v)
585+
!! print *, 'singular values = ',s
586+
!!```
587+
!!
588+
#:for rk,rt,ri in RC_KINDS_TYPES
589+
#:if rk!="xdp"
590+
module subroutine stdlib_linalg_svd_${ri}$(a,s,u,vt,overwrite_a,full_matrices,err)
591+
!!### Summary
592+
!! Compute singular value decomposition of a matrix \( A = U \cdot S \cdot \V^T \)
593+
!!
594+
!!### Description
595+
!!
596+
!! This function computes the singular value decomposition of a `real` or `complex` matrix \( A \),
597+
!! and returns the array of singular values, and optionally the left matrix \( U \) containing the
598+
!! left unitary singular vectors, and the right matrix \( V^T \), containing the right unitary
599+
!! singular vectors.
600+
!!
601+
!! param: a Input matrix of size [m,n].
602+
!! param: s Output `real` array of size [min(m,n)] returning a list of singular values.
603+
!! param: u [optional] Output left singular matrix of size [m,m] or [m,min(m,n)] (.not.full_matrices). Contains singular vectors as columns.
604+
!! param: vt [optional] Output right singular matrix of size [n,n] or [min(m,n),n] (.not.full_matrices). Contains singular vectors as rows.
605+
!! param: overwrite_a [optional] Flag indicating if the input matrix can be overwritten.
606+
!! 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)`.
607+
!! param: err [optional] State return flag.
608+
!!
609+
!> Input matrix A[m,n]
610+
${rt},ドル intent(inout), target :: a(:,:)
611+
!> Array of singular values
612+
real(${rk}$), intent(out) :: s(:)
613+
!> The columns of U contain the left singular vectors
614+
${rt},ドル optional, intent(out), target :: u(:,:)
615+
!> The rows of V^T contain the right singular vectors
616+
${rt},ドル optional, intent(out), target :: vt(:,:)
617+
!> [optional] Can A data be overwritten and destroyed?
618+
logical(lk), optional, intent(in) :: overwrite_a
619+
!> [optional] full matrices have shape(u)==[m,m], shape(vh)==[n,n] (default); otherwise
620+
!> they are shape(u)==[m,k] and shape(vh)==[k,n] with k=min(m,n)
621+
logical(lk), optional, intent(in) :: full_matrices
622+
!> [optional] state return flag. On error if not requested, the code will stop
623+
type(linalg_state_type), optional, intent(out) :: err
624+
end subroutine stdlib_linalg_svd_${ri}$
625+
#:endif
626+
#:endfor
627+
end interface svd
628+
629+
! Singular values
630+
interface svdvals
631+
!! version: experimental
632+
!!
633+
!! Computes the singular values of a `real` or `complex` 2d matrix.
634+
!! ([Specification](../page/specs/stdlib_linalg.html#svdvals-compute-the-singular-values-of-a-rank-2-array-matrix))
635+
!!
636+
!!### Summary
637+
!!
638+
!! Function interface for computing the array of singular values from the singular value decomposition
639+
!! of a `real` or `complex` 2d matrix.
640+
!!
641+
!!### Description
642+
!!
643+
!! This interface provides methods for computing the singular values a 2d matrix.
644+
!! Supported data types include `real` and `complex`. The function returns a `real` array of
645+
!! singular values, with size [min(m,n)].
646+
!!
647+
!!@note The solution is based on LAPACK's singular value decomposition `*GESDD` methods.
648+
!!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``).
649+
!!
650+
!!### Example
651+
!!
652+
!!```fortran
653+
!! real(sp) :: a(2,3), s(2)
654+
!! a = reshape([3,2, 2,3, 2,-2],[2,3])
655+
!!
656+
!! s = svdvals(A)
657+
!! print *, 'singular values = ',s
658+
!!```
659+
!!
660+
#:for rk,rt,ri in RC_KINDS_TYPES
661+
#:if rk!="xdp"
662+
module function stdlib_linalg_svdvals_${ri}$(a,err) result(s)
663+
!!### Summary
664+
!! Compute singular values \(S \) from the singular-value decomposition of a matrix \( A = U \cdot S \cdot \V^T \).
665+
!!
666+
!!### Description
667+
!!
668+
!! This function returns the array of singular values from the singular value decomposition of a `real`
669+
!! or `complex` matrix \( A = U \cdot S \cdot V^T \).
670+
!!
671+
!! param: a Input matrix of size [m,n].
672+
!! param: err [optional] State return flag.
673+
!!
674+
!!### Return value
675+
!!
676+
!! param: s `real` array of size [min(m,n)] returning a list of singular values.
677+
!!
678+
!> Input matrix A[m,n]
679+
${rt},ドル intent(in), target :: a(:,:)
680+
!> [optional] state return flag. On error if not requested, the code will stop
681+
type(linalg_state_type), optional, intent(out) :: err
682+
!> Array of singular values
683+
real(${rk}$), allocatable :: s(:)
684+
end function stdlib_linalg_svdvals_${ri}$
685+
#:endif
686+
#:endfor
687+
end interface svdvals
688+
555689
contains
556690

557691

‎src/stdlib_linalg_lapack.fypp

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -19,16 +19,16 @@ module stdlib_linalg_lapack
1919
interface bbcsd
2020
!! BBCSD computes the CS decomposition of a unitary matrix in
2121
!! bidiagonal-block form,
22-
!! [ B11 | B12 0 0 ]
23-
!! [ 0 | 0 -I 0 ]
22+
!! [ B11 | B12 0 0 ]
23+
!! [ 0 | 0 -I 0 ]
2424
!! X = [----------------]
25-
!! [ B21 | B22 0 0 ]
26-
!! [ 0 | 0 0 I ]
27-
!! [ C | -S 0 0 ]
28-
!! [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H
29-
!! = [---------] [---------------] [---------] .
30-
!! [ | U2 ] [ S | C 0 0 ] [ | V2 ]
31-
!! [ 0 | 0 0 I ]
25+
!! [ B21 | B22 0 0 ]
26+
!! [ 0 | 0 0 I ]
27+
!! [ C | -S 0 0 ]
28+
!! [ U1 | ] [ 0 | 0 -I 0 ] [ V1 | ]**H
29+
!! = [---------] [---------------] [---------] .
30+
!! [ | U2 ] [ S | C 0 0 ] [ | V2 ]
31+
!! [ 0 | 0 0 I ]
3232
!! X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
3333
!! than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
3434
!! transposed and/or permuted. This can be done in constant time using

0 commit comments

Comments
(0)

AltStyle によって変換されたページ (->オリジナル) /