#: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<n .or. nvs/=n) then
err0 = linalg_state_type(this, LINALG_VALUE_ERROR, 'Schur vectors size=',[ldvs,nvs], &
'should be n=',n)
call linalg_error_handling(err0, err)
return
end if
else
vs => vs_dummy
ldvs = size(vs, 1, kind=ilp)
nvs = size(vs, 2, kind=ilp)
end if
!> User or self-allocated storage
if (present(storage)) then
work => storage
lwork = size(work, 1, kind=ilp)
else
! Query optimal workspace
call get_schur_${ri}$_workspace(a,lwork,err0)
if (err0%error()) then
call linalg_error_handling(err0, err)
return
else
allocate(work(lwork))
end if
end if
!> SORTING: no sorting options are currently supported
sort = gees_sort_eigs(.false.)
sdim = 0_ilp
if (sort/=GEES_NOT) then
allocate(bwork(n),source=.false.)
else
bwork => bwork_dummy
end if
!> User or self-allocated eigenvalue storage
if (present(eigvals)) then
lde = size(eigvals, 1, kind=ilp)
#:if rt.startswith('c')
eigs => 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) then
err0 = linalg_state_type(this, LINALG_VALUE_ERROR, &
'Insufficient eigenvalue array size=',lde, &
'should be >=',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
