Skip to content

Commit 0222737

Browse files
perazzperazz
committed
implement Schur
1 parent 3369013 commit 0222737

File tree

2 files changed

+214
-0
lines changed

2 files changed

+214
-0
lines changed

‎src/stdlib_linalg.fypp

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -598,6 +598,80 @@ module stdlib_linalg
598598
#:endfor
599599
end interface qr_space
600600

601+
! Schur decomposition of rank-2 array A
602+
interface schur
603+
!! version: experimental
604+
!!
605+
!! Computes the Schur decomposition of matrix \( A = Z T Z^H \).
606+
!! ([Specification](../page/specs/stdlib_linalg.html#schur-compute-the-schur-decomposition-of-a-matrix))
607+
!!
608+
!!### Summary
609+
!! Compute the Schur decomposition of a `real` or `complex` matrix: \( A = Z T Z^H \), where \( Z \) is
610+
!! orthonormal/unitary and \( T \) is upper-triangular or quasi-upper-triangular. Matrix \( A \) has size `[m,m]`.
611+
!!
612+
!!### Description
613+
!!
614+
!! This interface provides methods for computing the Schur decomposition of a matrix.
615+
!! Supported data types include `real` and `complex`. If a pre-allocated workspace is provided, no internal
616+
!! memory allocations take place when using this interface.
617+
!!
618+
!! The output matrix \( T \) is upper-triangular for `complex` input matrices and quasi-upper-triangular
619+
!! for `real` input matrices (with possible \( 2 \times 2 \) blocks on the diagonal).
620+
!! The user can optionally request sorting of eigenvalues based on conditions, or a custom sorting function.
621+
!!
622+
!!@note The solution is based on LAPACK's Schur decomposition routines (`*GEES`). Sorting options
623+
!! are implemented using LAPACK's eigenvalue sorting mechanism.
624+
!!
625+
#:for rk,rt,ri in RC_KINDS_TYPES
626+
pure module subroutine stdlib_linalg_${ri}$_schur(a,t,z,sdim,output,overwrite_a,sort,storage,err)
627+
!> Input matrix a[m,m]
628+
${rt}$, intent(inout), target :: a(:,:)
629+
!> Schur form of A: upper-triangular or quasi-upper-triangular matrix T
630+
${rt}$, intent(out), contiguous, target :: t(:,:)
631+
!> Unitary/orthonormal transformation matrix Z
632+
${rt}$, intent(out), contiguous, target :: z(:,:)
633+
!> [optional] Number of eigenvalues satisfying the sort condition
634+
integer(ilp), optional, intent(out) :: sdim
635+
!> [optional] Output type: 'real' or 'complex'
636+
character(*), optional, intent(in) :: output
637+
!> [optional] Can A data be overwritten and destroyed?
638+
logical(lk), optional, intent(in) :: overwrite_a
639+
!> [optional] Sorting criterion: callable or predefined values ('lhp', 'rhp', 'iuc', 'ouc', or None)
640+
class(*), optional, intent(in) :: sort
641+
!> [optional] Provide pre-allocated workspace, size to be checked with schur_space
642+
${rt}$, intent(out), optional, target :: storage(:)
643+
!> [optional] State return flag. On error if not requested, the code will stop
644+
type(linalg_state_type), optional, intent(out) :: err
645+
end subroutine stdlib_linalg_${ri}$_schur
646+
#:endfor
647+
end interface schur
648+
649+
! Return the working array space required by the Schur decomposition solver
650+
interface schur_space
651+
!! version: experimental
652+
!!
653+
!! Computes the working array space required by the Schur decomposition solver
654+
!! ([Specification](../page/specs/stdlib_linalg.html#schur-space-compute-internal-working-space-requirements-for-the-schur-decomposition))
655+
!!
656+
!!### Description
657+
!!
658+
!! This interface returns the size of the `real` or `complex` working storage required by the
659+
!! Schur decomposition solver. The working size only depends on the kind (`real` or `complex`) and size of
660+
!! the matrix being decomposed. Storage size can be used to pre-allocate a working array in case several
661+
!! repeated Schur decompositions of same-size matrices are sought. If pre-allocated working arrays
662+
!! are provided, no internal allocations will take place during the decomposition.
663+
!!
664+
#:for rk,rt,ri in RC_KINDS_TYPES
665+
pure module subroutine get_schur_${ri}$_workspace(a,lwork,err)
666+
!> Input matrix a[m,m]
667+
${rt}$, intent(in), target :: a(:,:)
668+
!> Minimum workspace size for the decomposition operation
669+
integer(ilp), intent(out) :: lwork
670+
!> State return flag. Returns an error if the query failed
671+
type(linalg_state_type), optional, intent(out) :: err
672+
end subroutine get_schur_${ri}$_workspace
673+
#:endfor
674+
end interface schur_space
601675

602676
interface det
603677
!! version: experimental

‎src/stdlib_linalg_schur.fypp

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
1+
#:include "common.fypp"
2+
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3+
submodule (stdlib_linalg) stdlib_linalg_schur
4+
use stdlib_linalg_constants
5+
use stdlib_linalg_lapack, only: gees
6+
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
7+
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR
8+
implicit none(type,external)
9+
10+
character(*), parameter :: this = 'schur'
11+
12+
contains
13+
14+
elemental subroutine handle_gees_info(info, m, sort, err)
15+
integer(ilp), intent(in) :: info, m
16+
logical, intent(in) :: sort
17+
type(linalg_state_type), intent(out) :: err
18+
19+
! Process GEES output
20+
select case (info)
21+
case (0)
22+
! Success
23+
case (-1)
24+
err = linalg_state_type(this, LINALG_VALUE_ERROR, 'invalid matrix size m=', m)
25+
case default
26+
if (sort .and. info > 0) then
27+
err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'sorting eigenvalues failed at index ', info)
28+
else
29+
err = linalg_state_type(this, LINALG_INTERNAL_ERROR, 'GEES catastrophic error: info=', info)
30+
end if
31+
end select
32+
end subroutine handle_gees_info
33+
34+
#:for rk, rt, ri in RC_KINDS_TYPES
35+
36+
! Schur decomposition subroutine
37+
pure module subroutine stdlib_linalg_${ri}$_schur(a, t, z, lwork, overwrite_a, sort, err)
38+
!> Input matrix a[m, m]
39+
${rt}$, intent(inout), target :: a(:,:)
40+
!> Schur form of the matrix A
41+
${rt}$, intent(out), contiguous, target :: t(:,:)
42+
!> Unitary/orthogonal matrix Z
43+
${rt}$, intent(out), contiguous, target :: z(:,:)
44+
!> Workspace size (optional)
45+
integer(ilp), optional, intent(inout) :: lwork
46+
!> Overwrite input matrix A? (optional)
47+
logical(lk), optional, intent(in) :: overwrite_a
48+
!> Sorting eigenvalues? (optional)
49+
logical(lk), optional, intent(in) :: sort
50+
!> State return flag. On error if not requested, the code will stop
51+
type(linalg_state_type), optional, intent(out) :: err
52+
53+
! Local variables
54+
type(linalg_state_type) :: err0
55+
integer(ilp) :: m, lda, info, liwork
56+
logical(lk) :: overwrite_a_
57+
logical, pointer :: bwork(:)
58+
integer(ilp), allocatable :: iwork(:)
59+
${rt}$, pointer :: amat(:,:), tau(:), work(:)
60+
${rt}$ :: rwork_dummy(1) ! Dummy for real/complex cases
61+
${rt}$, allocatable :: tmat(:,:), zmat(:,:)
62+
character :: jobz
63+
64+
! Problem size
65+
m = size(a, 1, kind=ilp)
66+
lda = size(a, 1, kind=ilp)
67+
68+
! Validate dimensions
69+
if (size(a, 1, kind=ilp) /= size(a, 2, kind=ilp)) then
70+
err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Matrix A must be square: a=', [size(a,1), size(a,2)])
71+
call linalg_error_handling(err0, err)
72+
return
73+
end if
74+
75+
! Set default values
76+
overwrite_a_ = .false._lk
77+
if (present(overwrite_a)) overwrite_a_ = overwrite_a
78+
79+
! Job type based on sorting request
80+
if (present(sort) .and. sort) then
81+
jobz = 'S' ! Compute and sort eigenvalues
82+
else
83+
jobz = 'N' ! Compute Schur decomposition without sorting
84+
end if
85+
86+
! Prepare storage
87+
allocate(tmat(m, m), source=0.0_${rk}$)
88+
allocate(zmat(m, m), source=0.0_${rk}$)
89+
90+
if (overwrite_a_) then
91+
amat => a
92+
else
93+
allocate(amat(m, m), source=a)
94+
end if
95+
96+
! Allocate workspace
97+
liwork = -1
98+
if (present(lwork)) then
99+
allocate(work(lwork))
100+
else
101+
allocate(work(1)) ! Temporary workspace for querying size
102+
end if
103+
104+
! Workspace query
105+
call #{if rt.startswith('complex')}# zgees #{else}# gees #{endif}# &
106+
(jobz, 'N', nullify(bwork), m, amat, lda, tau, zmat, lda, work, liwork, iwork, info)
107+
call handle_gees_info(info, m, .false._lk, err0)
108+
if (err0%error()) then
109+
call linalg_error_handling(err0, err)
110+
return
111+
end if
112+
113+
! Update workspace size
114+
if (.not.present(lwork)) then
115+
liwork = ceiling(real(work(1), kind=${rk}$), kind=ilp)
116+
deallocate(work)
117+
allocate(work(liwork))
118+
end if
119+
120+
! Compute Schur decomposition
121+
call #{if rt.startswith('complex')}# zgees #{else}# gees #{endif}# &
122+
(jobz, 'N', nullify(bwork), m, amat, lda, tau, zmat, lda, work, liwork, iwork, info)
123+
call handle_gees_info(info, m, present(sort) .and. sort, err0)
124+
if (err0%error()) then
125+
call linalg_error_handling(err0, err)
126+
return
127+
end if
128+
129+
! Output results
130+
t = amat
131+
z = zmat
132+
133+
if (.not.overwrite_a_) deallocate(amat)
134+
if (.not.present(lwork)) deallocate(work)
135+
end subroutine stdlib_linalg_${ri}$_schur
136+
137+
#:endfor
138+
139+
end submodule stdlib_linalg_schur
140+

0 commit comments

Comments
 (0)