linalg: Schur decomposition

doc/specs/stdlib_linalg.md

@@ -1007,6 +1007,82 @@ This subroutine computes the internal working space requirements for the QR fact
 {!example/linalg/example_qr_space.f90!}
 ```
 
+## `schur` - Compute the Schur decomposition of a matrix
+
+### Status
+
+Experimental
+
+### Description
+
+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]`.
+
+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]`.
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):schur(interface)]] `(a, t [, z,] [, eigvals] [, overwrite_a] [, storage] [, err])`
+
+### Arguments
+
+- `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.
+
+- `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]`.
+
+- `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]`.
+
+- `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.
+
+- `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.
+
+- `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.
+
+- `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument. If not provided, exceptions trigger an `error stop`.
+
+### Return value
+
+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 \).
+
+Raises `LINALG_VALUE_ERROR` if any of the matrices have invalid or unsuitable size for the decomposition.  
+Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts.
+Raises `LINALG_ERROR` on insufficient user storage space.  
+If the state argument `err` is not present, exceptions trigger an `error stop`.
+
+### Example
+
+```fortran
+{!example/linalg/example_schur_eigvals.f90!}
+```
+
+---
+
+## `schur_space` - Compute internal working space requirements for the Schur decomposition
+
+### Status
+
+Experimental
+
+### Description
+
+This subroutine computes the internal working space requirements for the Schur decomposition, [[stdlib_linalg(module):schur(interface)]].
+
+### Syntax
+
+`call ` [[stdlib_linalg(module):schur_space(interface)]] `(a, lwork, [, err])`
+
+### Arguments
+
+- `a`: Shall be a rank-2 `real` or `complex` array containing the matrix to be decomposed. It is an `intent(in)` argument.
+
+- `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.
+
+- `err` (optional): Shall be a `type(linalg_state_type)` value. It is an `intent(out)` argument.
+
+### Return value
+
+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.
+
+
 ## `eig` - Eigenvalues and Eigenvectors of a Square Matrix
 
 ### Status

example/linalg/CMakeLists.txt

@@ -26,6 +26,9 @@ ADD_EXAMPLE(eigvalsh)
 ADD_EXAMPLE(trace)
 ADD_EXAMPLE(state1)
 ADD_EXAMPLE(state2)
+ADD_EXAMPLE(schur_real)
+ADD_EXAMPLE(schur_complex)
+ADD_EXAMPLE(schur_eigvals)
 ADD_EXAMPLE(blas_gemv)
 ADD_EXAMPLE(lapack_getrf)
 ADD_EXAMPLE(lstsq1)

example/linalg/example_schur_complex.f90

@@ -0,0 +1,30 @@
+! This example demonstrates the Schur decomposition for a complex-valued matrix.
+program example_schur_complex
+  use stdlib_linalg, only: schur
+  use stdlib_linalg_constants, only: dp
+  implicit none
+
+  integer, parameter :: n = 3
+  complex(dp), dimension(n,n) :: A, T, Z
+
+  ! Initialize a complex-valued square matrix
+  A = reshape([ (1, 2), (3,-1), (4, 1), &
+                (0,-1), (2, 0), (1,-2), &
+                (2, 3), (1, 1), (0,-1) ], shape=[n,n])
+
+  ! Compute the Schur decomposition: A = Z T Z^H
+  call schur(A, T, Z)
+
+  ! Output results
+  print *, "Original Matrix A:"
+  print *, A
+  print *, "Schur Form Matrix T:"
+  print *, T
+  print *, "Unitary Matrix Z:"
+  print *, Z
+
+  ! Test factorization: Z*T*Z^H = A
+  print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, conjg(transpose(Z)))) - A))
+
+end program example_schur_complex
+

example/linalg/example_schur_eigvals.f90

@@ -0,0 +1,32 @@
+! This example includes eigenvalue computation in addition to 
+! the Schur decomposition for a randomly generated matrix.
+program example_schur_eigenvalues
+  use stdlib_linalg, only: schur
+  use stdlib_linalg_constants, only: dp
+  implicit none
+
+  integer, parameter :: n = 5  
+  real(dp), dimension(n,n) :: A, T, Z
+  complex(dp), dimension(n) :: eigenvalues
+
+  ! Create a random real-valued square matrix
+  call random_number(A)
+
+  ! Compute the Schur decomposition and eigenvalues
+  call schur(A, T, Z, eigenvalues)
+
+  ! Output results
+  print *, "Random Matrix A:"
+  print *, A
+  print *, "Schur Form Matrix T:"
+  print *, T
+  print *, "Orthogonal Matrix Z:"
+  print *, Z
+  print *, "Eigenvalues:"
+  print *, eigenvalues
+
+  ! Test factorization: Z*T*Z^T = A
+  print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
+
+end program example_schur_eigenvalues
+

example/linalg/example_schur_real.f90

@@ -0,0 +1,29 @@
+! This example computes the Schur decomposition of a real-valued square matrix.
+program example_schur_real
+  use stdlib_linalg, only: schur
+  use stdlib_linalg_constants, only: dp
+  implicit none
+  integer, parameter :: n = 3
+  real(dp), dimension(n,n) :: A, T, Z
+
+  ! Initialize a real-valued square matrix
+  A = reshape([ 0, 2, 2, &
+                0, 1, 2, &
+                1, 0, 1], shape=[n,n])
+
+  ! Compute the Schur decomposition: A = Z T Z^T
+  call schur(A, T, Z)
+
+  ! Output results
+  print *, "Original Matrix A:"
+  print *, A
+  print *, "Schur Form Matrix T:"
+  print *, T
+  print *, "Orthogonal Matrix Z:"
+  print *, Z
+
+  ! Test factorization: Z*T*Z^T = A
+  print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
+
+end program example_schur_real
+

src/CMakeLists.txt

@@ -37,6 +37,7 @@ set(fppFiles
     stdlib_linalg_state.fypp
     stdlib_linalg_svd.fypp 
     stdlib_linalg_cholesky.fypp
+    stdlib_linalg_schur.fypp
     stdlib_optval.fypp
     stdlib_selection.fypp
     stdlib_sorting.fypp

src/stdlib_linalg.fypp

@@ -48,6 +48,8 @@ module stdlib_linalg
   public :: cross_product
   public :: qr
   public :: qr_space
+  public :: schur
+  public :: schur_space
   public :: is_square
   public :: is_diagonal
   public :: is_symmetric
@@ -638,6 +640,92 @@ module stdlib_linalg
     #:endfor
   end interface qr_space
 
+  ! Schur decomposition of rank-2 array A
+  interface schur
+    !! version: experimental
+    !!
+    !! Computes the Schur decomposition of matrix \( A = Z T Z^H \).
+    !! ([Specification](../page/specs/stdlib_linalg.html#schur-compute-the-schur-decomposition-of-a-matrix))
+    !!
+    !!### Summary
+    !! Compute the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is
+    !! orthonormal/unitary and \( T \) is upper-triangular or quasi-upper-triangular. Matrix \( A \) has size `[m,m]`.
+    !!
+    !!### Description
+    !! 
+    !! This interface provides methods for computing the Schur decomposition of a matrix.
+    !! Supported data types include `real` and `complex`. If a pre-allocated workspace is provided, no internal 
+    !! memory allocations take place when using this interface.
+    !!
+    !! The output matrix \( T \) is upper-triangular for `complex` input matrices and quasi-upper-triangular
+    !! for `real` input matrices (with possible \( 2 \times 2 \) blocks on the diagonal).
+    !!
+    !!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`). 
+    !! 
+    #:for rk,rt,ri in RC_KINDS_TYPES
+      module subroutine stdlib_linalg_${ri}$_schur(a, t, z, eigvals, overwrite_a, storage, err)
+         !> Input matrix a[m,m]
+         ${rt}$, intent(inout), target :: a(:,:)
+         !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
+         ${rt}$, intent(out), contiguous, target :: t(:,:)
+         !> Unitary/orthonormal transformation matrix Z
+         ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
+         !> [optional] Output eigenvalues that appear on the diagonal of T
+         complex(${rk}$), optional, intent(out), contiguous, target :: eigvals(:)
+         !> [optional] Can A data be overwritten and destroyed?
+         logical(lk), optional, intent(in) :: overwrite_a                 
+         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space         
+         ${rt}$, optional, intent(inout), target :: storage(:)
+         !> [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_${ri}$_schur
+
+      ! Schur decomposition subroutine: real eigenvalue interface
+      module subroutine stdlib_linalg_real_eig_${ri}$_schur(a,t,z,eigvals,overwrite_a,storage,err)
+         !> Input matrix a[m,m]
+         ${rt}$, intent(inout), target :: a(:,:)
+         !> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
+         ${rt}$, intent(out), contiguous, target :: t(:,:)
+         !> Unitary/orthonormal transformation matrix Z
+         ${rt}$, optional, intent(out), contiguous, target :: z(:,:)
+         !> Output real eigenvalues that appear on the diagonal of T
+         real(${rk}$), intent(out), contiguous, target :: eigvals(:)
+         !> [optional] Provide pre-allocated workspace, size to be checked with schur_space
+         ${rt}$, optional, intent(inout), target :: storage(:)
+         !> [optional] Can A data be overwritten and destroyed?
+         logical(lk), optional, intent(in) :: overwrite_a        
+         !> [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_real_eig_${ri}$_schur 
+    #:endfor
+  end interface schur
+
+  ! Return the working array space required by the Schur decomposition solver
+  interface schur_space
+    !! version: experimental
+    !!
+    !! Computes the working array space required by the Schur decomposition solver
+    !! ([Specification](../page/specs/stdlib_linalg.html#schur-space-compute-internal-working-space-requirements-for-the-schur-decomposition))
+    !!
+    !!### Description
+    !! 
+    !! This interface returns the size of the `real` or `complex` working storage required by the 
+    !! Schur decomposition solver. The working size only depends on the kind (`real` or `complex`) and size of
+    !! the matrix being decomposed. Storage size can be used to pre-allocate a working array in case several 
+    !! repeated Schur decompositions of same-size matrices are sought. If pre-allocated working arrays 
+    !! are provided, no internal allocations will take place during the decomposition.
+    !!     
+    #:for rk,rt,ri in RC_KINDS_TYPES
+      module subroutine get_schur_${ri}$_workspace(a,lwork,err)
+         !> Input matrix a[m,m]
+         ${rt}$, intent(in), target :: a(:,:)
+         !> Minimum workspace size for the decomposition operation
+         integer(ilp), intent(out) :: lwork
+         !> State return flag. Returns an error if the query failed
+         type(linalg_state_type), optional, intent(out) :: err
+      end subroutine get_schur_${ri}$_workspace
+    #:endfor
+  end interface schur_space
 
   interface det
     !! version: experimental