Skip to content

Commit 94f201d

Browse files
perazzperazz
authored
linalg: Schur decomposition (#892)
2 parents 191879f + f9a31cd commit 94f201d

10 files changed

+854
-0
lines changed

‎doc/specs/stdlib_linalg.md

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1007,6 +1007,82 @@ This subroutine computes the internal working space requirements for the QR fact
10071007
{!example/linalg/example_qr_space.f90!}
10081008
```
10091009

1010+
## `schur` - Compute the Schur decomposition of a matrix
1011+
1012+
### Status
1013+
1014+
Experimental
1015+
1016+
### Description
1017+
1018+
This subroutine computes the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is unitary (orthonormal) and \( T \) is upper-triangular (for complex matrices) or quasi-upper-triangular (for real matrices, with possible \( 2 \times 2 \) blocks on the diagonal). Matrix \( A \) has size `[n,n]`.
1019+
1020+
The results are returned in output matrices \( T \) and \( Z \). Matrix \( T \) is the Schur form, and matrix \( Z \) is the unitary transformation matrix such that \( A = Z T Z^H \). If requested, the eigenvalues of \( T \) can also be returned as a `complex` array of size `[n]`.
1021+
1022+
### Syntax
1023+
1024+
`call ` [[stdlib_linalg(module):schur(interface)]] `(a, t [, z,] [, eigvals] [, overwrite_a] [, storage] [, err])`
1025+
1026+
### Arguments
1027+
1028+
- `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(inout)` argument if `overwrite_a = .true.`; otherwise, it is an `intent(in)` argument.
1029+
1030+
- `t`: Shall be a rank-2 array of the same kind as `a`, containing the Schur form \( T \) of the matrix. It is an `intent(out)` argument and should have a shape of `[n,n]`.
1031+
1032+
- `z`: Shall be a rank-2 array of the same kind as `a`, containing the unitary matrix \( Z \). It is an `intent(out)` argument and is optional. If provided, it should have the shape `[n,n]`.
1033+
1034+
- `eigvals` (optional): Shall be a rank-1 `complex` or `real` array of the same kind as `a`, containing the eigenvalues of \( A \) (the diagonal elements of \( T \)), or their `real` component only. The array must be of size `[n]`. If not provided, the eigenvalues are not returned. It is an `intent(out)` argument.
1035+
1036+
- `overwrite_a` (optional): Shall be a `logical` flag (default: `.false.`). If `.true.`, the input matrix `a` will be overwritten and destroyed upon return, avoiding internal data allocation. It is an `intent(in)` argument.
1037+
1038+
- `storage` (optional): Shall be a rank-1 array of the same type and kind as `a`, providing working storage for the solver. Its minimum size can be determined with a call to [[stdlib_linalg(module):schur_space(interface)]]. It is an `intent(inout)` argument.
1039+
1040+
- `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. If not provided, exceptions trigger an `error stop`.
1041+
1042+
### Return value
1043+
1044+
Returns the Schur decomposition matrices into the \( T \) and \( Z \) arguments. If `eigvals` is provided, it will also return the eigenvalues of the matrix \( A \).
1045+
1046+
Raises `LINALG_VALUE_ERROR` if any of the matrices have invalid or unsuitable size for the decomposition.
1047+
Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts.
1048+
Raises `LINALG_ERROR` on insufficient user storage space.
1049+
If the state argument `err` is not present, exceptions trigger an `error stop`.
1050+
1051+
### Example
1052+
1053+
```fortran
1054+
{!example/linalg/example_schur_eigvals.f90!}
1055+
```
1056+
1057+
---
1058+
1059+
## `schur_space` - Compute internal working space requirements for the Schur decomposition
1060+
1061+
### Status
1062+
1063+
Experimental
1064+
1065+
### Description
1066+
1067+
This subroutine computes the internal working space requirements for the Schur decomposition, [[stdlib_linalg(module):schur(interface)]].
1068+
1069+
### Syntax
1070+
1071+
`call ` [[stdlib_linalg(module):schur_space(interface)]] `(a, lwork, [, err])`
1072+
1073+
### Arguments
1074+
1075+
- `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(in)` argument.
1076+
1077+
- `lwork`: Shall be an `integer` scalar that returns the minimum array size required for the working storage in [[stdlib_linalg(module):schur(interface)]] to decompose `a`. It is an `intent(out)` argument.
1078+
1079+
- `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
1080+
1081+
### Return value
1082+
1083+
Returns the required working storage size `lwork` for the Schur decomposition. This value can be used to pre-allocate a workspace array in case multiple Schur decompositions of the same matrix size are needed. If pre-allocated working arrays are provided, no internal memory allocations will take place during the decomposition.
1084+
1085+
10101086
## `eig` - Eigenvalues and Eigenvectors of a Square Matrix
10111087

10121088
### Status

‎example/linalg/CMakeLists.txt

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,9 @@ ADD_EXAMPLE(eigvalsh)
2626
ADD_EXAMPLE(trace)
2727
ADD_EXAMPLE(state1)
2828
ADD_EXAMPLE(state2)
29+
ADD_EXAMPLE(schur_real)
30+
ADD_EXAMPLE(schur_complex)
31+
ADD_EXAMPLE(schur_eigvals)
2932
ADD_EXAMPLE(blas_gemv)
3033
ADD_EXAMPLE(lapack_getrf)
3134
ADD_EXAMPLE(lstsq1)

‎example/linalg/example_schur_complex.f90

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
! This example demonstrates the Schur decomposition for a complex-valued matrix.
2+
program example_schur_complex
3+
use stdlib_linalg, only: schur
4+
use stdlib_linalg_constants, only: dp
5+
implicit none
6+
7+
integer, parameter :: n = 3
8+
complex(dp), dimension(n,n) :: A, T, Z
9+
10+
! Initialize a complex-valued square matrix
11+
A = reshape([ (1, 2), (3,-1), (4, 1), &
12+
(0,-1), (2, 0), (1,-2), &
13+
(2, 3), (1, 1), (0,-1) ], shape=[n,n])
14+
15+
! Compute the Schur decomposition: A = Z T Z^H
16+
call schur(A, T, Z)
17+
18+
! Output results
19+
print *, "Original Matrix A:"
20+
print *, A
21+
print *, "Schur Form Matrix T:"
22+
print *, T
23+
print *, "Unitary Matrix Z:"
24+
print *, Z
25+
26+
! Test factorization: Z*T*Z^H = A
27+
print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, conjg(transpose(Z)))) - A))
28+
29+
end program example_schur_complex
30+

‎example/linalg/example_schur_eigvals.f90

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
! This example includes eigenvalue computation in addition to
2+
! the Schur decomposition for a randomly generated matrix.
3+
program example_schur_eigenvalues
4+
use stdlib_linalg, only: schur
5+
use stdlib_linalg_constants, only: dp
6+
implicit none
7+
8+
integer, parameter :: n = 5
9+
real(dp), dimension(n,n) :: A, T, Z
10+
complex(dp), dimension(n) :: eigenvalues
11+
12+
! Create a random real-valued square matrix
13+
call random_number(A)
14+
15+
! Compute the Schur decomposition and eigenvalues
16+
call schur(A, T, Z, eigenvalues)
17+
18+
! Output results
19+
print *, "Random Matrix A:"
20+
print *, A
21+
print *, "Schur Form Matrix T:"
22+
print *, T
23+
print *, "Orthogonal Matrix Z:"
24+
print *, Z
25+
print *, "Eigenvalues:"
26+
print *, eigenvalues
27+
28+
! Test factorization: Z*T*Z^T = A
29+
print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
30+
31+
end program example_schur_eigenvalues
32+

‎example/linalg/example_schur_real.f90

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
! This example computes the Schur decomposition of a real-valued square matrix.
2+
program example_schur_real
3+
use stdlib_linalg, only: schur
4+
use stdlib_linalg_constants, only: dp
5+
implicit none
6+
integer, parameter :: n = 3
7+
real(dp), dimension(n,n) :: A, T, Z
8+
9+
! Initialize a real-valued square matrix
10+
A = reshape([ 0, 2, 2, &
11+
0, 1, 2, &
12+
1, 0, 1], shape=[n,n])
13+
14+
! Compute the Schur decomposition: A = Z T Z^T
15+
call schur(A, T, Z)
16+
17+
! Output results
18+
print *, "Original Matrix A:"
19+
print *, A
20+
print *, "Schur Form Matrix T:"
21+
print *, T
22+
print *, "Orthogonal Matrix Z:"
23+
print *, Z
24+
25+
! Test factorization: Z*T*Z^T = A
26+
print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
27+
28+
end program example_schur_real
29+

‎src/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@ set(fppFiles
3737
stdlib_linalg_state.fypp
3838
stdlib_linalg_svd.fypp
3939
stdlib_linalg_cholesky.fypp
40+
stdlib_linalg_schur.fypp
4041
stdlib_optval.fypp
4142
stdlib_selection.fypp
4243
stdlib_sorting.fypp

‎src/stdlib_linalg.fypp

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,8 @@ module stdlib_linalg
4848
public :: cross_product
4949
public :: qr
5050
public :: qr_space
51+
public :: schur
52+
public :: schur_space
5153
public :: is_square
5254
public :: is_diagonal
5355
public :: is_symmetric
@@ -638,6 +640,92 @@ module stdlib_linalg
638640
#:endfor
639641
end interface qr_space
640642

643+
! Schur decomposition of rank-2 array A
644+
interface schur
645+
!! version: experimental
646+
!!
647+
!! Computes the Schur decomposition of matrix \( A = Z T Z^H \).
648+
!! ([Specification](../page/specs/stdlib_linalg.html#schur-compute-the-schur-decomposition-of-a-matrix))
649+
!!
650+
!!### Summary
651+
!! Compute the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is
652+
!! orthonormal/unitary and \( T \) is upper-triangular or quasi-upper-triangular. Matrix \( A \) has size `[m,m]`.
653+
!!
654+
!!### Description
655+
!!
656+
!! This interface provides methods for computing the Schur decomposition of a matrix.
657+
!! Supported data types include `real` and `complex`. If a pre-allocated workspace is provided, no internal
658+
!! memory allocations take place when using this interface.
659+
!!
660+
!! The output matrix \( T \) is upper-triangular for `complex` input matrices and quasi-upper-triangular
661+
!! for `real` input matrices (with possible \( 2 \times 2 \) blocks on the diagonal).
662+
!!
663+
!!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`).
664+
!!
665+
#:for rk,rt,ri in RC_KINDS_TYPES
666+
module subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, overwrite_a, storage, err)
667+
!> Input matrix a[m,m]
668+
${rt}$, intent(inout), target :: a(:,:)
669+
!> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
670+
${rt}$, intent(out), contiguous, target :: t(:,:)
671+
!> Unitary/orthonormal transformation matrix Z
672+
${rt}$, optional, intent(out), contiguous, target :: z(:,:)
673+
!> [optional] Output eigenvalues that appear on the diagonal of T
674+
complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
675+
!> [optional] Can A data be overwritten and destroyed?
676+
logical(lk), optional, intent(in) :: overwrite_a
677+
!> [optional] Provide pre-allocated workspace, size to be checked with schur_space
678+
${rt}$, optional, intent(inout), target :: storage(:)
679+
!> [optional] State return flag. On error if not requested, the code will stop
680+
type(linalg_state_type), optional, intent(out) :: err
681+
end subroutine stdlib_linalg_${ri}$_schur
682+
683+
! Schur decomposition subroutine: real eigenvalue interface
684+
module subroutine stdlib_linalg_real_eig_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
685+
!> Input matrix a[m,m]
686+
${rt}$, intent(inout), target :: a(:,:)
687+
!> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
688+
${rt}$, intent(out), contiguous, target :: t(:,:)
689+
!> Unitary/orthonormal transformation matrix Z
690+
${rt}$, optional, intent(out), contiguous, target :: z(:,:)
691+
!> Output real eigenvalues that appear on the diagonal of T
692+
real(${rk}$), intent(out), contiguous, target :: eigvals(:)
693+
!> [optional] Provide pre-allocated workspace, size to be checked with schur_space
694+
${rt}$, optional, intent(inout), target :: storage(:)
695+
!> [optional] Can A data be overwritten and destroyed?
696+
logical(lk), optional, intent(in) :: overwrite_a
697+
!> [optional] State return flag. On error if not requested, the code will stop
698+
type(linalg_state_type), optional, intent(out) :: err
699+
end subroutine stdlib_linalg_real_eig_${ri}$_schur
700+
#:endfor
701+
end interface schur
702+
703+
! Return the working array space required by the Schur decomposition solver
704+
interface schur_space
705+
!! version: experimental
706+
!!
707+
!! Computes the working array space required by the Schur decomposition solver
708+
!! ([Specification](../page/specs/stdlib_linalg.html#schur-space-compute-internal-working-space-requirements-for-the-schur-decomposition))
709+
!!
710+
!!### Description
711+
!!
712+
!! This interface returns the size of the `real` or `complex` working storage required by the
713+
!! Schur decomposition solver. The working size only depends on the kind (`real` or `complex`) and size of
714+
!! the matrix being decomposed. Storage size can be used to pre-allocate a working array in case several
715+
!! repeated Schur decompositions of same-size matrices are sought. If pre-allocated working arrays
716+
!! are provided, no internal allocations will take place during the decomposition.
717+
!!
718+
#:for rk,rt,ri in RC_KINDS_TYPES
719+
module subroutine get_schur_${ri}$_workspace(a,lwork,err)
720+
!> Input matrix a[m,m]
721+
${rt}$, intent(in), target :: a(:,:)
722+
!> Minimum workspace size for the decomposition operation
723+
integer(ilp), intent(out) :: lwork
724+
!> State return flag. Returns an error if the query failed
725+
type(linalg_state_type), optional, intent(out) :: err
726+
end subroutine get_schur_${ri}$_workspace
727+
#:endfor
728+
end interface schur_space
641729

642730
interface det
643731
!! version: experimental

0 commit comments

Comments
 (0)