Merge branch 'master' into matrix_inverse

Author: perazz

doc/specs/stdlib_linalg.md

@@ -777,7 +777,8 @@ Result vector `x` returns the approximate solution that minimizes the 2-norm \(
 
 `cond` (optional): Shall be a scalar `real` value cut-off threshold for rank evaluation: `s_i >= cond*maxval(s), i=1:rank`. Shall be a scalar, `intent(in)` argument.
 
-`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `minval(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` argument.
+`singvals` (optional): Shall be a `real` rank-1 array of the same kind `a` and size at least `m
+al(shape(a))`, returning the list of singular values `s(i)>=cond*maxval(s)`, in descending order of magnitude. It is an `intent(out)` argument.
 
 `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `A` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument.
 
@@ -899,6 +900,96 @@ Exceptions trigger an `error stop`.
 {!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!}
+```
 ## `.inv.` - Inverse operator of a square matrix
 
 ### Status

doc/specs/stdlib_sorting.md

@@ -24,17 +24,18 @@ module's `string_type` type.
 
 ## Overview of the module
 
-The module `stdlib_sorting` defines several public entities, one
-default integer parameter, `int_index`, and four overloaded
+The module `stdlib_sorting` defines several public entities, two
+default integer parameters, `int_index` and `int_index_low`, and four overloaded
 subroutines: `ORD_SORT`, `SORT`, `RADIX_SORT` and `SORT_INDEX`. The
 overloaded subroutines also each have several specific names for
 versions corresponding to different types of array arguments.
 
-### The `int_index` parameter
+### The parameters `int_index` and `int_index_low`
 
-The `int_index` parameter is used to specify the kind of integer used
-in indexing the various arrays. Currently the module sets `int_index`
-to the value of `int64` from the `stdlib_kinds` module.
+The parameters `int_index` and `int_index_low` are used to specify the kind of integer used
+in indexing the various arrays. Currently the module sets `int_index` and
+`int_index_low`
+to the value of `int64` and `int32` from the `stdlib_kinds` module, respectively.
 
 ### The module subroutines
 
@@ -414,7 +415,7 @@ It is an `intent(inout)` argument. On input it
 will be an array whose sorting indices are to be determined. On return
 it will be the sorted array.
 
-`index`: shall be a rank one integer array of kind `int_index` and of
+`index`: shall be a rank one integer array of kind `int_index` or `int_index_low` and of
 the size of `array`. It is an `intent(out)` argument. On return it
 shall have values that are the indices needed to sort the original
 array in the desired direction.
@@ -426,8 +427,8 @@ memory for internal record keeping. If associated with an array in
 static storage, its use can significantly reduce the stack memory
 requirements for the code. Its contents on return are undefined.
 
-`iwork` (optional): shall be a rank one integer array of kind
-`int_index`, and shall have at least `size(array)/2` elements. It
+`iwork` (optional): shall be a rank one integer array of the same kind
+of the array `index`, and shall have at least `size(array)/2` elements. It
 is an `intent(out)` argument.  It is intended to be used as "scratch"
 memory for internal record keeping. If associated with an array in
 static storage, its use can significantly reduce the stack memory
@@ -457,6 +458,12 @@ different on return
 
 ##### Examples
 
+Sorting a rank one array with `sort_index`:
+
+```Fortran
+{!example/sorting/example_sort_index.f90!}
+```
+
 Sorting a related rank one array:
 
 ```Fortran
@@ -504,7 +511,7 @@ Sorting an array of a derived type based on the data in one component
 
 ```fortran
     subroutine sort_a_data( a_data, a, work, index, iwork )
-        ! Sort `a_data` in terms or its component `a`
+        ! Sort `a_data` in terms of its component `a`
         type(a_type), intent(inout)      :: a_data(:)
         integer(int32), intent(inout)    :: a(:)
         integer(int32), intent(out)    :: work(:)

example/linalg/CMakeLists.txt

@@ -26,5 +26,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

example/sorting/CMakeLists.txt

@@ -1,4 +1,5 @@
 ADD_EXAMPLE(ord_sort)
 ADD_EXAMPLE(sort)
+ADD_EXAMPLE(sort_index)
 ADD_EXAMPLE(radix_sort)
-ADD_EXAMPLE(sort_bitset)
+ADD_EXAMPLE(sort_bitset)

example/sorting/example_sort_index.f90

@@ -0,0 +1,15 @@
+program example_sort_index
+  use stdlib_sorting, only: sort_index
+  implicit none
+  integer, allocatable :: array(:)
+  integer, allocatable :: index(:)
+
+  array = [5, 4, 3, 1, 10, 4, 9]
+  allocate(index, mold=array)
+
+  call sort_index(array, index)
+
+  print *, array   !print [1, 3, 4, 4, 5, 9, 10]
+  print *, index   !print [4, 3, 2, 6, 1, 7, 5]
+
+end program example_sort_index

src/CMakeLists.txt

@@ -30,7 +30,8 @@ set(fppFiles
     stdlib_linalg_solve.fypp
     stdlib_linalg_determinant.fypp
     stdlib_linalg_inverse.fypp
-    stdlib_linalg_state.fypp 
+    stdlib_linalg_state.fypp
+    stdlib_linalg_svd.fypp 
     stdlib_optval.fypp
     stdlib_selection.fypp
     stdlib_sorting.fypp