#:include "common.fypp"
submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_gen
implicit none
contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
module subroutine stdlib${ii}$_sgeev( jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr,ldvr, work, lwork, &
!! SGEEV computes for an N-by-N real nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate-transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobvl, jobvr
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
! Array Arguments
real(sp), intent(inout) :: a(lda,*)
real(sp), intent(out) :: vl(ldvl,*), vr(ldvr,*), wi(*), work(*), wr(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr
character :: side
integer(${ik}$) :: hswork, i, ibal, ierr, ihi, ilo, itau, iwrk, k, lwork_trevc, maxwrk, &
minwrk, nout
real(sp) :: anrm, bignum, cs, cscale, eps, r, scl, smlnum, sn
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(sp) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -1_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -9_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -11_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! nb refers to the optimal block size for the immediately
! following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_shseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = 2_${ik}$*n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'SGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
if( wantvl ) then
minwrk = 4_${ik}$*n
maxwrk = max( maxwrk, 2_${ik}$*n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$,'SORGHR', ' ', n, 1_${ik}$, n,&
-1_${ik}$ ) )
call stdlib${ii}$_shseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vl, ldvl,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
call stdlib${ii}$_strevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr, n, nout,&
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
maxwrk = max( maxwrk, 4_${ik}$*n )
else if( wantvr ) then
minwrk = 4_${ik}$*n
maxwrk = max( maxwrk, 2_${ik}$*n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$,'SORGHR', ' ', n, 1_${ik}$, n,&
-1_${ik}$ ) )
call stdlib${ii}$_shseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
call stdlib${ii}$_strevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr, n, nout,&
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
maxwrk = max( maxwrk, 4_${ik}$*n )
else
minwrk = 3_${ik}$*n
call stdlib${ii}$_shseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
end if
maxwrk = max( maxwrk, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -13_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGEEV ', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_slamch( 'P' )
smlnum = stdlib${ii}$_slamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_slabad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
anrm = stdlib${ii}$_slange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix
! (workspace: need n)
ibal = 1_${ik}$
call stdlib${ii}$_sgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
! reduce to upper hessenberg form
! (workspace: need 3*n, prefer 2*n+n*nb)
itau = ibal + n
iwrk = itau + n
call stdlib${ii}$_sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_slacpy( 'L', n, n, a, lda, vl, ldvl )
! generate orthogonal matrix in vl
! (workspace: need 3*n-1, prefer 2*n+(n-1)*nb)
call stdlib${ii}$_sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,work( iwrk ), &
lwork-iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_slacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_slacpy( 'L', n, n, a, lda, vr, ldvr )
! generate orthogonal matrix in vr
! (workspace: need 3*n-1, prefer 2*n+(n-1)*nb)
call stdlib${ii}$_sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
else
! compute eigenvalues only
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_shseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_shseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (workspace: need 4*n, prefer n + n + 2*n*nb)
call stdlib${ii}$_strevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1, ierr )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
! (workspace: need n)
call stdlib${ii}$_sgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_snrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, &
vl( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_sscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vl( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( iwrk+k-1 ) = vl( k, i )**2_${ik}$ + vl( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_isamax( n, work( iwrk ), 1_${ik}$ )
call stdlib${ii}$_slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
call stdlib${ii}$_srot( n, vl( 1_${ik}$, i ), 1_${ik}$, vl( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vl( k, i+1 ) = zero
end if
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
! (workspace: need n)
call stdlib${ii}$_sgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_snrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, &
vr( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_sscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vr( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( iwrk+k-1 ) = vr( k, i )**2_${ik}$ + vr( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_isamax( n, work( iwrk ), 1_${ik}$ )
call stdlib${ii}$_slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
call stdlib${ii}$_srot( n, vr( 1_${ik}$, i ), 1_${ik}$, vr( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vr( k, i+1 ) = zero
end if
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wr( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wi( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
if( info>0_${ik}$ ) then
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wr, n,ierr )
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wi, n,ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_sgeev
module subroutine stdlib${ii}$_dgeev( jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr,ldvr, work, lwork, &
!! DGEEV computes for an N-by-N real nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate-transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobvl, jobvr
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
! Array Arguments
real(dp), intent(inout) :: a(lda,*)
real(dp), intent(out) :: vl(ldvl,*), vr(ldvr,*), wi(*), work(*), wr(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr
character :: side
integer(${ik}$) :: hswork, i, ibal, ierr, ihi, ilo, itau, iwrk, k, lwork_trevc, maxwrk, &
minwrk, nout
real(dp) :: anrm, bignum, cs, cscale, eps, r, scl, smlnum, sn
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(dp) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -1_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -9_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -11_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! nb refers to the optimal block size for the immediately
! following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_dhseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = 2_${ik}$*n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
if( wantvl ) then
minwrk = 4_${ik}$*n
maxwrk = max( maxwrk, 2_${ik}$*n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$,'DORGHR', ' ', n, 1_${ik}$, n,&
-1_${ik}$ ) )
call stdlib${ii}$_dhseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vl, ldvl,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
call stdlib${ii}$_dtrevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr, n, nout,&
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
maxwrk = max( maxwrk, 4_${ik}$*n )
else if( wantvr ) then
minwrk = 4_${ik}$*n
maxwrk = max( maxwrk, 2_${ik}$*n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$,'DORGHR', ' ', n, 1_${ik}$, n,&
-1_${ik}$ ) )
call stdlib${ii}$_dhseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
call stdlib${ii}$_dtrevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr, n, nout,&
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
maxwrk = max( maxwrk, 4_${ik}$*n )
else
minwrk = 3_${ik}$*n
call stdlib${ii}$_dhseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
end if
maxwrk = max( maxwrk, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -13_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGEEV ', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_dlamch( 'P' )
smlnum = stdlib${ii}$_dlamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_dlabad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
anrm = stdlib${ii}$_dlange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix
! (workspace: need n)
ibal = 1_${ik}$
call stdlib${ii}$_dgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
! reduce to upper hessenberg form
! (workspace: need 3*n, prefer 2*n+n*nb)
itau = ibal + n
iwrk = itau + n
call stdlib${ii}$_dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_dlacpy( 'L', n, n, a, lda, vl, ldvl )
! generate orthogonal matrix in vl
! (workspace: need 3*n-1, prefer 2*n+(n-1)*nb)
call stdlib${ii}$_dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,work( iwrk ), &
lwork-iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_dlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_dlacpy( 'L', n, n, a, lda, vr, ldvr )
! generate orthogonal matrix in vr
! (workspace: need 3*n-1, prefer 2*n+(n-1)*nb)
call stdlib${ii}$_dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
else
! compute eigenvalues only
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_dhseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_dhseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (workspace: need 4*n, prefer n + n + 2*n*nb)
call stdlib${ii}$_dtrevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1, ierr )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
! (workspace: need n)
call stdlib${ii}$_dgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, &
vl( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_dscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vl( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( iwrk+k-1 ) = vl( k, i )**2_${ik}$ + vl( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_idamax( n, work( iwrk ), 1_${ik}$ )
call stdlib${ii}$_dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
call stdlib${ii}$_drot( n, vl( 1_${ik}$, i ), 1_${ik}$, vl( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vl( k, i+1 ) = zero
end if
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
! (workspace: need n)
call stdlib${ii}$_dgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, &
vr( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_dscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vr( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( iwrk+k-1 ) = vr( k, i )**2_${ik}$ + vr( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_idamax( n, work( iwrk ), 1_${ik}$ )
call stdlib${ii}$_dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
call stdlib${ii}$_drot( n, vr( 1_${ik}$, i ), 1_${ik}$, vr( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vr( k, i+1 ) = zero
end if
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wr( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wi( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
if( info>0_${ik}$ ) then
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wr, n,ierr )
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wi, n,ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_dgeev
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$geev( jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr,ldvr, work, lwork, &
!! DGEEV: computes for an N-by-N real nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate-transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobvl, jobvr
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
! Array Arguments
real(${rk}$), intent(inout) :: a(lda,*)
real(${rk}$), intent(out) :: vl(ldvl,*), vr(ldvr,*), wi(*), work(*), wr(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr
character :: side
integer(${ik}$) :: hswork, i, ibal, ierr, ihi, ilo, itau, iwrk, k, lwork_trevc, maxwrk, &
minwrk, nout
real(${rk}$) :: anrm, bignum, cs, cscale, eps, r, scl, smlnum, sn
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(${rk}$) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -1_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -9_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -11_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! nb refers to the optimal block size for the immediately
! following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_${ri}$hseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = 2_${ik}$*n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
if( wantvl ) then
minwrk = 4_${ik}$*n
maxwrk = max( maxwrk, 2_${ik}$*n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$,'DORGHR', ' ', n, 1_${ik}$, n,&
-1_${ik}$ ) )
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vl, ldvl,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
call stdlib${ii}$_${ri}$trevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr, n, nout,&
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
maxwrk = max( maxwrk, 4_${ik}$*n )
else if( wantvr ) then
minwrk = 4_${ik}$*n
maxwrk = max( maxwrk, 2_${ik}$*n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$,'DORGHR', ' ', n, 1_${ik}$, n,&
-1_${ik}$ ) )
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
call stdlib${ii}$_${ri}$trevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr, n, nout,&
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
maxwrk = max( maxwrk, 4_${ik}$*n )
else
minwrk = 3_${ik}$*n
call stdlib${ii}$_${ri}$hseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + 1_${ik}$, n + hswork )
end if
maxwrk = max( maxwrk, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -13_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGEEV ', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_${ri}$lamch( 'P' )
smlnum = stdlib${ii}$_${ri}$lamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_${ri}$labad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
anrm = stdlib${ii}$_${ri}$lange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_${ri}$lascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix
! (workspace: need n)
ibal = 1_${ik}$
call stdlib${ii}$_${ri}$gebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
! reduce to upper hessenberg form
! (workspace: need 3*n, prefer 2*n+n*nb)
itau = ibal + n
iwrk = itau + n
call stdlib${ii}$_${ri}$gehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_${ri}$lacpy( 'L', n, n, a, lda, vl, ldvl )
! generate orthogonal matrix in vl
! (workspace: need 3*n-1, prefer 2*n+(n-1)*nb)
call stdlib${ii}$_${ri}$orghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,work( iwrk ), &
lwork-iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_${ri}$lacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_${ri}$lacpy( 'L', n, n, a, lda, vr, ldvr )
! generate orthogonal matrix in vr
! (workspace: need 3*n-1, prefer 2*n+(n-1)*nb)
call stdlib${ii}$_${ri}$orghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
else
! compute eigenvalues only
! (workspace: need n+1, prefer n+hswork (see comments) )
iwrk = itau
call stdlib${ii}$_${ri}$hseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_${ri}$hseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (workspace: need 4*n, prefer n + n + 2*n*nb)
call stdlib${ii}$_${ri}$trevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1, ierr )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
! (workspace: need n)
call stdlib${ii}$_${ri}$gebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_${ri}$nrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_${ri}$scal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_${ri}$lapy2( stdlib${ii}$_${ri}$nrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_${ri}$nrm2( n, &
vl( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_${ri}$scal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_${ri}$scal( n, scl, vl( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( iwrk+k-1 ) = vl( k, i )**2_${ik}$ + vl( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_i${ri}$amax( n, work( iwrk ), 1_${ik}$ )
call stdlib${ii}$_${ri}$lartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
call stdlib${ii}$_${ri}$rot( n, vl( 1_${ik}$, i ), 1_${ik}$, vl( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vl( k, i+1 ) = zero
end if
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
! (workspace: need n)
call stdlib${ii}$_${ri}$gebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_${ri}$nrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_${ri}$scal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_${ri}$lapy2( stdlib${ii}$_${ri}$nrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_${ri}$nrm2( n, &
vr( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_${ri}$scal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_${ri}$scal( n, scl, vr( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( iwrk+k-1 ) = vr( k, i )**2_${ik}$ + vr( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_i${ri}$amax( n, work( iwrk ), 1_${ik}$ )
call stdlib${ii}$_${ri}$lartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
call stdlib${ii}$_${ri}$rot( n, vr( 1_${ik}$, i ), 1_${ik}$, vr( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vr( k, i+1 ) = zero
end if
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_${ri}$lascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wr( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
call stdlib${ii}$_${ri}$lascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wi( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
if( info>0_${ik}$ ) then
call stdlib${ii}$_${ri}$lascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wr, n,ierr )
call stdlib${ii}$_${ri}$lascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wi, n,ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_${ri}$geev
#:endif
#:endfor
module subroutine stdlib${ii}$_cgeev( jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr,work, lwork, rwork, &
!! CGEEV computes for an N-by-N complex nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobvl, jobvr
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
! Array Arguments
real(sp), intent(out) :: rwork(*)
complex(sp), intent(inout) :: a(lda,*)
complex(sp), intent(out) :: vl(ldvl,*), vr(ldvr,*), w(*), work(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr
character :: side
integer(${ik}$) :: hswork, i, ibal, ierr, ihi, ilo, irwork, itau, iwrk, k, lwork_trevc, &
maxwrk, minwrk, nout
real(sp) :: anrm, bignum, cscale, eps, scl, smlnum
complex(sp) :: tmp
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(sp) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -1_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -8_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -10_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! cworkspace refers to complex workspace, and rworkspace to real
! workspace. nb refers to the optimal block size for the
! immediately following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_chseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'CGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
minwrk = 2_${ik}$*n
if( wantvl ) then
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'CUNGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
call stdlib${ii}$_ctrevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, rwork, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_chseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, w, vl, ldvl,work, -1_${ik}$, info )
else if( wantvr ) then
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'CUNGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
call stdlib${ii}$_ctrevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, rwork, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_chseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, w, vr, ldvr,work, -1_${ik}$, info )
else
call stdlib${ii}$_chseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, w, vr, ldvr,work, -1_${ik}$, info )
end if
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, hswork, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGEEV ', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_slamch( 'P' )
smlnum = stdlib${ii}$_slamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_slabad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
anrm = stdlib${ii}$_clange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_clascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix
! (cworkspace: none)
! (rworkspace: need n)
ibal = 1_${ik}$
call stdlib${ii}$_cgebal( 'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
! reduce to upper hessenberg form
! (cworkspace: need 2*n, prefer n+n*nb)
! (rworkspace: none)
itau = 1_${ik}$
iwrk = itau + n
call stdlib${ii}$_cgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_clacpy( 'L', n, n, a, lda, vl, ldvl )
! generate unitary matrix in vl
! (cworkspace: need 2*n-1, prefer n+(n-1)*nb)
! (rworkspace: none)
call stdlib${ii}$_cunghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_chseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vl, ldvl,work( iwrk ), lwork-&
iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_clacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_clacpy( 'L', n, n, a, lda, vr, ldvr )
! generate unitary matrix in vr
! (cworkspace: need 2*n-1, prefer n+(n-1)*nb)
! (rworkspace: none)
call stdlib${ii}$_cunghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_chseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vr, ldvr,work( iwrk ), lwork-&
iwrk+1, info )
else
! compute eigenvalues only
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_chseqr( 'E', 'N', n, ilo, ihi, a, lda, w, vr, ldvr,work( iwrk ), lwork-&
iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_chseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (cworkspace: need 2*n, prefer n + 2*n*nb)
! (rworkspace: need 2*n)
irwork = ibal + n
call stdlib${ii}$_ctrevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1,rwork( irwork ), n, ierr )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
! (cworkspace: none)
! (rworkspace: need n)
call stdlib${ii}$_cgebak( 'B', 'L', n, ilo, ihi, rwork( ibal ), n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
scl = one / stdlib${ii}$_scnrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_csscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
do k = 1, n
rwork( irwork+k-1 ) = real( vl( k, i ),KIND=sp)**2_${ik}$ +aimag( vl( k, i ) )&
**2_${ik}$
end do
k = stdlib${ii}$_isamax( n, rwork( irwork ), 1_${ik}$ )
tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
call stdlib${ii}$_cscal( n, tmp, vl( 1_${ik}$, i ), 1_${ik}$ )
vl( k, i ) = cmplx( real( vl( k, i ),KIND=sp), zero,KIND=sp)
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
! (cworkspace: none)
! (rworkspace: need n)
call stdlib${ii}$_cgebak( 'B', 'R', n, ilo, ihi, rwork( ibal ), n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
scl = one / stdlib${ii}$_scnrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_csscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
do k = 1, n
rwork( irwork+k-1 ) = real( vr( k, i ),KIND=sp)**2_${ik}$ +aimag( vr( k, i ) )&
**2_${ik}$
end do
k = stdlib${ii}$_isamax( n, rwork( irwork ), 1_${ik}$ )
tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
call stdlib${ii}$_cscal( n, tmp, vr( 1_${ik}$, i ), 1_${ik}$ )
vr( k, i ) = cmplx( real( vr( k, i ),KIND=sp), zero,KIND=sp)
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_clascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, w( info+1 ),max( n-info, 1_${ik}$ )&
, ierr )
if( info>0_${ik}$ ) then
call stdlib${ii}$_clascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, w, n, ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_cgeev
module subroutine stdlib${ii}$_zgeev( jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr,work, lwork, rwork, &
!! ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobvl, jobvr
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
! Array Arguments
real(dp), intent(out) :: rwork(*)
complex(dp), intent(inout) :: a(lda,*)
complex(dp), intent(out) :: vl(ldvl,*), vr(ldvr,*), w(*), work(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr
character :: side
integer(${ik}$) :: hswork, i, ibal, ierr, ihi, ilo, irwork, itau, iwrk, k, lwork_trevc, &
maxwrk, minwrk, nout
real(dp) :: anrm, bignum, cscale, eps, scl, smlnum
complex(dp) :: tmp
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(dp) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -1_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -8_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -10_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! cworkspace refers to complex workspace, and rworkspace to real
! workspace. nb refers to the optimal block size for the
! immediately following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_zhseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
minwrk = 2_${ik}$*n
if( wantvl ) then
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
call stdlib${ii}$_ztrevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, rwork, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_zhseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, w, vl, ldvl,work, -1_${ik}$, info )
else if( wantvr ) then
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
call stdlib${ii}$_ztrevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, rwork, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_zhseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, w, vr, ldvr,work, -1_${ik}$, info )
else
call stdlib${ii}$_zhseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, w, vr, ldvr,work, -1_${ik}$, info )
end if
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, hswork, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGEEV ', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_dlamch( 'P' )
smlnum = stdlib${ii}$_dlamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_dlabad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
anrm = stdlib${ii}$_zlange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_zlascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix
! (cworkspace: none)
! (rworkspace: need n)
ibal = 1_${ik}$
call stdlib${ii}$_zgebal( 'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
! reduce to upper hessenberg form
! (cworkspace: need 2*n, prefer n+n*nb)
! (rworkspace: none)
itau = 1_${ik}$
iwrk = itau + n
call stdlib${ii}$_zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_zlacpy( 'L', n, n, a, lda, vl, ldvl )
! generate unitary matrix in vl
! (cworkspace: need 2*n-1, prefer n+(n-1)*nb)
! (rworkspace: none)
call stdlib${ii}$_zunghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_zhseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vl, ldvl,work( iwrk ), lwork-&
iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_zlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_zlacpy( 'L', n, n, a, lda, vr, ldvr )
! generate unitary matrix in vr
! (cworkspace: need 2*n-1, prefer n+(n-1)*nb)
! (rworkspace: none)
call stdlib${ii}$_zunghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_zhseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vr, ldvr,work( iwrk ), lwork-&
iwrk+1, info )
else
! compute eigenvalues only
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_zhseqr( 'E', 'N', n, ilo, ihi, a, lda, w, vr, ldvr,work( iwrk ), lwork-&
iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_zhseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (cworkspace: need 2*n, prefer n + 2*n*nb)
! (rworkspace: need 2*n)
irwork = ibal + n
call stdlib${ii}$_ztrevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1,rwork( irwork ), n, ierr )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
! (cworkspace: none)
! (rworkspace: need n)
call stdlib${ii}$_zgebak( 'B', 'L', n, ilo, ihi, rwork( ibal ), n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
scl = one / stdlib${ii}$_dznrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_zdscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
do k = 1, n
rwork( irwork+k-1 ) = real( vl( k, i ),KIND=dp)**2_${ik}$ +aimag( vl( k, i ) )&
**2_${ik}$
end do
k = stdlib${ii}$_idamax( n, rwork( irwork ), 1_${ik}$ )
tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
call stdlib${ii}$_zscal( n, tmp, vl( 1_${ik}$, i ), 1_${ik}$ )
vl( k, i ) = cmplx( real( vl( k, i ),KIND=dp), zero,KIND=dp)
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
! (cworkspace: none)
! (rworkspace: need n)
call stdlib${ii}$_zgebak( 'B', 'R', n, ilo, ihi, rwork( ibal ), n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
scl = one / stdlib${ii}$_dznrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_zdscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
do k = 1, n
rwork( irwork+k-1 ) = real( vr( k, i ),KIND=dp)**2_${ik}$ +aimag( vr( k, i ) )&
**2_${ik}$
end do
k = stdlib${ii}$_idamax( n, rwork( irwork ), 1_${ik}$ )
tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
call stdlib${ii}$_zscal( n, tmp, vr( 1_${ik}$, i ), 1_${ik}$ )
vr( k, i ) = cmplx( real( vr( k, i ),KIND=dp), zero,KIND=dp)
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_zlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, w( info+1 ),max( n-info, 1_${ik}$ )&
, ierr )
if( info>0_${ik}$ ) then
call stdlib${ii}$_zlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, w, n, ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_zgeev
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$geev( jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr,work, lwork, rwork, &
!! ZGEEV: computes for an N-by-N complex nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: jobvl, jobvr
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
! Array Arguments
real(${ck}$), intent(out) :: rwork(*)
complex(${ck}$), intent(inout) :: a(lda,*)
complex(${ck}$), intent(out) :: vl(ldvl,*), vr(ldvr,*), w(*), work(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr
character :: side
integer(${ik}$) :: hswork, i, ibal, ierr, ihi, ilo, irwork, itau, iwrk, k, lwork_trevc, &
maxwrk, minwrk, nout
real(${ck}$) :: anrm, bignum, cscale, eps, scl, smlnum
complex(${ck}$) :: tmp
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(${ck}$) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -1_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -8_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -10_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! cworkspace refers to complex workspace, and rworkspace to real
! workspace. nb refers to the optimal block size for the
! immediately following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_${ci}$hseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
minwrk = 2_${ik}$*n
if( wantvl ) then
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
call stdlib${ii}$_${ci}$trevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, rwork, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_${ci}$hseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, w, vl, ldvl,work, -1_${ik}$, info )
else if( wantvr ) then
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
call stdlib${ii}$_${ci}$trevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, rwork, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_${ci}$hseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, w, vr, ldvr,work, -1_${ik}$, info )
else
call stdlib${ii}$_${ci}$hseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, w, vr, ldvr,work, -1_${ik}$, info )
end if
hswork = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, hswork, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGEEV ', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_${c2ri(ci)}$lamch( 'P' )
smlnum = stdlib${ii}$_${c2ri(ci)}$lamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_${c2ri(ci)}$labad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
anrm = stdlib${ii}$_${ci}$lange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_${ci}$lascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix
! (cworkspace: none)
! (rworkspace: need n)
ibal = 1_${ik}$
call stdlib${ii}$_${ci}$gebal( 'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
! reduce to upper hessenberg form
! (cworkspace: need 2*n, prefer n+n*nb)
! (rworkspace: none)
itau = 1_${ik}$
iwrk = itau + n
call stdlib${ii}$_${ci}$gehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_${ci}$lacpy( 'L', n, n, a, lda, vl, ldvl )
! generate unitary matrix in vl
! (cworkspace: need 2*n-1, prefer n+(n-1)*nb)
! (rworkspace: none)
call stdlib${ii}$_${ci}$unghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_${ci}$hseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vl, ldvl,work( iwrk ), lwork-&
iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_${ci}$lacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_${ci}$lacpy( 'L', n, n, a, lda, vr, ldvr )
! generate unitary matrix in vr
! (cworkspace: need 2*n-1, prefer n+(n-1)*nb)
! (rworkspace: none)
call stdlib${ii}$_${ci}$unghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_${ci}$hseqr( 'S', 'V', n, ilo, ihi, a, lda, w, vr, ldvr,work( iwrk ), lwork-&
iwrk+1, info )
else
! compute eigenvalues only
! (cworkspace: need 1, prefer hswork (see comments) )
! (rworkspace: none)
iwrk = itau
call stdlib${ii}$_${ci}$hseqr( 'E', 'N', n, ilo, ihi, a, lda, w, vr, ldvr,work( iwrk ), lwork-&
iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_${ci}$hseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (cworkspace: need 2*n, prefer n + 2*n*nb)
! (rworkspace: need 2*n)
irwork = ibal + n
call stdlib${ii}$_${ci}$trevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1,rwork( irwork ), n, ierr )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
! (cworkspace: none)
! (rworkspace: need n)
call stdlib${ii}$_${ci}$gebak( 'B', 'L', n, ilo, ihi, rwork( ibal ), n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
scl = one / stdlib${ii}$_${c2ri(ci)}$znrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_${ci}$dscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
do k = 1, n
rwork( irwork+k-1 ) = real( vl( k, i ),KIND=${ck}$)**2_${ik}$ +aimag( vl( k, i ) )&
**2_${ik}$
end do
k = stdlib${ii}$_i${c2ri(ci)}$amax( n, rwork( irwork ), 1_${ik}$ )
tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
call stdlib${ii}$_${ci}$scal( n, tmp, vl( 1_${ik}$, i ), 1_${ik}$ )
vl( k, i ) = cmplx( real( vl( k, i ),KIND=${ck}$), zero,KIND=${ck}$)
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
! (cworkspace: none)
! (rworkspace: need n)
call stdlib${ii}$_${ci}$gebak( 'B', 'R', n, ilo, ihi, rwork( ibal ), n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
scl = one / stdlib${ii}$_${c2ri(ci)}$znrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_${ci}$dscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
do k = 1, n
rwork( irwork+k-1 ) = real( vr( k, i ),KIND=${ck}$)**2_${ik}$ +aimag( vr( k, i ) )&
**2_${ik}$
end do
k = stdlib${ii}$_i${c2ri(ci)}$amax( n, rwork( irwork ), 1_${ik}$ )
tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
call stdlib${ii}$_${ci}$scal( n, tmp, vr( 1_${ik}$, i ), 1_${ik}$ )
vr( k, i ) = cmplx( real( vr( k, i ),KIND=${ck}$), zero,KIND=${ck}$)
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_${ci}$lascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, w( info+1 ),max( n-info, 1_${ik}$ )&
, ierr )
if( info>0_${ik}$ ) then
call stdlib${ii}$_${ci}$lascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, w, n, ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_${ci}$geev
#:endif
#:endfor
module subroutine stdlib${ii}$_sgeevx( balanc, jobvl, jobvr, sense, n, a, lda, wr, wi,vl, ldvl, vr, ldvr, &
!! SGEEVX computes for an N-by-N real nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! Optionally also, it computes a balancing transformation to improve
!! the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
!! SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
!! (RCONDE), and reciprocal condition numbers for the right
!! eigenvectors (RCONDV).
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate-transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
!! Balancing a matrix means permuting the rows and columns to make it
!! more nearly upper triangular, and applying a diagonal similarity
!! transformation D * A * D**(-1), where D is a diagonal matrix, to
!! make its rows and columns closer in norm and the condition numbers
!! of its eigenvalues and eigenvectors smaller. The computed
!! reciprocal condition numbers correspond to the balanced matrix.
!! Permuting rows and columns will not change the condition numbers
!! (in exact arithmetic) but diagonal scaling will. For further
!! explanation of balancing, see section 4.10.2_sp of the LAPACK
!! Users' Guide.
ilo, ihi, scale, abnrm,rconde, rcondv, work, lwork, iwork, info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: balanc, jobvl, jobvr, sense
integer(${ik}$), intent(out) :: ihi, ilo, info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
real(sp), intent(out) :: abnrm
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(inout) :: a(lda,*)
real(sp), intent(out) :: rconde(*), rcondv(*), scale(*), vl(ldvl,*), vr(ldvr,*), wi(*),&
work(*), wr(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr, wntsnb, wntsne, wntsnn, wntsnv
character :: job, side
integer(${ik}$) :: hswork, i, icond, ierr, itau, iwrk, k, lwork_trevc, maxwrk, minwrk, &
nout
real(sp) :: anrm, bignum, cs, cscale, eps, r, scl, smlnum, sn
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(sp) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
wntsnn = stdlib_lsame( sense, 'N' )
wntsne = stdlib_lsame( sense, 'E' )
wntsnv = stdlib_lsame( sense, 'V' )
wntsnb = stdlib_lsame( sense, 'B' )
if( .not.( stdlib_lsame( balanc, 'N' ) .or. stdlib_lsame( balanc, 'S' ).or. &
stdlib_lsame( balanc, 'P' ) .or. stdlib_lsame( balanc, 'B' ) ) )then
info = -1_${ik}$
else if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -2_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -3_${ik}$
else if( .not.( wntsnn .or. wntsne .or. wntsnb .or. wntsnv ) .or.( ( wntsne .or. &
wntsnb ) .and. .not.( wantvl .and.wantvr ) ) ) then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -11_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -13_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! nb refers to the optimal block size for the immediately
! following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_shseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'SGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
if( wantvl ) then
call stdlib${ii}$_strevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_shseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vl, ldvl,work, -1_${ik}$, &
info )
else if( wantvr ) then
call stdlib${ii}$_strevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_shseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
else
if( wntsnn ) then
call stdlib${ii}$_shseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr,ldvr, work, -1_${ik}$, &
info )
else
call stdlib${ii}$_shseqr( 'S', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr,ldvr, work, -1_${ik}$, &
info )
end if
end if
hswork = int( work(1_${ik}$),KIND=${ik}$)
if( ( .not.wantvl ) .and. ( .not.wantvr ) ) then
minwrk = 2_${ik}$*n
if( .not.wntsnn )minwrk = max( minwrk, n*n+6*n )
maxwrk = max( maxwrk, hswork )
if( .not.wntsnn )maxwrk = max( maxwrk, n*n + 6_${ik}$*n )
else
minwrk = 3_${ik}$*n
if( ( .not.wntsnn ) .and. ( .not.wntsne ) )minwrk = max( minwrk, n*n + 6_${ik}$*n )
maxwrk = max( maxwrk, hswork )
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'SORGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
if( ( .not.wntsnn ) .and. ( .not.wntsne ) )maxwrk = max( maxwrk, n*n + 6_${ik}$*n )
maxwrk = max( maxwrk, 3_${ik}$*n )
end if
maxwrk = max( maxwrk, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -21_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGEEVX', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_slamch( 'P' )
smlnum = stdlib${ii}$_slamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_slabad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
icond = 0_${ik}$
anrm = stdlib${ii}$_slange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix and compute abnrm
call stdlib${ii}$_sgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
abnrm = stdlib${ii}$_slange( '1', n, n, a, lda, dum )
if( scalea ) then
dum( 1_${ik}$ ) = abnrm
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, 1_${ik}$, 1_${ik}$, dum, 1_${ik}$, ierr )
abnrm = dum( 1_${ik}$ )
end if
! reduce to upper hessenberg form
! (workspace: need 2*n, prefer n+n*nb)
itau = 1_${ik}$
iwrk = itau + n
call stdlib${ii}$_sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_slacpy( 'L', n, n, a, lda, vl, ldvl )
! generate orthogonal matrix in vl
! (workspace: need 2*n-1, prefer n+(n-1)*nb)
call stdlib${ii}$_sorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,work( iwrk ), &
lwork-iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_slacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_slacpy( 'L', n, n, a, lda, vr, ldvr )
! generate orthogonal matrix in vr
! (workspace: need 2*n-1, prefer n+(n-1)*nb)
call stdlib${ii}$_sorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_shseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
else
! compute eigenvalues only
! if condition numbers desired, compute schur form
if( wntsnn ) then
job = 'E'
else
job = 'S'
end if
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_shseqr( job, 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_shseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (workspace: need 3*n, prefer n + 2*n*nb)
call stdlib${ii}$_strevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1, ierr )
end if
! compute condition numbers if desired
! (workspace: need n*n+6*n unless sense = 'e')
if( .not.wntsnn ) then
call stdlib${ii}$_strsna( sense, 'A', select, n, a, lda, vl, ldvl, vr, ldvr,rconde, &
rcondv, n, nout, work( iwrk ), n, iwork,icond )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
call stdlib${ii}$_sgebak( balanc, 'L', n, ilo, ihi, scale, n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_snrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, &
vl( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_sscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vl( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( k ) = vl( k, i )**2_${ik}$ + vl( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_isamax( n, work, 1_${ik}$ )
call stdlib${ii}$_slartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
call stdlib${ii}$_srot( n, vl( 1_${ik}$, i ), 1_${ik}$, vl( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vl( k, i+1 ) = zero
end if
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
call stdlib${ii}$_sgebak( balanc, 'R', n, ilo, ihi, scale, n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_snrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, &
vr( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_sscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_sscal( n, scl, vr( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( k ) = vr( k, i )**2_${ik}$ + vr( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_isamax( n, work, 1_${ik}$ )
call stdlib${ii}$_slartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
call stdlib${ii}$_srot( n, vr( 1_${ik}$, i ), 1_${ik}$, vr( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vr( k, i+1 ) = zero
end if
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wr( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wi( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
if( info==0_${ik}$ ) then
if( ( wntsnv .or. wntsnb ) .and. icond==0_${ik}$ )call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale,&
anrm, n, 1_${ik}$, rcondv, n,ierr )
else
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wr, n,ierr )
call stdlib${ii}$_slascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wi, n,ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_sgeevx
module subroutine stdlib${ii}$_dgeevx( balanc, jobvl, jobvr, sense, n, a, lda, wr, wi,vl, ldvl, vr, ldvr, &
!! DGEEVX computes for an N-by-N real nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! Optionally also, it computes a balancing transformation to improve
!! the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
!! SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
!! (RCONDE), and reciprocal condition numbers for the right
!! eigenvectors (RCONDV).
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate-transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
!! Balancing a matrix means permuting the rows and columns to make it
!! more nearly upper triangular, and applying a diagonal similarity
!! transformation D * A * D**(-1), where D is a diagonal matrix, to
!! make its rows and columns closer in norm and the condition numbers
!! of its eigenvalues and eigenvectors smaller. The computed
!! reciprocal condition numbers correspond to the balanced matrix.
!! Permuting rows and columns will not change the condition numbers
!! (in exact arithmetic) but diagonal scaling will. For further
!! explanation of balancing, see section 4.10.2_dp of the LAPACK
!! Users' Guide.
ilo, ihi, scale, abnrm,rconde, rcondv, work, lwork, iwork, info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: balanc, jobvl, jobvr, sense
integer(${ik}$), intent(out) :: ihi, ilo, info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
real(dp), intent(out) :: abnrm
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(inout) :: a(lda,*)
real(dp), intent(out) :: rconde(*), rcondv(*), scale(*), vl(ldvl,*), vr(ldvr,*), wi(*),&
work(*), wr(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr, wntsnb, wntsne, wntsnn, wntsnv
character :: job, side
integer(${ik}$) :: hswork, i, icond, ierr, itau, iwrk, k, lwork_trevc, maxwrk, minwrk, &
nout
real(dp) :: anrm, bignum, cs, cscale, eps, r, scl, smlnum, sn
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(dp) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
wntsnn = stdlib_lsame( sense, 'N' )
wntsne = stdlib_lsame( sense, 'E' )
wntsnv = stdlib_lsame( sense, 'V' )
wntsnb = stdlib_lsame( sense, 'B' )
if( .not.( stdlib_lsame( balanc, 'N' ) .or. stdlib_lsame( balanc, 'S' ).or. &
stdlib_lsame( balanc, 'P' ) .or. stdlib_lsame( balanc, 'B' ) ) )then
info = -1_${ik}$
else if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -2_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -3_${ik}$
else if( .not.( wntsnn .or. wntsne .or. wntsnb .or. wntsnv ) .or.( ( wntsne .or. &
wntsnb ) .and. .not.( wantvl .and.wantvr ) ) ) then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -11_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -13_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! nb refers to the optimal block size for the immediately
! following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_dhseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
if( wantvl ) then
call stdlib${ii}$_dtrevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_dhseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vl, ldvl,work, -1_${ik}$, &
info )
else if( wantvr ) then
call stdlib${ii}$_dtrevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_dhseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
else
if( wntsnn ) then
call stdlib${ii}$_dhseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr,ldvr, work, -1_${ik}$, &
info )
else
call stdlib${ii}$_dhseqr( 'S', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr,ldvr, work, -1_${ik}$, &
info )
end if
end if
hswork = int( work(1_${ik}$),KIND=${ik}$)
if( ( .not.wantvl ) .and. ( .not.wantvr ) ) then
minwrk = 2_${ik}$*n
if( .not.wntsnn )minwrk = max( minwrk, n*n+6*n )
maxwrk = max( maxwrk, hswork )
if( .not.wntsnn )maxwrk = max( maxwrk, n*n + 6_${ik}$*n )
else
minwrk = 3_${ik}$*n
if( ( .not.wntsnn ) .and. ( .not.wntsne ) )minwrk = max( minwrk, n*n + 6_${ik}$*n )
maxwrk = max( maxwrk, hswork )
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'DORGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
if( ( .not.wntsnn ) .and. ( .not.wntsne ) )maxwrk = max( maxwrk, n*n + 6_${ik}$*n )
maxwrk = max( maxwrk, 3_${ik}$*n )
end if
maxwrk = max( maxwrk, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -21_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGEEVX', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_dlamch( 'P' )
smlnum = stdlib${ii}$_dlamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_dlabad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
icond = 0_${ik}$
anrm = stdlib${ii}$_dlange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix and compute abnrm
call stdlib${ii}$_dgebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
abnrm = stdlib${ii}$_dlange( '1', n, n, a, lda, dum )
if( scalea ) then
dum( 1_${ik}$ ) = abnrm
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, 1_${ik}$, 1_${ik}$, dum, 1_${ik}$, ierr )
abnrm = dum( 1_${ik}$ )
end if
! reduce to upper hessenberg form
! (workspace: need 2*n, prefer n+n*nb)
itau = 1_${ik}$
iwrk = itau + n
call stdlib${ii}$_dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_dlacpy( 'L', n, n, a, lda, vl, ldvl )
! generate orthogonal matrix in vl
! (workspace: need 2*n-1, prefer n+(n-1)*nb)
call stdlib${ii}$_dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,work( iwrk ), &
lwork-iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_dlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_dlacpy( 'L', n, n, a, lda, vr, ldvr )
! generate orthogonal matrix in vr
! (workspace: need 2*n-1, prefer n+(n-1)*nb)
call stdlib${ii}$_dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
else
! compute eigenvalues only
! if condition numbers desired, compute schur form
if( wntsnn ) then
job = 'E'
else
job = 'S'
end if
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_dhseqr( job, 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
end if
! if info /= 0 from stdlib${ii}$_dhseqr, then quit
if( info/=0 )go to 50
if( wantvl .or. wantvr ) then
! compute left and/or right eigenvectors
! (workspace: need 3*n, prefer n + 2*n*nb)
call stdlib${ii}$_dtrevc3( side, 'B', select, n, a, lda, vl, ldvl, vr, ldvr,n, nout, work(&
iwrk ), lwork-iwrk+1, ierr )
end if
! compute condition numbers if desired
! (workspace: need n*n+6*n unless sense = 'e')
if( .not.wntsnn ) then
call stdlib${ii}$_dtrsna( sense, 'A', select, n, a, lda, vl, ldvl, vr, ldvr,rconde, &
rcondv, n, nout, work( iwrk ), n, iwork,icond )
end if
if( wantvl ) then
! undo balancing of left eigenvectors
call stdlib${ii}$_dgebak( balanc, 'L', n, ilo, ihi, scale, n, vl, ldvl,ierr )
! normalize left eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, &
vl( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_dscal( n, scl, vl( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vl( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( k ) = vl( k, i )**2_${ik}$ + vl( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_idamax( n, work, 1_${ik}$ )
call stdlib${ii}$_dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
call stdlib${ii}$_drot( n, vl( 1_${ik}$, i ), 1_${ik}$, vl( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vl( k, i+1 ) = zero
end if
end do
end if
if( wantvr ) then
! undo balancing of right eigenvectors
call stdlib${ii}$_dgebak( balanc, 'R', n, ilo, ihi, scale, n, vr, ldvr,ierr )
! normalize right eigenvectors and make largest component real
do i = 1, n
if( wi( i )==zero ) then
scl = one / stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
else if( wi( i )>zero ) then
scl = one / stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, i ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, &
vr( 1_${ik}$, i+1 ), 1_${ik}$ ) )
call stdlib${ii}$_dscal( n, scl, vr( 1_${ik}$, i ), 1_${ik}$ )
call stdlib${ii}$_dscal( n, scl, vr( 1_${ik}$, i+1 ), 1_${ik}$ )
do k = 1, n
work( k ) = vr( k, i )**2_${ik}$ + vr( k, i+1 )**2_${ik}$
end do
k = stdlib${ii}$_idamax( n, work, 1_${ik}$ )
call stdlib${ii}$_dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
call stdlib${ii}$_drot( n, vr( 1_${ik}$, i ), 1_${ik}$, vr( 1_${ik}$, i+1 ), 1_${ik}$, cs, sn )
vr( k, i+1 ) = zero
end if
end do
end if
! undo scaling if necessary
50 continue
if( scalea ) then
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wr( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, n-info, 1_${ik}$, wi( info+1 ),max( n-info, 1_${ik}$ &
), ierr )
if( info==0_${ik}$ ) then
if( ( wntsnv .or. wntsnb ) .and. icond==0_${ik}$ )call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale,&
anrm, n, 1_${ik}$, rcondv, n,ierr )
else
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wr, n,ierr )
call stdlib${ii}$_dlascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, ilo-1, 1_${ik}$, wi, n,ierr )
end if
end if
work( 1_${ik}$ ) = maxwrk
return
end subroutine stdlib${ii}$_dgeevx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$geevx( balanc, jobvl, jobvr, sense, n, a, lda, wr, wi,vl, ldvl, vr, ldvr, &
!! DGEEVX: computes for an N-by-N real nonsymmetric matrix A, the
!! eigenvalues and, optionally, the left and/or right eigenvectors.
!! Optionally also, it computes a balancing transformation to improve
!! the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
!! SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues
!! (RCONDE), and reciprocal condition numbers for the right
!! eigenvectors (RCONDV).
!! The right eigenvector v(j) of A satisfies
!! A * v(j) = lambda(j) * v(j)
!! where lambda(j) is its eigenvalue.
!! The left eigenvector u(j) of A satisfies
!! u(j)**H * A = lambda(j) * u(j)**H
!! where u(j)**H denotes the conjugate-transpose of u(j).
!! The computed eigenvectors are normalized to have Euclidean norm
!! equal to 1 and largest component real.
!! Balancing a matrix means permuting the rows and columns to make it
!! more nearly upper triangular, and applying a diagonal similarity
!! transformation D * A * D**(-1), where D is a diagonal matrix, to
!! make its rows and columns closer in norm and the condition numbers
!! of its eigenvalues and eigenvectors smaller. The computed
!! reciprocal condition numbers correspond to the balanced matrix.
!! Permuting rows and columns will not change the condition numbers
!! (in exact arithmetic) but diagonal scaling will. For further
!! explanation of balancing, see section 4.10.2_${rk}$ of the LAPACK
!! Users' Guide.
ilo, ihi, scale, abnrm,rconde, rcondv, work, lwork, iwork, info )
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: balanc, jobvl, jobvr, sense
integer(${ik}$), intent(out) :: ihi, ilo, info
integer(${ik}$), intent(in) :: lda, ldvl, ldvr, lwork, n
real(${rk}$), intent(out) :: abnrm
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(inout) :: a(lda,*)
real(${rk}$), intent(out) :: rconde(*), rcondv(*), scale(*), vl(ldvl,*), vr(ldvr,*), wi(*),&
work(*), wr(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery, scalea, wantvl, wantvr, wntsnb, wntsne, wntsnn, wntsnv
character :: job, side
integer(${ik}$) :: hswork, i, icond, ierr, itau, iwrk, k, lwork_trevc, maxwrk, minwrk, &
nout
real(${rk}$) :: anrm, bignum, cs, cscale, eps, r, scl, smlnum, sn
! Local Arrays
logical(lk) :: select(1_${ik}$)
real(${rk}$) :: dum(1_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input arguments
info = 0_${ik}$
lquery = ( lwork==-1_${ik}$ )
wantvl = stdlib_lsame( jobvl, 'V' )
wantvr = stdlib_lsame( jobvr, 'V' )
wntsnn = stdlib_lsame( sense, 'N' )
wntsne = stdlib_lsame( sense, 'E' )
wntsnv = stdlib_lsame( sense, 'V' )
wntsnb = stdlib_lsame( sense, 'B' )
if( .not.( stdlib_lsame( balanc, 'N' ) .or. stdlib_lsame( balanc, 'S' ).or. &
stdlib_lsame( balanc, 'P' ) .or. stdlib_lsame( balanc, 'B' ) ) )then
info = -1_${ik}$
else if( ( .not.wantvl ) .and. ( .not.stdlib_lsame( jobvl, 'N' ) ) ) then
info = -2_${ik}$
else if( ( .not.wantvr ) .and. ( .not.stdlib_lsame( jobvr, 'N' ) ) ) then
info = -3_${ik}$
else if( .not.( wntsnn .or. wntsne .or. wntsnb .or. wntsnv ) .or.( ( wntsne .or. &
wntsnb ) .and. .not.( wantvl .and.wantvr ) ) ) then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldvl<1_${ik}$ .or. ( wantvl .and. ldvl<n ) ) then
info = -11_${ik}$
else if( ldvr<1_${ik}$ .or. ( wantvr .and. ldvr<n ) ) then
info = -13_${ik}$
end if
! compute workspace
! (note: comments in the code beginning "workspace:" describe the
! minimal amount of workspace needed at that point in the code,
! as well as the preferred amount for good performance.
! nb refers to the optimal block size for the immediately
! following subroutine, as returned by stdlib${ii}$_ilaenv.
! hswork refers to the workspace preferred by stdlib${ii}$_${ri}$hseqr, as
! calculated below. hswork is computed assuming ilo=1 and ihi=n,
! the worst case.)
if( info==0_${ik}$ ) then
if( n==0_${ik}$ ) then
minwrk = 1_${ik}$
maxwrk = 1_${ik}$
else
maxwrk = n + n*stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEHRD', ' ', n, 1_${ik}$, n, 0_${ik}$ )
if( wantvl ) then
call stdlib${ii}$_${ri}$trevc3( 'L', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vl, ldvl,work, -1_${ik}$, &
info )
else if( wantvr ) then
call stdlib${ii}$_${ri}$trevc3( 'R', 'B', select, n, a, lda,vl, ldvl, vr, ldvr,n, nout, &
work, -1_${ik}$, ierr )
lwork_trevc = int( work(1_${ik}$),KIND=${ik}$)
maxwrk = max( maxwrk, n + lwork_trevc )
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, 1_${ik}$, n, a, lda, wr, wi, vr, ldvr,work, -1_${ik}$, &
info )
else
if( wntsnn ) then
call stdlib${ii}$_${ri}$hseqr( 'E', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr,ldvr, work, -1_${ik}$, &
info )
else
call stdlib${ii}$_${ri}$hseqr( 'S', 'N', n, 1_${ik}$, n, a, lda, wr, wi, vr,ldvr, work, -1_${ik}$, &
info )
end if
end if
hswork = int( work(1_${ik}$),KIND=${ik}$)
if( ( .not.wantvl ) .and. ( .not.wantvr ) ) then
minwrk = 2_${ik}$*n
if( .not.wntsnn )minwrk = max( minwrk, n*n+6*n )
maxwrk = max( maxwrk, hswork )
if( .not.wntsnn )maxwrk = max( maxwrk, n*n + 6_${ik}$*n )
else
minwrk = 3_${ik}$*n
if( ( .not.wntsnn ) .and. ( .not.wntsne ) )minwrk = max( minwrk, n*n + 6_${ik}$*n )
maxwrk = max( maxwrk, hswork )
maxwrk = max( maxwrk, n + ( n - 1_${ik}$ )*stdlib${ii}$_ilaenv( 1_${ik}$, 'DORGHR',' ', n, 1_${ik}$, n, -&
1_${ik}$ ) )
if( ( .not.wntsnn ) .and. ( .not.wntsne ) )maxwrk = max( maxwrk, n*n + 6_${ik}$*n )
maxwrk = max( maxwrk, 3_${ik}$*n )
end if
maxwrk = max( maxwrk, minwrk )
end if
work( 1_${ik}$ ) = maxwrk
if( lwork<minwrk .and. .not.lquery ) then
info = -21_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGEEVX', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! get machine constants
eps = stdlib${ii}$_${ri}$lamch( 'P' )
smlnum = stdlib${ii}$_${ri}$lamch( 'S' )
bignum = one / smlnum
call stdlib${ii}$_${ri}$labad( smlnum, bignum )
smlnum = sqrt( smlnum ) / eps
bignum = one / smlnum
! scale a if max element outside range [smlnum,bignum]
icond = 0_${ik}$
anrm = stdlib${ii}$_${ri}$lange( 'M', n, n, a, lda, dum )
scalea = .false.
if( anrm>zero .and. anrm<smlnum ) then
scalea = .true.
cscale = smlnum
else if( anrm>bignum ) then
scalea = .true.
cscale = bignum
end if
if( scalea )call stdlib${ii}$_${ri}$lascl( 'G', 0_${ik}$, 0_${ik}$, anrm, cscale, n, n, a, lda, ierr )
! balance the matrix and compute abnrm
call stdlib${ii}$_${ri}$gebal( balanc, n, a, lda, ilo, ihi, scale, ierr )
abnrm = stdlib${ii}$_${ri}$lange( '1', n, n, a, lda, dum )
if( scalea ) then
dum( 1_${ik}$ ) = abnrm
call stdlib${ii}$_${ri}$lascl( 'G', 0_${ik}$, 0_${ik}$, cscale, anrm, 1_${ik}$, 1_${ik}$, dum, 1_${ik}$, ierr )
abnrm = dum( 1_${ik}$ )
end if
! reduce to upper hessenberg form
! (workspace: need 2*n, prefer n+n*nb)
itau = 1_${ik}$
iwrk = itau + n
call stdlib${ii}$_${ri}$gehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),lwork-iwrk+1, ierr &
)
if( wantvl ) then
! want left eigenvectors
! copy householder vectors to vl
side = 'L'
call stdlib${ii}$_${ri}$lacpy( 'L', n, n, a, lda, vl, ldvl )
! generate orthogonal matrix in vl
! (workspace: need 2*n-1, prefer n+(n-1)*nb)
call stdlib${ii}$_${ri}$orghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vl
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,work( iwrk ), &
lwork-iwrk+1, info )
if( wantvr ) then
! want left and right eigenvectors
! copy schur vectors to vr
side = 'B'
call stdlib${ii}$_${ri}$lacpy( 'F', n, n, vl, ldvl, vr, ldvr )
end if
else if( wantvr ) then
! want right eigenvectors
! copy householder vectors to vr
side = 'R'
call stdlib${ii}$_${ri}$lacpy( 'L', n, n, a, lda, vr, ldvr )
! generate orthogonal matrix in vr
! (workspace: need 2*n-1, prefer n+(n-1)*nb)
call stdlib${ii}$_${ri}$orghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),lwork-iwrk+1, &
ierr )
! perform qr iteration, accumulating schur vectors in vr
! (workspace: need 1, prefer hswork (see comments) )
iwrk = itau
call stdlib${ii}$_${ri}$hseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,work( iwrk ), &
lwork-iwrk+1, info )
else
! compute eigenvalues only
! if condition numbers desired, compute schur form
if(