diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 77298e24c..ead216f2d 100644 --- a/doc/specs/stdlib_linalg.md +++ b/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 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 0bdd73c98..d017bd241 100644 --- a/example/linalg/CMakeLists.txt +++ b/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) diff --git a/example/linalg/example_schur_complex.f90 b/example/linalg/example_schur_complex.f90 new file mode 100644 index 000000000..dc4f4e913 --- /dev/null +++ b/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 + diff --git a/example/linalg/example_schur_eigvals.f90 b/example/linalg/example_schur_eigvals.f90 new file mode 100644 index 000000000..c62887439 --- /dev/null +++ b/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 + diff --git a/example/linalg/example_schur_real.f90 b/example/linalg/example_schur_real.f90 new file mode 100644 index 000000000..c46c730db --- /dev/null +++ b/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 + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f5c9c96c7..75989dbc9 100644 --- a/src/CMakeLists.txt +++ b/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 diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index e2116fb8d..2f4419a22 100644 --- a/src/stdlib_linalg.fypp +++ b/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 diff --git a/src/stdlib_linalg_schur.fypp b/src/stdlib_linalg_schur.fypp new file mode 100644 index 000000000..3e477a319 --- /dev/null +++ b/src/stdlib_linalg_schur.fypp @@ -0,0 +1,355 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +submodule (stdlib_linalg) stdlib_linalg_schur + use stdlib_linalg_constants + use stdlib_linalg_lapack, only: gees + use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, & + LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR + implicit none(type,external) + + character(*), parameter :: this = 'schur' + + !> List of internal GEES tasks: + !> No task request + character, parameter :: GEES_NOT = 'N' + + !> Request Schur vectors to be computed + character, parameter :: GEES_WITH_VECTORS = 'V' + + !> Request Schur vectors to be sorted + character, parameter :: GEES_SORTED_VECTORS = 'S' + + contains + + !> Wrapper function for Schur vectors request + elemental character function gees_vectors(wanted) + !> Are Schur vectors wanted? + logical(lk), intent(in) :: wanted + gees_vectors = merge(GEES_WITH_VECTORS,GEES_NOT,wanted) + end function gees_vectors + + !> Wrapper function for Schur vectors request + elemental character function gees_sort_eigs(sorted) + !> Should the eigenvalues be sorted? + logical(lk), intent(in) :: sorted + gees_sort_eigs = merge(GEES_SORTED_VECTORS,GEES_NOT,sorted) + end function gees_sort_eigs + + !> Wrapper function to handle GEES error codes + elemental subroutine handle_gees_info(info, m, n, ldvs, err) + integer(ilp), intent(in) :: info, m, n, ldvs + type(linalg_state_type), intent(out) :: err + + ! Process GEES output + select case (info) + case (0_ilp) + ! Success + case (-1_ilp) + ! Vector not wanted, but task is wrong + err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid Schur vector task request') + case (-2_ilp) + ! Vector not wanted, but task is wrong + err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Invalid sorting task request') + case (-4_ilp,-6_ilp) + ! Vector not wanted, but task is wrong + err = linalg_state_type(this, LINALG_VALUE_ERROR,'Invalid/non-square input matrix size:',[m,n]) + case (-11_ilp) + err = linalg_state_type(this, LINALG_VALUE_ERROR,'Schur vector matrix has insufficient size',[ldvs,n]) + case (-13_ilp) + err = linalg_state_type(this, LINALG_INTERNAL_ERROR,'Insufficient working storage size') + case (1_ilp:) + + if (info==n+2) then + err = linalg_state_type(this, LINALG_ERROR, 'Ill-conditioned problem: could not sort eigenvalues') + elseif (info==n+1) then + err = linalg_state_type(this, LINALG_ERROR, 'Some selected eigenvalues lost property due to sorting') + elseif (info==n) then + err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure: no converged eigenvalues') + else + err = linalg_state_type(this, LINALG_ERROR, 'Convergence failure; converged range is',[info,n]) + end if + + case default + + err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info) + + end select + + end subroutine handle_gees_info + + #:for rk, rt, ri in RC_KINDS_TYPES + !> Workspace query + 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 + + integer(ilp) :: m,n,sdim,info + character :: jobvs,sort + logical(lk) :: bwork_dummy(1) + ${rt}$, pointer :: amat(:,:) + real(${rk}$) :: rwork_dummy(1) + ${rt}$ :: wr_dummy(1),wi_dummy(1),vs_dummy(1,1),work_dummy(1) + type(linalg_state_type) :: err0 + + !> Initialize problem + lwork = -1_ilp + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + + !> Create a dummy intent(inout) argument + amat => a + + !> Select dummy task + jobvs = gees_vectors(.true.) + sort = gees_sort_eigs(.false.) + sdim = 0_ilp + + !> Get Schur workspace + call gees(jobvs,sort,do_not_select,n,amat,m,sdim,wr_dummy,#{if rt.startswith('r')}#wi_dummy, #{endif}#& + vs_dummy,m,work_dummy,lwork,#{if rt.startswith('c')}#rwork_dummy,#{endif}#bwork_dummy,info) + if (info==0) lwork = nint(real(work_dummy(1),kind=${rk}$),kind=ilp) + call handle_gees_info(info,m,n,m,err0) + call linalg_error_handling(err0,err) + + contains + + ! Interface to dummy select routine + pure logical(lk) function do_not_select(alpha#{if rt.startswith('r')}#r,alphai#{endif}#) + ${rt}$, intent(in) :: alpha#{if rt.startswith('r')}#r,alphai#{endif}# + do_not_select = .false. + end function do_not_select + + end subroutine get_schur_${ri}$_workspace + + ! Schur decomposition subroutine + 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] 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 + + ! Local variables + integer(ilp) :: m,n,mt,nt,ldvs,nvs,lde,lwork,sdim,info + logical(lk) :: overwrite_a_ + logical(lk), target :: bwork_dummy(1),local_eigs + logical(lk), pointer :: bwork(:) + real(${rk}$), allocatable :: rwork(:) + ${rt}$, target :: vs_dummy(1,1) + ${rt}$, pointer :: vs(:,:),work(:),eigs(:)#{if rt.startswith('r')}#,eigi(:)#{endif}# + character :: jobvs,sort + type(linalg_state_type) :: err0 + + ! Problem size + m = size(a, 1, kind=ilp) + n = size(a, 2, kind=ilp) + mt = size(t, 1, kind=ilp) + nt = size(t, 2, kind=ilp) + + ! Validate dimensions + if (m/=n .or. m<=0 .or. n<=0) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix A must be square: size(a)=',[m,n]) + call linalg_error_handling(err0, err) + return + end if + if (mt/=nt .or. mt/=n .or. nt/=n) then + err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix T must be square: size(T)=',[mt,nt], & + 'should be',[m,n]) + call linalg_error_handling(err0, err) + return + end if + + !> Copy data into the output array + t = a + + ! Can A be overwritten? By default, do not overwrite + overwrite_a_ = .false._lk + if (present(overwrite_a)) overwrite_a_ = overwrite_a .and. n>=2 + + !> Schur vectors + jobvs = gees_vectors(present(z)) + if (present(z)) then + vs => z + + ldvs = size(vs, 1, kind=ilp) + nvs = size(vs, 2, kind=ilp) + + if (ldvs eigvals + local_eigs = .false. + #:else + local_eigs = .true. + #:endif + else + local_eigs = .true. + lde = n + end if + + if (local_eigs) then + ! Use A storage if possible + if (overwrite_a_) then + eigs => a(:,1) + #:if rt.startswith('r') + eigi => a(:,2) + #:endif + else + allocate(eigs(n)#{if rt.startswith('r')}#,eigi(n)#{endif}#) + end if + endif + + #:if rt.startswith('c') + allocate(rwork(n)) + #:endif + + if (lde=',n) + + else + + ! Compute Schur decomposition + call gees(jobvs,sort,eig_select,nt,t,mt,sdim,eigs,#{if rt.startswith('r')}#eigi,#{endif}# & + vs,ldvs,work,lwork,#{if rt.startswith('c')}#rwork,#{endif}#bwork,info) + call handle_gees_info(info,m,n,m,err0) + + + end if + + eigenvalue_output: if (local_eigs) then + #:if rt.startswith('r') + ! Build complex eigenvalues + if (present(eigvals)) eigvals = cmplx(eigs,eigi,kind=${rk}$) + #:endif + if (.not.overwrite_a_) deallocate(eigs#{if rt.startswith('r')}#,eigi#{endif}#) + endif eigenvalue_output + if (.not.present(storage)) deallocate(work) + if (sort/=GEES_NOT) deallocate(bwork) + call linalg_error_handling(err0,err) + + contains + + ! Dummy select routine: currently, no sorting options are offered + pure logical(lk) function eig_select(alpha#{if rt.startswith('r')}#r,alphai#{endif}#) + #:if rt.startswith('r') + ${rt}$, intent(in) :: alphar,alphai + complex(${rk}$) :: alpha + alpha = cmplx(alphar,alphai,kind=${rk}$) + #:else + ${rt}$, intent(in) :: alpha + #:endif + eig_select = .false. + end function eig_select + + 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 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 + + type(linalg_state_type) :: err0 + integer(ilp) :: n + complex(${rk}$), allocatable :: ceigvals(:) + real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$) + real(${rk}$), parameter :: atol = tiny(0.0_${rk}$) + + n = size(eigvals,dim=1,kind=ilp) + allocate(ceigvals(n)) + + !> Compute Schur decomposition with complex eigenvalues + call stdlib_linalg_${ri}$_schur(a,t,z,ceigvals,overwrite_a,storage,err0) + + ! Check that no eigenvalues have meaningful imaginary part + if (err0%ok() .and. any(aimag(ceigvals)>atol+rtol*abs(abs(ceigvals)))) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR, & + 'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(ceigvals))) + endif + + ! Return real components only + eigvals(:n) = real(ceigvals,kind=${rk}$) + + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_real_eig_${ri}$_schur + + #:endfor + +end submodule stdlib_linalg_schur + diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 49f493924..c496bd2b3 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -12,6 +12,7 @@ set( "test_linalg_mnorm.fypp" "test_linalg_determinant.fypp" "test_linalg_qr.fypp" + "test_linalg_schur.fypp" "test_linalg_svd.fypp" "test_linalg_matrix_property_checks.fypp" "test_linalg_sparse.fypp" @@ -30,6 +31,7 @@ ADDTEST(linalg_mnorm) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(linalg_qr) +ADDTEST(linalg_schur) ADDTEST(linalg_svd) ADDTEST(blas_lapack) ADDTEST(linalg_sparse) diff --git a/test/linalg/test_linalg_schur.fypp b/test/linalg/test_linalg_schur.fypp new file mode 100644 index 000000000..139fdf360 --- /dev/null +++ b/test/linalg/test_linalg_schur.fypp @@ -0,0 +1,238 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +! Test Schur decomposition +module test_linalg_schur + use testdrive, only: error_type, check, new_unittest, unittest_type + use stdlib_linalg_constants + use stdlib_linalg_state, only: LINALG_VALUE_ERROR,linalg_state_type + use stdlib_linalg, only: schur,schur_space + use ieee_arithmetic, only: ieee_value,ieee_quiet_nan + + implicit none (type,external) + + public :: test_schur_decomposition + + contains + + !> schur decomposition tests + subroutine test_schur_decomposition(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in RC_KINDS_TYPES + tests = [tests,new_unittest("schur_api_${ri}$",test_schur_api_${ri}$), & + new_unittest("schur_random_${ri}$",test_schur_random_${ri}$), & + new_unittest("schur_symmetric_${ri}$",test_schur_symmetric_${ri}$)] + #:endfor + + end subroutine test_schur_decomposition + + !> schur decomposition of a random matrix + #:for rk,rt,ri in RC_KINDS_TYPES + subroutine test_schur_api_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 15_ilp + integer(ilp) :: lwork + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) + complex(${rk}$) :: eigs(n) + ${rt}$, dimension(n,n) :: a,t,z + ${rt}$, allocatable :: storage(:) + #:if 'complex' in rt + real(${rk}$) :: rea(n,n),ima(n,n) + #:endif + type(linalg_state_type) :: state + + #:if 'complex' in rt + call random_number(rea) + call random_number(ima) + a = cmplx(rea,ima,kind=${rk}$) + #:else + call random_number(a) + #:endif + + ! Test simple API + call schur(a,t,err=state) + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Test output transformation matrix + call schur(a,t,z,err=state) + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Test output eigenvalues + call schur(a,t,eigvals=eigs,err=state) + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Test storage query + call schur_space(a,lwork,err=state) + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Test with user-defined storage + allocate(storage(lwork)) + call schur(a,t,eigvals=eigs,storage=storage,err=state) + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + end subroutine test_schur_api_${ri}$ + + subroutine test_schur_random_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + real(${rk}$), parameter :: rtol = 1.0e-4_${rk}$ + real(${rk}$), parameter :: eps = sqrt(epsilon(0.0_${rk}$)) + integer(ilp) :: lwork + ${rt}$, allocatable :: storage(:) + ${rt}$, dimension(n,n) :: a,t,z,aorig + #:if 'complex' in rt + real(${rk}$), dimension(n,n) :: a_re,a_im + #:endif + type(linalg_state_type) :: state + + #:if 'complex' in rt + call random_number(a_re) + call random_number(a_im) + a = cmplx(a_re,a_im,kind=${rk}$) + #:else + call random_number(a) + #:endif + aorig = a + + ! 1) Run schur (standard) + call schur(a,t,z,err=state) + + ! Check return code + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(schur_error(a,z,t)<=max(rtol*abs(a),eps)), & + 'converged solution (${rt}$)') + if (allocated(error)) return + + ! 2) Run schur (overwrite A) + call schur(a,t,z,overwrite_a=.true.,err=state) + + ! Check return code + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(schur_error(aorig,z,t)<=max(rtol*abs(aorig),eps)), & + 'converged solution (${rt}$ - overwrite A)') + if (allocated(error)) return + + ! 3) Use working storage + a = aorig + call schur_space(a,lwork,err=state) + + ! Check return code + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + allocate(storage(lwork)) + call schur(a,t,z,storage=storage,err=state) + + ! Check return code + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(schur_error(a,z,t)<=max(rtol*abs(a),eps)), & + 'converged solution (${rt}$ - external storage)') + if (allocated(error)) return + + contains + + pure function schur_error(a,z,t) result(err) + ${rt}$, intent(in), dimension(:,:) :: a,z,t + real(${rk}$), dimension(size(a,1),size(a,2)) :: err + + #:if rt.startswith('real') + err = abs(matmul(matmul(z,t),transpose(z)) - a) + #:else + err = abs(matmul(matmul(z,t),conjg(transpose(z))) - a) + #:endif + end function schur_error + + end subroutine test_schur_random_${ri}$ + + !> Test symmetric matrix (real eigenvalues) + subroutine test_schur_symmetric_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + integer(ilp), parameter :: n = 3_ilp + real(${rk}$), parameter :: rtol = 1.0e-4_${rk}$ + real(${rk}$), parameter :: eps = sqrt(epsilon(0.0_${rk}$)) + real(${rk}$) :: reigs(n) + ${rt}$, dimension(n,n) :: a, t, z + type(linalg_state_type) :: state + + ! Define a symmetric 3x3 matrix with real eigenvalues + a = reshape([ 3, 1, 0, & + 1, 3, 1, & + 0, 1, 3], shape=[n, n]) + + ! Return real eigenvalues (Should trigger an error if they have an imaginary part) + call schur(a, t, z, eigvals=reigs, err=state) + + ! Check return code + call check(error, state%ok(), state%print()) + if (allocated(error)) return + + ! Check solution + call check(error, all(schur_error(a, z, t) <= max(rtol * abs(a), eps)), & + 'converged solution (real symmetric, real eigs)') + if (allocated(error)) return + + contains + + pure function schur_error(a,z,t) result(err) + ${rt}$, intent(in), dimension(:,:) :: a,z,t + real(${rk}$), dimension(size(a,1),size(a,2)) :: err + + #:if rt.startswith('real') + err = abs(matmul(matmul(z,t),transpose(z)) - a) + #:else + err = abs(matmul(matmul(z,t),conjg(transpose(z))) - a) + #:endif + end function schur_error + + end subroutine test_schur_symmetric_${ri}$ + + + #:endfor + +end module test_linalg_schur + +program test_schur + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_schur, only : test_schur_decomposition + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_schur", test_schur_decomposition) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_schur