#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
#:set EIG_PROBLEM = ["standard", "generalized"]
#:set EIG_FUNCTION = ["geev","ggev"]
#:set EIG_PROBLEM_LIST = list(zip(EIG_PROBLEM, EIG_FUNCTION))
submodule (stdlib_linalg) stdlib_linalg_eigenvalues
!! Compute eigenvalues and eigenvectors
use stdlib_linalg_constants
use stdlib_linalg_lapack, only: geev, ggev, heev, syev
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_positive_inf, ieee_quiet_nan
implicit none
character(*), parameter :: this = 'eigenvalues'
!> Utility function: Scale generalized eigenvalue
interface scale_general_eig
#:for rk,rt,ri in RC_KINDS_TYPES
module procedure scale_general_eig_${ri}$
#:endfor
end interface scale_general_eig
contains
!> Request for eigenvector calculation
elemental character function eigenvectors_task(required)
logical(lk), intent(in) :: required
eigenvectors_task = merge('V','N',required)
end function eigenvectors_task
!> Request for symmetry side (default: lower)
elemental character function symmetric_triangle_task(upper)
logical(lk), optional, intent(in) :: upper
symmetric_triangle_task = 'L'
if (present(upper)) symmetric_triangle_task = merge('U','L',upper)
end function symmetric_triangle_task
!> Process GEEV output flags
pure subroutine handle_geev_info(err,info,shapea)
!> Error handler
type(linalg_state_type), intent(inout) :: err
!> GEEV return flag
integer(ilp), intent(in) :: info
!> Input matrix size
integer(ilp), intent(in) :: shapea(2)
select case (info)
case (0)
! Success!
err%state = LINALG_SUCCESS
case (-1)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.')
case (-2)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.')
case (-5,-3)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea)
case (-9)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.')
case (-11)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.')
case (-13)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.')
case (1:)
err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
case default
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.')
end select
end subroutine handle_geev_info
!> Process GGEV output flags
pure subroutine handle_ggev_info(err,info,shapea,shapeb)
!> Error handler
type(linalg_state_type), intent(inout) :: err
!> GEEV return flag
integer(ilp), intent(in) :: info
!> Input matrix size
integer(ilp), intent(in) :: shapea(2),shapeb(2)
select case (info)
case (0)
! Success!
err%state = LINALG_SUCCESS
case (-1)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.')
case (-2)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.')
case (-5,-3)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',shapea)
case (-7)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: b=',shapeb)
case (-12)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.')
case (-14)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.')
case (-16)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.')
case (1:)
err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
case default
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by ggev.')
end select
end subroutine handle_ggev_info
!> Process SYEV/HEEV output flags
elemental subroutine handle_heev_info(err,info,m,n)
!> Error handler
type(linalg_state_type), intent(inout) :: err
!> SYEV/HEEV return flag
integer(ilp), intent(in) :: info
!> Input matrix size
integer(ilp), intent(in) :: m,n
select case (info)
case (0)
! Success!
err%state = LINALG_SUCCESS
case (-1)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.')
case (-2)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.')
case (-5,-3)
err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n])
case (-8)
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.')
case (1:)
err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.')
case default
err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.')
end select
end subroutine handle_heev_info
#:for rk,rt,ri in RC_KINDS_TYPES
#:for ep,ei in EIG_PROBLEM_LIST
module function stdlib_linalg_eigvals_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,err) result(lambda)
!! Return an array of eigenvalues of matrix A.
!> Input matrix A[m,n]
${rt}$, intent(in), dimension(:,:), target :: a
#:if ei=='ggev'
!> Generalized problem matrix B[n,n]
${rt}$, intent(inout), dimension(:,:), target :: b
#:endif
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out) :: err
!> Array of eigenvalues
complex(${rk}$), allocatable :: lambda(:)
!> Create
${rt}$, pointer, dimension(:,:) :: amat#{if ei=='ggev'}#, bmat #{endif}#
integer(ilp) :: m,n,k
!> Create an internal pointer so the intent of A won't affect the next call
amat => a
#{if ei=='ggev'}#bmat => b#{endif}#
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)
!> Allocate return storage
allocate(lambda(k))
!> Compute eigenvalues only
call stdlib_linalg_eig_${ep}$_${ri}$(amat#{if ei=='ggev'}#,bmat#{endif}#,lambda,err=err)
end function stdlib_linalg_eigvals_${ep}$_${ri}$
module function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#) result(lambda)
!! Return an array of eigenvalues of matrix A.
!> Input matrix A[m,n]
${rt}$, intent(in), dimension(:,:), target :: a
#:if ei=='ggev'
!> Generalized problem matrix B[n,n]
${rt}$, intent(inout), dimension(:,:), target :: b
#:endif
!> Array of eigenvalues
complex(${rk}$), allocatable :: lambda(:)
!> Create
${rt}$, pointer, dimension(:,:) :: amat#{if ei=='ggev'}#, bmat #{endif}#
integer(ilp) :: m,n,k
!> Create an internal pointer so the intent of A won't affect the next call
amat => a
#{if ei=='ggev'}#bmat => b#{endif}#
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)
!> Allocate return storage
allocate(lambda(k))
!> Compute eigenvalues only
call stdlib_linalg_eig_${ep}$_${ri}$(amat#{if ei=='ggev'}#,bmat#{endif}#,lambda,overwrite_a=.false.)
end function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$
module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left,&
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
!! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues,
!! and optionally right or left eigenvectors.
!> Input matrix A[m,n]
${rt}$, intent(inout), dimension(:,:), target :: a
#:if ei=='ggev'
!> Generalized problem matrix B[n,n]
${rt}$, intent(inout), dimension(:,:), target :: b
#:endif
!> Array of eigenvalues
complex(${rk}$), intent(out) :: lambda(:)
!> [optional] RIGHT eigenvectors of A (as columns)
complex(${rk}$), optional, intent(out), target :: right(:,:)
!> [optional] LEFT eigenvectors of A (as columns)
complex(${rk}$), optional, intent(out), target :: left(:,:)
!> [optional] Can A data be overwritten and destroyed? (default: no)
logical(lk), optional, intent(in) :: overwrite_a
#:if ei=='ggev'
!> [optional] Can B data be overwritten and destroyed? (default: no)
logical(lk), optional, intent(in) :: overwrite_b
#:endif
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
!> Local variables
type(linalg_state_type) :: err0
integer(ilp) :: m,n,lda,ldu,ldv,info,k,lwork,neig#{if ei=='ggev'}#,ldb,nb#{endif}#
logical(lk) :: copy_a#{if ei=='ggev'}#,copy_b#{endif}#
character :: task_u,task_v
${rt}$, target :: work_dummy(1),u_dummy(1,1),v_dummy(1,1)
${rt}$, allocatable :: work(:)
${rt}$, dimension(:,:), pointer :: amat,umat,vmat#{if ei=='ggev'}#,bmat#{endif}#
#:if rt.startswith('complex')
real(${rk}$), allocatable :: rwork(:)
#:else
${rt}$, pointer :: lreal(:),limag(:)
#:endif
#:if ei=='ggev'
${rt}$, allocatable :: beta(:)
#:endif
!> Matrix size
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)
neig = size(lambda,kind=ilp)
lda = m
if (k<=0 .or. m/=n) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,&
'invalid or matrix size a=',[m,n],', must be nonempty square.')
call linalg_error_handling(err0,err)
return
elseif (neig<k) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,&
'eigenvalue array has insufficient size:',&
' lambda=',neig,', n=',n)
call linalg_error_handling(err0,err)
return
endif
#:if ep=='generalized'
ldb = size(b,1,kind=ilp)
nb = size(b,2,kind=ilp)
if (ldb/=n .or. nb/=n) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,&
'invalid or matrix size b=',[ldb,nb],', must be same as a=',[m,n])
call linalg_error_handling(err0,err)
return
end if
#:endif
! Can A be overwritten? By default, do not overwrite
copy_a = .true._lk
if (present(overwrite_a)) copy_a = .not.overwrite_a
! Initialize a matrix temporary
if (copy_a) then
allocate(amat(m,n),source=a)
else
amat => a
endif
#:if ep=='generalized'
! Can B be overwritten? By default, do not overwrite
copy_b = .true._lk
if (present(overwrite_b)) copy_b = .not.overwrite_b
! Initialize a matrix temporary
if (copy_b) then
allocate(bmat,source=b)
else
bmat => b
endif
allocate(beta(n))
#:endif
! Decide if U, V eigenvectors
task_u = eigenvectors_task(present(left))
task_v = eigenvectors_task(present(right))
if (present(right)) then
#:if rt.startswith('complex')
! For a complex matrix, GEEV returns complex arrays.
! Point directly to output.
vmat => right
#:else
! For a real matrix, GEEV returns real arrays.
! Allocate temporary reals, will be converted to complex vectors at the end.
allocate(vmat(n,n))
#:endif
if (size(vmat,2,kind=ilp)<n) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,&
'right eigenvector matrix has insufficient size: ',&
shape(vmat),', with n=',n)
endif
else
vmat => v_dummy
endif
if (present(left)) then
#:if rt.startswith('complex')
! For a complex matrix, GEEV returns complex arrays.
! Point directly to output.
umat => left
#:else
! For a real matrix, GEEV returns real arrays.
! Allocate temporary reals, will be converted to complex vectors at the end.
allocate(umat(n,n))
#:endif
if (size(umat,2,kind=ilp)<n) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,&
'left eigenvector matrix has insufficient size: ',&
shape(umat),', with n=',n)
endif
else
umat => u_dummy
endif
get_${ei}$: if (err0%ok()) then
ldu = size(umat,1,kind=ilp)
ldv = size(vmat,1,kind=ilp)
! Compute workspace size
#:if rt.startswith('complex')
allocate(rwork(2*n))
#:else
allocate(lreal(n),limag(n))
#:endif
lwork = -1_ilp
call ${ei}$(task_u,task_v,n,amat,lda,&
#:if ep=='generalized'
bmat,ldb, &
#:endif
#{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# &
#:if ep=='generalized'
beta, &
#:endif
umat,ldu,vmat,ldv,&
work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info)
call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
! Compute eigenvalues
if (info==0) then
!> Prepare working storage
lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp)
allocate(work(lwork))
!> Compute eigensystem
call ${ei}$(task_u,task_v,n,amat,lda,&
#:if ep=='generalized'
bmat,ldb, &
#:endif
#{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# &
#:if ep=='generalized'
beta, &
#:endif
umat,ldu,vmat,ldv,&
work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info)
call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
endif
! Finalize storage and process output flag
#:if not rt.startswith('complex')
lambda(:n) = cmplx(lreal(:n),limag(:n),kind=${rk}$)
#:endif
#:if ep=='generalized'
! Scale generalized eigenvalues
lambda(:n) = scale_general_eig(lambda(:n),beta)
#:endif
#:if not rt.startswith('complex')
! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair,
! ${ei}$ returns reals as:
! u(j) = VL(:,j) + i*VL(:,j+1) and
! u(j+1) = VL(:,j) - i*VL(:,j+1).
! Convert these to complex numbers here.
if (present(right)) call assign_real_eigenvectors_${rk}$(n,lambda,vmat,right)
if (present(left)) call assign_real_eigenvectors_${rk}$(n,lambda,umat,left)
#:endif
endif get_${ei}$
if (copy_a) deallocate(amat)
#:if ep=='generalized'
if (copy_b) deallocate(bmat)
#:endif
#:if not rt.startswith('complex')
if (present(right)) deallocate(vmat)
if (present(left)) deallocate(umat)
#:endif
call linalg_error_handling(err0,err)
end subroutine stdlib_linalg_eig_${ep}$_${ri}$
#:endfor
module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda)
!! Return an array of eigenvalues of real symmetric / complex hermitian A
!> Input matrix A[m,n]
${rt}$, intent(in), target :: a(:,:)
!> [optional] Should the upper/lower half of A be used? Default: lower
logical(lk), optional, intent(in) :: upper_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), intent(out) :: err
!> Array of eigenvalues
real(${rk}$), allocatable :: lambda(:)
${rt}$, pointer :: amat(:,:)
integer(ilp) :: m,n,k
!> Create an internal pointer so the intent of A won't affect the next call
amat => a
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)
!> Allocate return storage
allocate(lambda(k))
!> Compute eigenvalues only
call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,err=err)
end function stdlib_linalg_eigvalsh_${ri}$
module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda)
!! Return an array of eigenvalues of real symmetric / complex hermitian A
!> Input matrix A[m,n]
${rt}$, intent(in), target :: a(:,:)
!> [optional] Should the upper/lower half of A be used? Default: use lower
logical(lk), optional, intent(in) :: upper_a
!> Array of eigenvalues
real(${rk}$), allocatable :: lambda(:)
${rt}$, pointer :: amat(:,:)
integer(ilp) :: m,n,k
!> Create an internal pointer so the intent of A won't affect the next call
amat => a
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)
!> Allocate return storage
allocate(lambda(k))
!> Compute eigenvalues only
call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.)
end function stdlib_linalg_eigvalsh_noerr_${ri}$
module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err)
!! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda`
!! of eigenvalues, and optionally right or left eigenvectors.
!> Input matrix A[m,n]
${rt}$, intent(inout), target :: a(:,:)
!> Array of eigenvalues
real(${rk}$), intent(out) :: lambda(:)
!> The columns of vectors contain the orthonormal eigenvectors of A
${rt}$, optional, intent(out), target :: vectors(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
!> [optional] Should the upper/lower half of A be used? Default: use lower
logical(lk), optional, intent(in) :: upper_a
!> [optional] state return flag. On error if not requested, the code will stop
type(linalg_state_type), optional, intent(out) :: err
!> Local variables
type(linalg_state_type) :: err0
integer(ilp) :: m,n,lda,info,k,lwork,neig
logical(lk) :: copy_a
character :: triangle,task
${rt}$, target :: work_dummy(1)
${rt}$, allocatable :: work(:)
#:if rt.startswith('complex')
real(${rk}$), allocatable :: rwork(:)
#:endif
${rt}$, pointer :: amat(:,:)
!> Matrix size
m = size(a,1,kind=ilp)
n = size(a,2,kind=ilp)
k = min(m,n)
neig = size(lambda,kind=ilp)
if (k<=0 .or. m/=n) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], &
', must be non-empty square.')
call linalg_error_handling(err0,err)
return
elseif (neig<k) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',&
' lambda=',neig,' must be >= n=',n)
call linalg_error_handling(err0,err)
return
endif
! Check if input A can be overwritten
copy_a = .true._lk
if (present(vectors)) then
! No need to copy A anyways
copy_a = .false.
elseif (present(overwrite_a)) then
copy_a = .not.overwrite_a
endif
! Should we use the upper or lower half of the matrix?
triangle = symmetric_triangle_task(upper_a)
! Request for eigenvectors
task = eigenvectors_task(present(vectors))
if (present(vectors)) then
! Check size
if (any(shape(vectors,kind=ilp)<n)) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,&
'eigenvector matrix has insufficient size: ',&
shape(vectors),', with n=',n)
call linalg_error_handling(err0,err)
return
endif
! The input matrix will be overwritten by the vectors.
! So, use this one as storage for syev/heev
amat => vectors
! Copy data in
amat(:n,:n) = a(:n,:n)
elseif (copy_a) then
! Initialize a matrix temporary
allocate(amat(m,n),source=a)
else
! Overwrite A
amat => a
endif
lda = size(amat,1,kind=ilp)
! Request workspace size
lwork = -1_ilp
#:if rt.startswith('complex')
allocate(rwork(max(1,3*n-2)))
call heev(task,triangle,n,amat,lda,lambda,work_dummy,lwork,rwork,info)
#:else
call syev(task,triangle,n,amat,lda,lambda,work_dummy,lwork,info)
#:endif
call handle_heev_info(err0,info,m,n)
! Compute eigenvalues
if (info==0) then
!> Prepare working storage
lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp)
allocate(work(lwork))
!> Compute eigensystem
#:if rt.startswith('complex')
call heev(task,triangle,n,amat,lda,lambda,work,lwork,rwork,info)
#:else
call syev(task,triangle,n,amat,lda,lambda,work,lwork,info)
#:endif
call handle_heev_info(err0,info,m,n)
endif
! Finalize storage and process output flag
if (copy_a) deallocate(amat)
call linalg_error_handling(err0,err)
end subroutine stdlib_linalg_eigh_${ri}$
#:endfor
#:for rk,rt,ri in REAL_KINDS_TYPES
pure subroutine assign_real_eigenvectors_${rk}$(n,lambda,lmat,out_mat)
!! GEEV for real matrices returns complex eigenvalues in real arrays, where two consecutive
!! reals at [j,j+1] locations represent the real and imaginary parts of two complex conjugate
!! eigenvalues. Convert them to complex here, following the GEEV logic.
!> Problem size
integer(ilp), intent(in) :: n
!> Array of eigenvalues
complex(${rk}$), intent(in) :: lambda(:)
!> Real matrix as returned by geev
${rt}$, intent(in) :: lmat(:,:)
!> Complex matrix as returned by eig
complex(${rk}$), intent(out) :: out_mat(:,:)
integer(ilp) :: i,j
! Copy matrix
do concurrent(i=1:n,j=1:n)
out_mat(i,j) = lmat(i,j)
end do
! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair,
! geev returns them as reals as:
! u(j) = VL(:,j) + i*VL(:,j+1) and
! u(j+1) = VL(:,j) - i*VL(:,j+1).
! Convert these to complex numbers here.
do j=1,n-1
if (lambda(j)==conjg(lambda(j+1))) then
out_mat(:, j) = cmplx(lmat(:,j), lmat(:,j+1),kind=${rk}$)
out_mat(:,j+1) = cmplx(lmat(:,j),-lmat(:,j+1),kind=${rk}$)
endif
end do
end subroutine assign_real_eigenvectors_${rk}$
#:for ep,ei in EIG_PROBLEM_LIST
module subroutine stdlib_linalg_real_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
!! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues,
!! and optionally right or left eigenvectors. Returns an error if the eigenvalues had
!! non-trivial imaginary parts.
!> Input matrix A[m,n]
${rt}$, intent(inout), target :: a(:,:)
#:if ei=='ggev'
!> Generalized problem matrix B[n,n]
${rt}$, intent(inout), target :: b(:,:)
#:endif
!> Array of real eigenvalues
real(${rk}$), intent(out) :: lambda(:)
!> The columns of RIGHT contain the right eigenvectors of A
complex(${rk}$), optional, intent(out), target :: right(:,:)
!> The columns of LEFT contain the left eigenvectors of A
complex(${rk}$), optional, intent(out), target :: left(:,:)
!> [optional] Can A data be overwritten and destroyed?
logical(lk), optional, intent(in) :: overwrite_a
#:if ei=='ggev'
!> [optional] Can B data be overwritten and destroyed? (default: no)
logical(lk), optional, intent(in) :: overwrite_b
#:endif
!> [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 :: clambda(:)
real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$)
real(${rk}$), parameter :: atol = tiny(0.0_${rk}$)
n = size(lambda,dim=1,kind=ilp)
allocate(clambda(n))
call stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,clambda,right,left, &
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err0)
! Check that no eigenvalues have meaningful imaginary part
if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then
err0 = linalg_state_type(this,LINALG_VALUE_ERROR, &
'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(clambda)))
endif
! Return real components only
lambda(:n) = real(clambda,kind=${rk}$)
call linalg_error_handling(err0,err)
end subroutine stdlib_linalg_real_eig_${ep}$_${ri}$
#:endfor
#:endfor
#:for rk,rt,ri in RC_KINDS_TYPES
!> Utility function: Scale generalized eigenvalue
elemental complex(${rk}$) function scale_general_eig_${ri}$(alpha,beta) result(lambda)
!! A generalized eigenvalue for a pair of matrices (A,B) is a scalar lambda or a ratio
!! alpha/beta = lambda, such that A - lambda*B is singular. It is usually represented as the
!! pair (alpha,beta), there is a reasonable interpretation for beta=0, and even for both
!! being zero.
complex(${rk}$), intent(in) :: alpha
${rt}$, intent(in) :: beta
real (${rk}$), parameter :: rzero = 0.0_${rk}$
complex(${rk}$), parameter :: czero = (0.0_${rk}$,0.0_${rk}$)
if (beta==#{if rt.startswith('real')}#r#{else}#c#{endif}#zero) then
if (alpha/=czero) then
lambda = cmplx(ieee_value(1.0_${rk}$, ieee_positive_inf), &
ieee_value(1.0_${rk}$, ieee_positive_inf), kind=${rk}$)
else
lambda = ieee_value(1.0_${rk}$, ieee_quiet_nan)
end if
else
lambda = alpha/beta
end if
end function scale_general_eig_${ri}$
#:endfor
end submodule stdlib_linalg_eigenvalues
