#:include "common.fypp" submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_gen2 implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES module subroutine stdlib${ii}$_shseqr( job, compz, n, ilo, ihi, h, ldh, wr, wi, z,ldz, work, lwork, info ) !! SHSEQR computes the eigenvalues of a Hessenberg matrix H !! and, optionally, the matrices T and Z from the Schur decomposition !! H = Z T Z**T, where T is an upper quasi-triangular matrix (the !! Schur form), and Z is the orthogonal matrix of Schur vectors. !! Optionally Z may be postmultiplied into an input orthogonal !! matrix Q so that this routine can give the Schur factorization !! of a matrix A which has been reduced to the Hessenberg form H !! by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. ! -- lapack computational 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 integer(${ik}$), intent(in) :: ihi, ilo, ldh, ldz, lwork, n integer(${ik}$), intent(out) :: info character, intent(in) :: compz, job ! Array Arguments real(sp), intent(inout) :: h(ldh,*), z(ldz,*) real(sp), intent(out) :: wi(*), work(*), wr(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: ntiny = 15_${ik}$ integer(${ik}$), parameter :: nl = 49_${ik}$ ! ==== matrices of order ntiny or smaller must be processed by ! . stdlib${ii}$_slahqr because of insufficient subdiagonal scratch space. ! . (this is a hard limit.) ==== ! ==== nl allocates some local workspace to help small matrices ! . through a rare stdlib${ii}$_slahqr failure. nl > ntiny = 15 is ! . required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom- ! . mended. (the default value of nmin is 75.) using nl = 49 ! . allows up to six simultaneous shifts and a 16-by-16 ! . deflation window. ==== ! Local Arrays real(sp) :: hl(nl,nl), workl(nl) ! Local Scalars integer(${ik}$) :: i, kbot, nmin logical(lk) :: initz, lquery, wantt, wantz ! Intrinsic Functions ! Executable Statements ! ==== decode and check the input parameters. ==== wantt = stdlib_lsame( job, 'S' ) initz = stdlib_lsame( compz, 'I' ) wantz = initz .or. stdlib_lsame( compz, 'V' ) work( 1_${ik}$ ) = real( max( 1_${ik}$, n ),KIND=sp) lquery = lwork==-1_${ik}$ info = 0_${ik}$ if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then info = -1_${ik}$ else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( ihi<min( ilo, n ) .or. ihi>n ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then info = -11_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -13_${ik}$ end if if( info/=0_${ik}$ ) then ! ==== quick return in case of invalid argument. ==== call stdlib${ii}$_xerbla( 'SHSEQR', -info ) return else if( n==0_${ik}$ ) then ! ==== quick return in case n = 0; nothing to do. ==== return else if( lquery ) then ! ==== quick return in case of a workspace query ==== call stdlib${ii}$_slaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & work, lwork, info ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=sp), work( 1_${ik}$ ) ) return else ! ==== copy eigenvalues isolated by stdlib${ii}$_sgebal ==== do i = 1, ilo - 1 wr( i ) = h( i, i ) wi( i ) = zero end do do i = ihi + 1, n wr( i ) = h( i, i ) wi( i ) = zero end do ! ==== initialize z, if requested ==== if( initz )call stdlib${ii}$_slaset( 'A', n, n, zero, one, z, ldz ) ! ==== quick return if possible ==== if( ilo==ihi ) then wr( ilo ) = h( ilo, ilo ) wi( ilo ) = zero return end if ! ==== stdlib${ii}$_slahqr/stdlib${ii}$_slaqr0 crossover point ==== nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'SHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork ) nmin = max( ntiny, nmin ) ! ==== stdlib${ii}$_slaqr0 for big matrices; stdlib${ii}$_slahqr for small ones ==== if( n>nmin ) then call stdlib${ii}$_slaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & work, lwork, info ) else ! ==== small matrix ==== call stdlib${ii}$_slahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & info ) if( info>0_${ik}$ ) then ! ==== a rare stdlib${ii}$_slahqr failure! stdlib${ii}$_slaqr0 sometimes succeeds ! . when stdlib${ii}$_slahqr fails. ==== kbot = info if( n>=nl ) then ! ==== larger matrices have enough subdiagonal scratch ! . space to call stdlib${ii}$_slaqr0 directly. ==== call stdlib${ii}$_slaqr0( wantt, wantz, n, ilo, kbot, h, ldh, wr,wi, ilo, ihi, z,& ldz, work, lwork, info ) else ! ==== tiny matrices don't have enough subdiagonal ! . scratch space to benefit from stdlib${ii}$_slaqr0. hence, ! . tiny matrices must be copied into a larger ! . array before calling stdlib${ii}$_slaqr0. ==== call stdlib${ii}$_slacpy( 'A', n, n, h, ldh, hl, nl ) hl( n+1, n ) = zero call stdlib${ii}$_slaset( 'A', nl, nl-n, zero, zero, hl( 1_${ik}$, n+1 ),nl ) call stdlib${ii}$_slaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, wr,wi, ilo, ihi, & z, ldz, workl, nl, info ) if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_slacpy( 'A', n, n, hl, nl, h, ldh ) end if end if end if ! ==== clear out the trash, if necessary. ==== if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_slaset( 'L', n-2, n-2, zero, zero,& h( 3_${ik}$, 1_${ik}$ ), ldh ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=sp), work( 1_${ik}$ ) ) end if end subroutine stdlib${ii}$_shseqr module subroutine stdlib${ii}$_dhseqr( job, compz, n, ilo, ihi, h, ldh, wr, wi, z,ldz, work, lwork, info ) !! DHSEQR computes the eigenvalues of a Hessenberg matrix H !! and, optionally, the matrices T and Z from the Schur decomposition !! H = Z T Z**T, where T is an upper quasi-triangular matrix (the !! Schur form), and Z is the orthogonal matrix of Schur vectors. !! Optionally Z may be postmultiplied into an input orthogonal !! matrix Q so that this routine can give the Schur factorization !! of a matrix A which has been reduced to the Hessenberg form H !! by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. ! -- lapack computational 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 integer(${ik}$), intent(in) :: ihi, ilo, ldh, ldz, lwork, n integer(${ik}$), intent(out) :: info character, intent(in) :: compz, job ! Array Arguments real(dp), intent(inout) :: h(ldh,*), z(ldz,*) real(dp), intent(out) :: wi(*), work(*), wr(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: ntiny = 15_${ik}$ integer(${ik}$), parameter :: nl = 49_${ik}$ ! ==== matrices of order ntiny or smaller must be processed by ! . stdlib${ii}$_dlahqr because of insufficient subdiagonal scratch space. ! . (this is a hard limit.) ==== ! ==== nl allocates some local workspace to help small matrices ! . through a rare stdlib${ii}$_dlahqr failure. nl > ntiny = 15 is ! . required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom- ! . mended. (the default value of nmin is 75.) using nl = 49 ! . allows up to six simultaneous shifts and a 16-by-16 ! . deflation window. ==== ! Local Arrays real(dp) :: hl(nl,nl), workl(nl) ! Local Scalars integer(${ik}$) :: i, kbot, nmin logical(lk) :: initz, lquery, wantt, wantz ! Intrinsic Functions ! Executable Statements ! ==== decode and check the input parameters. ==== wantt = stdlib_lsame( job, 'S' ) initz = stdlib_lsame( compz, 'I' ) wantz = initz .or. stdlib_lsame( compz, 'V' ) work( 1_${ik}$ ) = real( max( 1_${ik}$, n ),KIND=dp) lquery = lwork==-1_${ik}$ info = 0_${ik}$ if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then info = -1_${ik}$ else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( ihi<min( ilo, n ) .or. ihi>n ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then info = -11_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -13_${ik}$ end if if( info/=0_${ik}$ ) then ! ==== quick return in case of invalid argument. ==== call stdlib${ii}$_xerbla( 'DHSEQR', -info ) return else if( n==0_${ik}$ ) then ! ==== quick return in case n = 0; nothing to do. ==== return else if( lquery ) then ! ==== quick return in case of a workspace query ==== call stdlib${ii}$_dlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & work, lwork, info ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=dp), work( 1_${ik}$ ) ) return else ! ==== copy eigenvalues isolated by stdlib${ii}$_dgebal ==== do i = 1, ilo - 1 wr( i ) = h( i, i ) wi( i ) = zero end do do i = ihi + 1, n wr( i ) = h( i, i ) wi( i ) = zero end do ! ==== initialize z, if requested ==== if( initz )call stdlib${ii}$_dlaset( 'A', n, n, zero, one, z, ldz ) ! ==== quick return if possible ==== if( ilo==ihi ) then wr( ilo ) = h( ilo, ilo ) wi( ilo ) = zero return end if ! ==== stdlib${ii}$_dlahqr/stdlib${ii}$_dlaqr0 crossover point ==== nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'DHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork ) nmin = max( ntiny, nmin ) ! ==== stdlib${ii}$_dlaqr0 for big matrices; stdlib${ii}$_dlahqr for small ones ==== if( n>nmin ) then call stdlib${ii}$_dlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & work, lwork, info ) else ! ==== small matrix ==== call stdlib${ii}$_dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & info ) if( info>0_${ik}$ ) then ! ==== a rare stdlib${ii}$_dlahqr failure! stdlib${ii}$_dlaqr0 sometimes succeeds ! . when stdlib${ii}$_dlahqr fails. ==== kbot = info if( n>=nl ) then ! ==== larger matrices have enough subdiagonal scratch ! . space to call stdlib${ii}$_dlaqr0 directly. ==== call stdlib${ii}$_dlaqr0( wantt, wantz, n, ilo, kbot, h, ldh, wr,wi, ilo, ihi, z,& ldz, work, lwork, info ) else ! ==== tiny matrices don't have enough subdiagonal ! . scratch space to benefit from stdlib${ii}$_dlaqr0. hence, ! . tiny matrices must be copied into a larger ! . array before calling stdlib${ii}$_dlaqr0. ==== call stdlib${ii}$_dlacpy( 'A', n, n, h, ldh, hl, nl ) hl( n+1, n ) = zero call stdlib${ii}$_dlaset( 'A', nl, nl-n, zero, zero, hl( 1_${ik}$, n+1 ),nl ) call stdlib${ii}$_dlaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, wr,wi, ilo, ihi, & z, ldz, workl, nl, info ) if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_dlacpy( 'A', n, n, hl, nl, h, ldh ) end if end if end if ! ==== clear out the trash, if necessary. ==== if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_dlaset( 'L', n-2, n-2, zero, zero,& h( 3_${ik}$, 1_${ik}$ ), ldh ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=dp), work( 1_${ik}$ ) ) end if end subroutine stdlib${ii}$_dhseqr #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] module subroutine stdlib${ii}$_${ri}$hseqr( job, compz, n, ilo, ihi, h, ldh, wr, wi, z,ldz, work, lwork, info ) !! DHSEQR: computes the eigenvalues of a Hessenberg matrix H !! and, optionally, the matrices T and Z from the Schur decomposition !! H = Z T Z**T, where T is an upper quasi-triangular matrix (the !! Schur form), and Z is the orthogonal matrix of Schur vectors. !! Optionally Z may be postmultiplied into an input orthogonal !! matrix Q so that this routine can give the Schur factorization !! of a matrix A which has been reduced to the Hessenberg form H !! by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. ! -- lapack computational 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 integer(${ik}$), intent(in) :: ihi, ilo, ldh, ldz, lwork, n integer(${ik}$), intent(out) :: info character, intent(in) :: compz, job ! Array Arguments real(${rk}$), intent(inout) :: h(ldh,*), z(ldz,*) real(${rk}$), intent(out) :: wi(*), work(*), wr(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: ntiny = 15_${ik}$ integer(${ik}$), parameter :: nl = 49_${ik}$ ! ==== matrices of order ntiny or smaller must be processed by ! . stdlib${ii}$_${ri}$lahqr because of insufficient subdiagonal scratch space. ! . (this is a hard limit.) ==== ! ==== nl allocates some local workspace to help small matrices ! . through a rare stdlib${ii}$_${ri}$lahqr failure. nl > ntiny = 15 is ! . required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom- ! . mended. (the default value of nmin is 75.) using nl = 49 ! . allows up to six simultaneous shifts and a 16-by-16 ! . deflation window. ==== ! Local Arrays real(${rk}$) :: hl(nl,nl), workl(nl) ! Local Scalars integer(${ik}$) :: i, kbot, nmin logical(lk) :: initz, lquery, wantt, wantz ! Intrinsic Functions ! Executable Statements ! ==== decode and check the input parameters. ==== wantt = stdlib_lsame( job, 'S' ) initz = stdlib_lsame( compz, 'I' ) wantz = initz .or. stdlib_lsame( compz, 'V' ) work( 1_${ik}$ ) = real( max( 1_${ik}$, n ),KIND=${rk}$) lquery = lwork==-1_${ik}$ info = 0_${ik}$ if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then info = -1_${ik}$ else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( ihi<min( ilo, n ) .or. ihi>n ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then info = -11_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -13_${ik}$ end if if( info/=0_${ik}$ ) then ! ==== quick return in case of invalid argument. ==== call stdlib${ii}$_xerbla( 'DHSEQR', -info ) return else if( n==0_${ik}$ ) then ! ==== quick return in case n = 0; nothing to do. ==== return else if( lquery ) then ! ==== quick return in case of a workspace query ==== call stdlib${ii}$_${ri}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & work, lwork, info ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=${rk}$), work( 1_${ik}$ ) ) return else ! ==== copy eigenvalues isolated by stdlib${ii}$_${ri}$gebal ==== do i = 1, ilo - 1 wr( i ) = h( i, i ) wi( i ) = zero end do do i = ihi + 1, n wr( i ) = h( i, i ) wi( i ) = zero end do ! ==== initialize z, if requested ==== if( initz )call stdlib${ii}$_${ri}$laset( 'A', n, n, zero, one, z, ldz ) ! ==== quick return if possible ==== if( ilo==ihi ) then wr( ilo ) = h( ilo, ilo ) wi( ilo ) = zero return end if ! ==== stdlib${ii}$_${ri}$lahqr/stdlib${ii}$_${ri}$laqr0 crossover point ==== nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'DHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork ) nmin = max( ntiny, nmin ) ! ==== stdlib${ii}$_${ri}$laqr0 for big matrices; stdlib${ii}$_${ri}$lahqr for small ones ==== if( n>nmin ) then call stdlib${ii}$_${ri}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & work, lwork, info ) else ! ==== small matrix ==== call stdlib${ii}$_${ri}$lahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, ilo,ihi, z, ldz, & info ) if( info>0_${ik}$ ) then ! ==== a rare stdlib${ii}$_${ri}$lahqr failure! stdlib${ii}$_${ri}$laqr0 sometimes succeeds ! . when stdlib${ii}$_${ri}$lahqr fails. ==== kbot = info if( n>=nl ) then ! ==== larger matrices have enough subdiagonal scratch ! . space to call stdlib${ii}$_${ri}$laqr0 directly. ==== call stdlib${ii}$_${ri}$laqr0( wantt, wantz, n, ilo, kbot, h, ldh, wr,wi, ilo, ihi, z,& ldz, work, lwork, info ) else ! ==== tiny matrices don't have enough subdiagonal ! . scratch space to benefit from stdlib${ii}$_${ri}$laqr0. hence, ! . tiny matrices must be copied into a larger ! . array before calling stdlib${ii}$_${ri}$laqr0. ==== call stdlib${ii}$_${ri}$lacpy( 'A', n, n, h, ldh, hl, nl ) hl( n+1, n ) = zero call stdlib${ii}$_${ri}$laset( 'A', nl, nl-n, zero, zero, hl( 1_${ik}$, n+1 ),nl ) call stdlib${ii}$_${ri}$laqr0( wantt, wantz, nl, ilo, kbot, hl, nl, wr,wi, ilo, ihi, & z, ldz, workl, nl, info ) if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_${ri}$lacpy( 'A', n, n, hl, nl, h, ldh ) end if end if end if ! ==== clear out the trash, if necessary. ==== if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_${ri}$laset( 'L', n-2, n-2, zero, zero,& h( 3_${ik}$, 1_${ik}$ ), ldh ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = max( real( max( 1_${ik}$, n ),KIND=${rk}$), work( 1_${ik}$ ) ) end if end subroutine stdlib${ii}$_${ri}$hseqr #:endif #:endfor pure module subroutine stdlib${ii}$_chseqr( job, compz, n, ilo, ihi, h, ldh, w, z, ldz,work, lwork, info ) !! CHSEQR computes the eigenvalues of a Hessenberg matrix H !! and, optionally, the matrices T and Z from the Schur decomposition !! H = Z T Z**H, where T is an upper triangular matrix (the !! Schur form), and Z is the unitary matrix of Schur vectors. !! Optionally Z may be postmultiplied into an input unitary !! matrix Q so that this routine can give the Schur factorization !! of a matrix A which has been reduced to the Hessenberg form H !! by the unitary matrix Q: A = Q*H*Q**H = (QZ)*T*(QZ)**H. ! -- lapack computational 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 integer(${ik}$), intent(in) :: ihi, ilo, ldh, ldz, lwork, n integer(${ik}$), intent(out) :: info character, intent(in) :: compz, job ! Array Arguments complex(sp), intent(inout) :: h(ldh,*), z(ldz,*) complex(sp), intent(out) :: w(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: ntiny = 15_${ik}$ integer(${ik}$), parameter :: nl = 49_${ik}$ ! ==== matrices of order ntiny or smaller must be processed by ! . stdlib${ii}$_clahqr because of insufficient subdiagonal scratch space. ! . (this is a hard limit.) ==== ! ==== nl allocates some local workspace to help small matrices ! . through a rare stdlib${ii}$_clahqr failure. nl > ntiny = 15 is ! . required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom- ! . mended. (the default value of nmin is 75.) using nl = 49 ! . allows up to six simultaneous shifts and a 16-by-16 ! . deflation window. ==== ! Local Arrays complex(sp) :: hl(nl,nl), workl(nl) ! Local Scalars integer(${ik}$) :: kbot, nmin logical(lk) :: initz, lquery, wantt, wantz ! Intrinsic Functions ! Executable Statements ! ==== decode and check the input parameters. ==== wantt = stdlib_lsame( job, 'S' ) initz = stdlib_lsame( compz, 'I' ) wantz = initz .or. stdlib_lsame( compz, 'V' ) work( 1_${ik}$ ) = cmplx( real( max( 1_${ik}$, n ),KIND=sp), zero,KIND=sp) lquery = lwork==-1_${ik}$ info = 0_${ik}$ if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then info = -1_${ik}$ else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( ihi<min( ilo, n ) .or. ihi>n ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then info = -10_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -12_${ik}$ end if if( info/=0_${ik}$ ) then ! ==== quick return in case of invalid argument. ==== call stdlib${ii}$_xerbla( 'CHSEQR', -info ) return else if( n==0_${ik}$ ) then ! ==== quick return in case n = 0; nothing to do. ==== return else if( lquery ) then ! ==== quick return in case of a workspace query ==== call stdlib${ii}$_claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,ldz, work, & lwork, info ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = cmplx( max( real( work( 1_${ik}$ ),KIND=sp), real( max( 1_${ik}$,n ),KIND=sp) ), & zero,KIND=sp) return else ! ==== copy eigenvalues isolated by stdlib${ii}$_cgebal ==== if( ilo>1_${ik}$ )call stdlib${ii}$_ccopy( ilo-1, h, ldh+1, w, 1_${ik}$ ) if( ihi<n )call stdlib${ii}$_ccopy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1_${ik}$ ) ! ==== initialize z, if requested ==== if( initz )call stdlib${ii}$_claset( 'A', n, n, czero, cone, z, ldz ) ! ==== quick return if possible ==== if( ilo==ihi ) then w( ilo ) = h( ilo, ilo ) return end if ! ==== stdlib${ii}$_clahqr/stdlib${ii}$_claqr0 crossover point ==== nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'CHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork ) nmin = max( ntiny, nmin ) ! ==== stdlib${ii}$_claqr0 for big matrices; stdlib${ii}$_clahqr for small ones ==== if( n>nmin ) then call stdlib${ii}$_claqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, work, & lwork, info ) else ! ==== small matrix ==== call stdlib${ii}$_clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, info ) if( info>0_${ik}$ ) then ! ==== a rare stdlib${ii}$_clahqr failure! stdlib${ii}$_claqr0 sometimes succeeds ! . when stdlib${ii}$_clahqr fails. ==== kbot = info if( n>=nl ) then ! ==== larger matrices have enough subdiagonal scratch ! . space to call stdlib${ii}$_claqr0 directly. ==== call stdlib${ii}$_claqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,ilo, ihi, z, ldz,& work, lwork, info ) else ! ==== tiny matrices don't have enough subdiagonal ! . scratch space to benefit from stdlib${ii}$_claqr0. hence, ! . tiny matrices must be copied into a larger ! . array before calling stdlib${ii}$_claqr0. ==== call stdlib${ii}$_clacpy( 'A', n, n, h, ldh, hl, nl ) hl( n+1, n ) = czero call stdlib${ii}$_claset( 'A', nl, nl-n, czero, czero, hl( 1_${ik}$, n+1 ),nl ) call stdlib${ii}$_claqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,ilo, ihi, z, & ldz, workl, nl, info ) if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_clacpy( 'A', n, n, hl, nl, h, ldh ) end if end if end if ! ==== clear out the trash, if necessary. ==== if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_claset( 'L', n-2, n-2, czero, & czero, h( 3_${ik}$, 1_${ik}$ ), ldh ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = cmplx( max( real( max( 1_${ik}$, n ),KIND=sp),real( work( 1_${ik}$ ),KIND=sp) ), & zero,KIND=sp) end if end subroutine stdlib${ii}$_chseqr pure module subroutine stdlib${ii}$_zhseqr( job, compz, n, ilo, ihi, h, ldh, w, z, ldz,work, lwork, info ) !! ZHSEQR computes the eigenvalues of a Hessenberg matrix H !! and, optionally, the matrices T and Z from the Schur decomposition !! H = Z T Z**H, where T is an upper triangular matrix (the !! Schur form), and Z is the unitary matrix of Schur vectors. !! Optionally Z may be postmultiplied into an input unitary !! matrix Q so that this routine can give the Schur factorization !! of a matrix A which has been reduced to the Hessenberg form H !! by the unitary matrix Q: A = Q*H*Q**H = (QZ)*T*(QZ)**H. ! -- lapack computational 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 integer(${ik}$), intent(in) :: ihi, ilo, ldh, ldz, lwork, n integer(${ik}$), intent(out) :: info character, intent(in) :: compz, job ! Array Arguments complex(dp), intent(inout) :: h(ldh,*), z(ldz,*) complex(dp), intent(out) :: w(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: ntiny = 15_${ik}$ integer(${ik}$), parameter :: nl = 49_${ik}$ ! ==== matrices of order ntiny or smaller must be processed by ! . stdlib${ii}$_zlahqr because of insufficient subdiagonal scratch space. ! . (this is a hard limit.) ==== ! ==== nl allocates some local workspace to help small matrices ! . through a rare stdlib${ii}$_zlahqr failure. nl > ntiny = 15 is ! . required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom- ! . mended. (the default value of nmin is 75.) using nl = 49 ! . allows up to six simultaneous shifts and a 16-by-16 ! . deflation window. ==== ! Local Arrays complex(dp) :: hl(nl,nl), workl(nl) ! Local Scalars integer(${ik}$) :: kbot, nmin logical(lk) :: initz, lquery, wantt, wantz ! Intrinsic Functions ! Executable Statements ! ==== decode and check the input parameters. ==== wantt = stdlib_lsame( job, 'S' ) initz = stdlib_lsame( compz, 'I' ) wantz = initz .or. stdlib_lsame( compz, 'V' ) work( 1_${ik}$ ) = cmplx( real( max( 1_${ik}$, n ),KIND=dp), zero,KIND=dp) lquery = lwork==-1_${ik}$ info = 0_${ik}$ if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then info = -1_${ik}$ else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( ihi<min( ilo, n ) .or. ihi>n ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then info = -10_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -12_${ik}$ end if if( info/=0_${ik}$ ) then ! ==== quick return in case of invalid argument. ==== call stdlib${ii}$_xerbla( 'ZHSEQR', -info ) return else if( n==0_${ik}$ ) then ! ==== quick return in case n = 0; nothing to do. ==== return else if( lquery ) then ! ==== quick return in case of a workspace query ==== call stdlib${ii}$_zlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,ldz, work, & lwork, info ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = cmplx( max( real( work( 1_${ik}$ ),KIND=dp), real( max( 1_${ik}$,n ),KIND=dp) ), & zero,KIND=dp) return else ! ==== copy eigenvalues isolated by stdlib${ii}$_zgebal ==== if( ilo>1_${ik}$ )call stdlib${ii}$_zcopy( ilo-1, h, ldh+1, w, 1_${ik}$ ) if( ihi<n )call stdlib${ii}$_zcopy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1_${ik}$ ) ! ==== initialize z, if requested ==== if( initz )call stdlib${ii}$_zlaset( 'A', n, n, czero, cone, z, ldz ) ! ==== quick return if possible ==== if( ilo==ihi ) then w( ilo ) = h( ilo, ilo ) return end if ! ==== stdlib${ii}$_zlahqr/stdlib${ii}$_zlaqr0 crossover point ==== nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'ZHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork ) nmin = max( ntiny, nmin ) ! ==== stdlib${ii}$_zlaqr0 for big matrices; stdlib${ii}$_zlahqr for small ones ==== if( n>nmin ) then call stdlib${ii}$_zlaqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, work, & lwork, info ) else ! ==== small matrix ==== call stdlib${ii}$_zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, info ) if( info>0_${ik}$ ) then ! ==== a rare stdlib${ii}$_zlahqr failure! stdlib${ii}$_zlaqr0 sometimes succeeds ! . when stdlib${ii}$_zlahqr fails. ==== kbot = info if( n>=nl ) then ! ==== larger matrices have enough subdiagonal scratch ! . space to call stdlib${ii}$_zlaqr0 directly. ==== call stdlib${ii}$_zlaqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,ilo, ihi, z, ldz,& work, lwork, info ) else ! ==== tiny matrices don't have enough subdiagonal ! . scratch space to benefit from stdlib${ii}$_zlaqr0. hence, ! . tiny matrices must be copied into a larger ! . array before calling stdlib${ii}$_zlaqr0. ==== call stdlib${ii}$_zlacpy( 'A', n, n, h, ldh, hl, nl ) hl( n+1, n ) = czero call stdlib${ii}$_zlaset( 'A', nl, nl-n, czero, czero, hl( 1_${ik}$, n+1 ),nl ) call stdlib${ii}$_zlaqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,ilo, ihi, z, & ldz, workl, nl, info ) if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_zlacpy( 'A', n, n, hl, nl, h, ldh ) end if end if end if ! ==== clear out the trash, if necessary. ==== if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_zlaset( 'L', n-2, n-2, czero, & czero, h( 3_${ik}$, 1_${ik}$ ), ldh ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = cmplx( max( real( max( 1_${ik}$, n ),KIND=dp),real( work( 1_${ik}$ ),KIND=dp) ), & zero,KIND=dp) end if end subroutine stdlib${ii}$_zhseqr #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$hseqr( job, compz, n, ilo, ihi, h, ldh, w, z, ldz,work, lwork, info ) !! ZHSEQR: computes the eigenvalues of a Hessenberg matrix H !! and, optionally, the matrices T and Z from the Schur decomposition !! H = Z T Z**H, where T is an upper triangular matrix (the !! Schur form), and Z is the unitary matrix of Schur vectors. !! Optionally Z may be postmultiplied into an input unitary !! matrix Q so that this routine can give the Schur factorization !! of a matrix A which has been reduced to the Hessenberg form H !! by the unitary matrix Q: A = Q*H*Q**H = (QZ)*T*(QZ)**H. ! -- lapack computational 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 integer(${ik}$), intent(in) :: ihi, ilo, ldh, ldz, lwork, n integer(${ik}$), intent(out) :: info character, intent(in) :: compz, job ! Array Arguments complex(${ck}$), intent(inout) :: h(ldh,*), z(ldz,*) complex(${ck}$), intent(out) :: w(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: ntiny = 15_${ik}$ integer(${ik}$), parameter :: nl = 49_${ik}$ ! ==== matrices of order ntiny or smaller must be processed by ! . stdlib${ii}$_${ci}$lahqr because of insufficient subdiagonal scratch space. ! . (this is a hard limit.) ==== ! ==== nl allocates some local workspace to help small matrices ! . through a rare stdlib${ii}$_${ci}$lahqr failure. nl > ntiny = 15 is ! . required and nl <= nmin = stdlib${ii}$_ilaenv(ispec=12,...) is recom- ! . mended. (the default value of nmin is 75.) using nl = 49 ! . allows up to six simultaneous shifts and a 16-by-16 ! . deflation window. ==== ! Local Arrays complex(${ck}$) :: hl(nl,nl), workl(nl) ! Local Scalars integer(${ik}$) :: kbot, nmin logical(lk) :: initz, lquery, wantt, wantz ! Intrinsic Functions ! Executable Statements ! ==== decode and check the input parameters. ==== wantt = stdlib_lsame( job, 'S' ) initz = stdlib_lsame( compz, 'I' ) wantz = initz .or. stdlib_lsame( compz, 'V' ) work( 1_${ik}$ ) = cmplx( real( max( 1_${ik}$, n ),KIND=${ck}$), zero,KIND=${ck}$) lquery = lwork==-1_${ik}$ info = 0_${ik}$ if( .not.stdlib_lsame( job, 'E' ) .and. .not.wantt ) then info = -1_${ik}$ else if( .not.stdlib_lsame( compz, 'N' ) .and. .not.wantz ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ilo<1_${ik}$ .or. ilo>max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( ihi<min( ilo, n ) .or. ihi>n ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<max( 1_${ik}$, n ) ) ) then info = -10_${ik}$ else if( lwork<max( 1_${ik}$, n ) .and. .not.lquery ) then info = -12_${ik}$ end if if( info/=0_${ik}$ ) then ! ==== quick return in case of invalid argument. ==== call stdlib${ii}$_xerbla( 'ZHSEQR', -info ) return else if( n==0_${ik}$ ) then ! ==== quick return in case n = 0; nothing to do. ==== return else if( lquery ) then ! ==== quick return in case of a workspace query ==== call stdlib${ii}$_${ci}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi, z,ldz, work, & lwork, info ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = cmplx( max( real( work( 1_${ik}$ ),KIND=${ck}$), real( max( 1_${ik}$,n ),KIND=${ck}$) ), & zero,KIND=${ck}$) return else ! ==== copy eigenvalues isolated by stdlib${ii}$_${ci}$gebal ==== if( ilo>1_${ik}$ )call stdlib${ii}$_${ci}$copy( ilo-1, h, ldh+1, w, 1_${ik}$ ) if( ihi<n )call stdlib${ii}$_${ci}$copy( n-ihi, h( ihi+1, ihi+1 ), ldh+1, w( ihi+1 ), 1_${ik}$ ) ! ==== initialize z, if requested ==== if( initz )call stdlib${ii}$_${ci}$laset( 'A', n, n, czero, cone, z, ldz ) ! ==== quick return if possible ==== if( ilo==ihi ) then w( ilo ) = h( ilo, ilo ) return end if ! ==== stdlib${ii}$_${ci}$lahqr/stdlib${ii}$_${ci}$laqr0 crossover point ==== nmin = stdlib${ii}$_ilaenv( 12_${ik}$, 'ZHSEQR', job( : 1_${ik}$ ) // compz( : 1_${ik}$ ), n,ilo, ihi, lwork ) nmin = max( ntiny, nmin ) ! ==== stdlib${ii}$_${ci}$laqr0 for big matrices; stdlib${ii}$_${ci}$lahqr for small ones ==== if( n>nmin ) then call stdlib${ii}$_${ci}$laqr0( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, work, & lwork, info ) else ! ==== small matrix ==== call stdlib${ii}$_${ci}$lahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, ilo, ihi,z, ldz, info ) if( info>0_${ik}$ ) then ! ==== a rare stdlib${ii}$_${ci}$lahqr failure! stdlib${ii}$_${ci}$laqr0 sometimes succeeds ! . when stdlib${ii}$_${ci}$lahqr fails. ==== kbot = info if( n>=nl ) then ! ==== larger matrices have enough subdiagonal scratch ! . space to call stdlib${ii}$_${ci}$laqr0 directly. ==== call stdlib${ii}$_${ci}$laqr0( wantt, wantz, n, ilo, kbot, h, ldh, w,ilo, ihi, z, ldz,& work, lwork, info ) else ! ==== tiny matrices don't have enough subdiagonal ! . scratch space to benefit from stdlib${ii}$_${ci}$laqr0. hence, ! . tiny matrices must be copied into a larger ! . array before calling stdlib${ii}$_${ci}$laqr0. ==== call stdlib${ii}$_${ci}$lacpy( 'A', n, n, h, ldh, hl, nl ) hl( n+1, n ) = czero call stdlib${ii}$_${ci}$laset( 'A', nl, nl-n, czero, czero, hl( 1_${ik}$, n+1 ),nl ) call stdlib${ii}$_${ci}$laqr0( wantt, wantz, nl, ilo, kbot, hl, nl, w,ilo, ihi, z, & ldz, workl, nl, info ) if( wantt .or. info/=0_${ik}$ )call stdlib${ii}$_${ci}$lacpy( 'A', n, n, hl, nl, h, ldh ) end if end if end if ! ==== clear out the trash, if necessary. ==== if( ( wantt .or. info/=0_${ik}$ ) .and. n>2_${ik}$ )call stdlib${ii}$_${ci}$laset( 'L', n-2, n-2, czero, & czero, h( 3_${ik}$, 1_${ik}$ ), ldh ) ! ==== ensure reported workspace size is backward-compatible with ! . previous lapack versions. ==== work( 1_${ik}$ ) = cmplx( max( real( max( 1_${ik}$, n ),KIND=${ck}$),real( work( 1_${ik}$ ),KIND=${ck}$) ), & zero,KIND=${ck}$) end if end subroutine stdlib${ii}$_${ci}$hseqr #:endif #:endfor module subroutine stdlib${ii}$_shsein( side, eigsrc, initv, select, n, h, ldh, wr, wi,vl, ldvl, vr, ldvr, & !! SHSEIN uses inverse iteration to find specified right and/or left !! eigenvectors of a real upper Hessenberg matrix H. !! The right eigenvector x and the left eigenvector y of the matrix H !! corresponding to an eigenvalue w are defined by: !! H * x = w * x, y**h * H = w * y**h !! where y**h denotes the conjugate transpose of the vector y. mm, m, work, ifaill,ifailr, info ) ! -- lapack computational 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) :: eigsrc, initv, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(inout) :: select(*) integer(${ik}$), intent(out) :: ifaill(*), ifailr(*) real(sp), intent(in) :: h(ldh,*), wi(*) real(sp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), wr(*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: bothv, fromqr, leftv, noinit, pair, rightv integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ksi, ksr, ldwork real(sp) :: bignum, eps3, hnorm, smlnum, ulp, unfl, wki, wkr ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters. bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv fromqr = stdlib_lsame( eigsrc, 'Q' ) noinit = stdlib_lsame( initv, 'N' ) ! set m to the number of columns required to store the selected ! eigenvectors, and standardize the array select. m = 0_${ik}$ pair = .false. do k = 1, n if( pair ) then pair = .false. select( k ) = .false. else if( wi( k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) ) then select( k ) = .true. m = m + 2_${ik}$ end if end if end if end do info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then info = -2_${ik}$ else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then info = -3_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -11_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info = -13_${ik}$ else if( mm<m ) then info = -14_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SHSEIN', -info ) return end if ! quick return if possible. if( n==0 )return ! set machine-dependent constants. unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) ulp = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = unfl*( n / ulp ) bignum = ( one-ulp ) / smlnum ldwork = n + 1_${ik}$ kl = 1_${ik}$ kln = 0_${ik}$ if( fromqr ) then kr = 0_${ik}$ else kr = n end if ksr = 1_${ik}$ loop_120: do k = 1, n if( select( k ) ) then ! compute eigenvector(s) corresponding to w(k). if( fromqr ) then ! if affiliation of eigenvalues is known, check whether ! the matrix splits. ! determine kl and kr such that 1 <= kl <= k <= kr <= n ! and h(kl,kl-1) and h(kr+1,kr) are zero (or kl = 1 or ! kr = n). ! then inverse iteration can be performed with the ! submatrix h(kl:n,kl:n) for a left eigenvector, and with ! the submatrix h(1:kr,1:kr) for a right eigenvector. do i = k, kl + 1, -1 if( h( i, i-1 )==zero )go to 30 end do 30 continue kl = i if( k>kr ) then do i = k, n - 1 if( h( i+1, i )==zero )go to 50 end do 50 continue kr = i end if end if if( kl/=kln ) then kln = kl ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it ! has not ben computed before. hnorm = stdlib${ii}$_slanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work ) if( stdlib${ii}$_sisnan( hnorm ) ) then info = -6_${ik}$ return else if( hnorm>zero ) then eps3 = hnorm*ulp else eps3 = smlnum end if end if ! perturb eigenvalue if it is close to any previous ! selected eigenvalues affiliated to the submatrix ! h(kl:kr,kl:kr). close roots are modified by eps3. wkr = wr( k ) wki = wi( k ) 60 continue do i = k - 1, kl, -1 if( select( i ) .and. abs( wr( i )-wkr )+abs( wi( i )-wki )<eps3 ) & then wkr = wkr + eps3 go to 60 end if end do wr( k ) = wkr pair = wki/=zero if( pair ) then ksi = ksr + 1_${ik}$ else ksi = ksr end if if( leftv ) then ! compute left eigenvector. call stdlib${ii}$_slaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wkr, wki, vl( & kl, ksr ), vl( kl, ksi ),work, ldwork, work( n*n+n+1 ), eps3, smlnum,bignum, & iinfo ) if( iinfo>0_${ik}$ ) then if( pair ) then info = info + 2_${ik}$ else info = info + 1_${ik}$ end if ifaill( ksr ) = k ifaill( ksi ) = k else ifaill( ksr ) = 0_${ik}$ ifaill( ksi ) = 0_${ik}$ end if do i = 1, kl - 1 vl( i, ksr ) = zero end do if( pair ) then do i = 1, kl - 1 vl( i, ksi ) = zero end do end if end if if( rightv ) then ! compute right eigenvector. call stdlib${ii}$_slaein( .true., noinit, kr, h, ldh, wkr, wki,vr( 1_${ik}$, ksr ), vr( 1_${ik}$, & ksi ), work, ldwork,work( n*n+n+1 ), eps3, smlnum, bignum,iinfo ) if( iinfo>0_${ik}$ ) then if( pair ) then info = info + 2_${ik}$ else info = info + 1_${ik}$ end if ifailr( ksr ) = k ifailr( ksi ) = k else ifailr( ksr ) = 0_${ik}$ ifailr( ksi ) = 0_${ik}$ end if do i = kr + 1, n vr( i, ksr ) = zero end do if( pair ) then do i = kr + 1, n vr( i, ksi ) = zero end do end if end if if( pair ) then ksr = ksr + 2_${ik}$ else ksr = ksr + 1_${ik}$ end if end if end do loop_120 return end subroutine stdlib${ii}$_shsein module subroutine stdlib${ii}$_dhsein( side, eigsrc, initv, select, n, h, ldh, wr, wi,vl, ldvl, vr, ldvr, & !! DHSEIN uses inverse iteration to find specified right and/or left !! eigenvectors of a real upper Hessenberg matrix H. !! The right eigenvector x and the left eigenvector y of the matrix H !! corresponding to an eigenvalue w are defined by: !! H * x = w * x, y**h * H = w * y**h !! where y**h denotes the conjugate transpose of the vector y. mm, m, work, ifaill,ifailr, info ) ! -- lapack computational 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) :: eigsrc, initv, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(inout) :: select(*) integer(${ik}$), intent(out) :: ifaill(*), ifailr(*) real(dp), intent(in) :: h(ldh,*), wi(*) real(dp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), wr(*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: bothv, fromqr, leftv, noinit, pair, rightv integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ksi, ksr, ldwork real(dp) :: bignum, eps3, hnorm, smlnum, ulp, unfl, wki, wkr ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters. bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv fromqr = stdlib_lsame( eigsrc, 'Q' ) noinit = stdlib_lsame( initv, 'N' ) ! set m to the number of columns required to store the selected ! eigenvectors, and standardize the array select. m = 0_${ik}$ pair = .false. do k = 1, n if( pair ) then pair = .false. select( k ) = .false. else if( wi( k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) ) then select( k ) = .true. m = m + 2_${ik}$ end if end if end if end do info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then info = -2_${ik}$ else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then info = -3_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -11_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info = -13_${ik}$ else if( mm<m ) then info = -14_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DHSEIN', -info ) return end if ! quick return if possible. if( n==0 )return ! set machine-dependent constants. unfl = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) ulp = stdlib${ii}$_dlamch( 'PRECISION' ) smlnum = unfl*( n / ulp ) bignum = ( one-ulp ) / smlnum ldwork = n + 1_${ik}$ kl = 1_${ik}$ kln = 0_${ik}$ if( fromqr ) then kr = 0_${ik}$ else kr = n end if ksr = 1_${ik}$ loop_120: do k = 1, n if( select( k ) ) then ! compute eigenvector(s) corresponding to w(k). if( fromqr ) then ! if affiliation of eigenvalues is known, check whether ! the matrix splits. ! determine kl and kr such that 1 <= kl <= k <= kr <= n ! and h(kl,kl-1) and h(kr+1,kr) are zero (or kl = 1 or ! kr = n). ! then inverse iteration can be performed with the ! submatrix h(kl:n,kl:n) for a left eigenvector, and with ! the submatrix h(1:kr,1:kr) for a right eigenvector. do i = k, kl + 1, -1 if( h( i, i-1 )==zero )go to 30 end do 30 continue kl = i if( k>kr ) then do i = k, n - 1 if( h( i+1, i )==zero )go to 50 end do 50 continue kr = i end if end if if( kl/=kln ) then kln = kl ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it ! has not ben computed before. hnorm = stdlib${ii}$_dlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work ) if( stdlib${ii}$_disnan( hnorm ) ) then info = -6_${ik}$ return else if( hnorm>zero ) then eps3 = hnorm*ulp else eps3 = smlnum end if end if ! perturb eigenvalue if it is close to any previous ! selected eigenvalues affiliated to the submatrix ! h(kl:kr,kl:kr). close roots are modified by eps3. wkr = wr( k ) wki = wi( k ) 60 continue do i = k - 1, kl, -1 if( select( i ) .and. abs( wr( i )-wkr )+abs( wi( i )-wki )<eps3 ) & then wkr = wkr + eps3 go to 60 end if end do wr( k ) = wkr pair = wki/=zero if( pair ) then ksi = ksr + 1_${ik}$ else ksi = ksr end if if( leftv ) then ! compute left eigenvector. call stdlib${ii}$_dlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wkr, wki, vl( & kl, ksr ), vl( kl, ksi ),work, ldwork, work( n*n+n+1 ), eps3, smlnum,bignum, & iinfo ) if( iinfo>0_${ik}$ ) then if( pair ) then info = info + 2_${ik}$ else info = info + 1_${ik}$ end if ifaill( ksr ) = k ifaill( ksi ) = k else ifaill( ksr ) = 0_${ik}$ ifaill( ksi ) = 0_${ik}$ end if do i = 1, kl - 1 vl( i, ksr ) = zero end do if( pair ) then do i = 1, kl - 1 vl( i, ksi ) = zero end do end if end if if( rightv ) then ! compute right eigenvector. call stdlib${ii}$_dlaein( .true., noinit, kr, h, ldh, wkr, wki,vr( 1_${ik}$, ksr ), vr( 1_${ik}$, & ksi ), work, ldwork,work( n*n+n+1 ), eps3, smlnum, bignum,iinfo ) if( iinfo>0_${ik}$ ) then if( pair ) then info = info + 2_${ik}$ else info = info + 1_${ik}$ end if ifailr( ksr ) = k ifailr( ksi ) = k else ifailr( ksr ) = 0_${ik}$ ifailr( ksi ) = 0_${ik}$ end if do i = kr + 1, n vr( i, ksr ) = zero end do if( pair ) then do i = kr + 1, n vr( i, ksi ) = zero end do end if end if if( pair ) then ksr = ksr + 2_${ik}$ else ksr = ksr + 1_${ik}$ end if end if end do loop_120 return end subroutine stdlib${ii}$_dhsein #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] module subroutine stdlib${ii}$_${ri}$hsein( side, eigsrc, initv, select, n, h, ldh, wr, wi,vl, ldvl, vr, ldvr, & !! DHSEIN: uses inverse iteration to find specified right and/or left !! eigenvectors of a real upper Hessenberg matrix H. !! The right eigenvector x and the left eigenvector y of the matrix H !! corresponding to an eigenvalue w are defined by: !! H * x = w * x, y**h * H = w * y**h !! where y**h denotes the conjugate transpose of the vector y. mm, m, work, ifaill,ifailr, info ) ! -- lapack computational 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) :: eigsrc, initv, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(inout) :: select(*) integer(${ik}$), intent(out) :: ifaill(*), ifailr(*) real(${rk}$), intent(in) :: h(ldh,*), wi(*) real(${rk}$), intent(inout) :: vl(ldvl,*), vr(ldvr,*), wr(*) real(${rk}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: bothv, fromqr, leftv, noinit, pair, rightv integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ksi, ksr, ldwork real(${rk}$) :: bignum, eps3, hnorm, smlnum, ulp, unfl, wki, wkr ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters. bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv fromqr = stdlib_lsame( eigsrc, 'Q' ) noinit = stdlib_lsame( initv, 'N' ) ! set m to the number of columns required to store the selected ! eigenvectors, and standardize the array select. m = 0_${ik}$ pair = .false. do k = 1, n if( pair ) then pair = .false. select( k ) = .false. else if( wi( k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) ) then select( k ) = .true. m = m + 2_${ik}$ end if end if end if end do info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then info = -2_${ik}$ else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then info = -3_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -11_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info = -13_${ik}$ else if( mm<m ) then info = -14_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DHSEIN', -info ) return end if ! quick return if possible. if( n==0 )return ! set machine-dependent constants. unfl = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) ulp = stdlib${ii}$_${ri}$lamch( 'PRECISION' ) smlnum = unfl*( n / ulp ) bignum = ( one-ulp ) / smlnum ldwork = n + 1_${ik}$ kl = 1_${ik}$ kln = 0_${ik}$ if( fromqr ) then kr = 0_${ik}$ else kr = n end if ksr = 1_${ik}$ loop_120: do k = 1, n if( select( k ) ) then ! compute eigenvector(s) corresponding to w(k). if( fromqr ) then ! if affiliation of eigenvalues is known, check whether ! the matrix splits. ! determine kl and kr such that 1 <= kl <= k <= kr <= n ! and h(kl,kl-1) and h(kr+1,kr) are zero (or kl = 1 or ! kr = n). ! then inverse iteration can be performed with the ! submatrix h(kl:n,kl:n) for a left eigenvector, and with ! the submatrix h(1:kr,1:kr) for a right eigenvector. do i = k, kl + 1, -1 if( h( i, i-1 )==zero )go to 30 end do 30 continue kl = i if( k>kr ) then do i = k, n - 1 if( h( i+1, i )==zero )go to 50 end do 50 continue kr = i end if end if if( kl/=kln ) then kln = kl ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it ! has not ben computed before. hnorm = stdlib${ii}$_${ri}$lanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work ) if( stdlib${ii}$_${ri}$isnan( hnorm ) ) then info = -6_${ik}$ return else if( hnorm>zero ) then eps3 = hnorm*ulp else eps3 = smlnum end if end if ! perturb eigenvalue if it is close to any previous ! selected eigenvalues affiliated to the submatrix ! h(kl:kr,kl:kr). close roots are modified by eps3. wkr = wr( k ) wki = wi( k ) 60 continue do i = k - 1, kl, -1 if( select( i ) .and. abs( wr( i )-wkr )+abs( wi( i )-wki )<eps3 ) & then wkr = wkr + eps3 go to 60 end if end do wr( k ) = wkr pair = wki/=zero if( pair ) then ksi = ksr + 1_${ik}$ else ksi = ksr end if if( leftv ) then ! compute left eigenvector. call stdlib${ii}$_${ri}$laein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wkr, wki, vl( & kl, ksr ), vl( kl, ksi ),work, ldwork, work( n*n+n+1 ), eps3, smlnum,bignum, & iinfo ) if( iinfo>0_${ik}$ ) then if( pair ) then info = info + 2_${ik}$ else info = info + 1_${ik}$ end if ifaill( ksr ) = k ifaill( ksi ) = k else ifaill( ksr ) = 0_${ik}$ ifaill( ksi ) = 0_${ik}$ end if do i = 1, kl - 1 vl( i, ksr ) = zero end do if( pair ) then do i = 1, kl - 1 vl( i, ksi ) = zero end do end if end if if( rightv ) then ! compute right eigenvector. call stdlib${ii}$_${ri}$laein( .true., noinit, kr, h, ldh, wkr, wki,vr( 1_${ik}$, ksr ), vr( 1_${ik}$, & ksi ), work, ldwork,work( n*n+n+1 ), eps3, smlnum, bignum,iinfo ) if( iinfo>0_${ik}$ ) then if( pair ) then info = info + 2_${ik}$ else info = info + 1_${ik}$ end if ifailr( ksr ) = k ifailr( ksi ) = k else ifailr( ksr ) = 0_${ik}$ ifailr( ksi ) = 0_${ik}$ end if do i = kr + 1, n vr( i, ksr ) = zero end do if( pair ) then do i = kr + 1, n vr( i, ksi ) = zero end do end if end if if( pair ) then ksr = ksr + 2_${ik}$ else ksr = ksr + 1_${ik}$ end if end if end do loop_120 return end subroutine stdlib${ii}$_${ri}$hsein #:endif #:endfor module subroutine stdlib${ii}$_chsein( side, eigsrc, initv, select, n, h, ldh, w, vl,ldvl, vr, ldvr, mm, & !! CHSEIN uses inverse iteration to find specified right and/or left !! eigenvectors of a complex upper Hessenberg matrix H. !! The right eigenvector x and the left eigenvector y of the matrix H !! corresponding to an eigenvalue w are defined by: !! H * x = w * x, y**h * H = w * y**h !! where y**h denotes the conjugate transpose of the vector y. m, work, rwork, ifaill,ifailr, info ) ! -- lapack computational 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) :: eigsrc, initv, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: ifaill(*), ifailr(*) real(sp), intent(out) :: rwork(*) complex(sp), intent(in) :: h(ldh,*) complex(sp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), w(*) complex(sp), intent(out) :: work(*) ! ===================================================================== ! Parameters ! Local Scalars logical(lk) :: bothv, fromqr, leftv, noinit, rightv integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ks, ldwork real(sp) :: eps3, hnorm, smlnum, ulp, unfl complex(sp) :: cdum, wk ! Intrinsic Functions ! Statement Functions real(sp) :: cabs1 ! Statement Function Definitions cabs1( cdum ) = abs( real( cdum,KIND=sp) ) + abs( aimag( cdum ) ) ! Executable Statements ! decode and test the input parameters. bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv fromqr = stdlib_lsame( eigsrc, 'Q' ) noinit = stdlib_lsame( initv, 'N' ) ! set m to the number of columns required to store the selected ! eigenvectors. m = 0_${ik}$ do k = 1, n if( select( k ) )m = m + 1_${ik}$ end do info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then info = -2_${ik}$ else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then info = -3_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -10_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info = -12_${ik}$ else if( mm<m ) then info = -13_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CHSEIN', -info ) return end if ! quick return if possible. if( n==0 )return ! set machine-dependent constants. unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) ulp = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = unfl*( n / ulp ) ldwork = n kl = 1_${ik}$ kln = 0_${ik}$ if( fromqr ) then kr = 0_${ik}$ else kr = n end if ks = 1_${ik}$ loop_100: do k = 1, n if( select( k ) ) then ! compute eigenvector(s) corresponding to w(k). if( fromqr ) then ! if affiliation of eigenvalues is known, check whether ! the matrix splits. ! determine kl and kr such that 1 <= kl <= k <= kr <= n ! and h(kl,kl-1) and h(kr+1,kr) are czero (or kl = 1 or ! kr = n). ! then inverse iteration can be performed with the ! submatrix h(kl:n,kl:n) for a left eigenvector, and with ! the submatrix h(1:kr,1:kr) for a right eigenvector. do i = k, kl + 1, -1 if( h( i, i-1 )==czero )go to 30 end do 30 continue kl = i if( k>kr ) then do i = k, n - 1 if( h( i+1, i )==czero )go to 50 end do 50 continue kr = i end if end if if( kl/=kln ) then kln = kl ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it ! has not ben computed before. hnorm = stdlib${ii}$_clanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork ) if( stdlib${ii}$_sisnan( hnorm ) ) then info = -6_${ik}$ return else if( (hnorm>zero) ) then eps3 = hnorm*ulp else eps3 = smlnum end if end if ! perturb eigenvalue if it is close to any previous ! selected eigenvalues affiliated to the submatrix ! h(kl:kr,kl:kr). close roots are modified by eps3. wk = w( k ) 60 continue do i = k - 1, kl, -1 if( select( i ) .and. cabs1( w( i )-wk )<eps3 ) then wk = wk + eps3 go to 60 end if end do w( k ) = wk if( leftv ) then ! compute left eigenvector. call stdlib${ii}$_claein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wk, vl( kl, ks )& , work, ldwork, rwork, eps3,smlnum, iinfo ) if( iinfo>0_${ik}$ ) then info = info + 1_${ik}$ ifaill( ks ) = k else ifaill( ks ) = 0_${ik}$ end if do i = 1, kl - 1 vl( i, ks ) = czero end do end if if( rightv ) then ! compute right eigenvector. call stdlib${ii}$_claein( .true., noinit, kr, h, ldh, wk, vr( 1_${ik}$, ks ),work, ldwork, & rwork, eps3, smlnum, iinfo ) if( iinfo>0_${ik}$ ) then info = info + 1_${ik}$ ifailr( ks ) = k else ifailr( ks ) = 0_${ik}$ end if do i = kr + 1, n vr( i, ks ) = czero end do end if ks = ks + 1_${ik}$ end if end do loop_100 return end subroutine stdlib${ii}$_chsein module subroutine stdlib${ii}$_zhsein( side, eigsrc, initv, select, n, h, ldh, w, vl,ldvl, vr, ldvr, mm, & !! ZHSEIN uses inverse iteration to find specified right and/or left !! eigenvectors of a complex upper Hessenberg matrix H. !! The right eigenvector x and the left eigenvector y of the matrix H !! corresponding to an eigenvalue w are defined by: !! H * x = w * x, y**h * H = w * y**h !! where y**h denotes the conjugate transpose of the vector y. m, work, rwork, ifaill,ifailr, info ) ! -- lapack computational 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) :: eigsrc, initv, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: ifaill(*), ifailr(*) real(dp), intent(out) :: rwork(*) complex(dp), intent(in) :: h(ldh,*) complex(dp), intent(inout) :: vl(ldvl,*), vr(ldvr,*), w(*) complex(dp), intent(out) :: work(*) ! ===================================================================== ! Parameters ! Local Scalars logical(lk) :: bothv, fromqr, leftv, noinit, rightv integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ks, ldwork real(dp) :: eps3, hnorm, smlnum, ulp, unfl complex(dp) :: cdum, wk ! Intrinsic Functions ! Statement Functions real(dp) :: cabs1 ! Statement Function Definitions cabs1( cdum ) = abs( real( cdum,KIND=dp) ) + abs( aimag( cdum ) ) ! Executable Statements ! decode and test the input parameters. bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv fromqr = stdlib_lsame( eigsrc, 'Q' ) noinit = stdlib_lsame( initv, 'N' ) ! set m to the number of columns required to store the selected ! eigenvectors. m = 0_${ik}$ do k = 1, n if( select( k ) )m = m + 1_${ik}$ end do info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then info = -2_${ik}$ else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then info = -3_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -10_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info = -12_${ik}$ else if( mm<m ) then info = -13_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZHSEIN', -info ) return end if ! quick return if possible. if( n==0 )return ! set machine-dependent constants. unfl = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) ulp = stdlib${ii}$_dlamch( 'PRECISION' ) smlnum = unfl*( n / ulp ) ldwork = n kl = 1_${ik}$ kln = 0_${ik}$ if( fromqr ) then kr = 0_${ik}$ else kr = n end if ks = 1_${ik}$ loop_100: do k = 1, n if( select( k ) ) then ! compute eigenvector(s) corresponding to w(k). if( fromqr ) then ! if affiliation of eigenvalues is known, check whether ! the matrix splits. ! determine kl and kr such that 1 <= kl <= k <= kr <= n ! and h(kl,kl-1) and h(kr+1,kr) are czero (or kl = 1 or ! kr = n). ! then inverse iteration can be performed with the ! submatrix h(kl:n,kl:n) for a left eigenvector, and with ! the submatrix h(1:kr,1:kr) for a right eigenvector. do i = k, kl + 1, -1 if( h( i, i-1 )==czero )go to 30 end do 30 continue kl = i if( k>kr ) then do i = k, n - 1 if( h( i+1, i )==czero )go to 50 end do 50 continue kr = i end if end if if( kl/=kln ) then kln = kl ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it ! has not ben computed before. hnorm = stdlib${ii}$_zlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork ) if( stdlib${ii}$_disnan( hnorm ) ) then info = -6_${ik}$ return else if( hnorm>zero ) then eps3 = hnorm*ulp else eps3 = smlnum end if end if ! perturb eigenvalue if it is close to any previous ! selected eigenvalues affiliated to the submatrix ! h(kl:kr,kl:kr). close roots are modified by eps3. wk = w( k ) 60 continue do i = k - 1, kl, -1 if( select( i ) .and. cabs1( w( i )-wk )<eps3 ) then wk = wk + eps3 go to 60 end if end do w( k ) = wk if( leftv ) then ! compute left eigenvector. call stdlib${ii}$_zlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wk, vl( kl, ks )& , work, ldwork, rwork, eps3,smlnum, iinfo ) if( iinfo>0_${ik}$ ) then info = info + 1_${ik}$ ifaill( ks ) = k else ifaill( ks ) = 0_${ik}$ end if do i = 1, kl - 1 vl( i, ks ) = czero end do end if if( rightv ) then ! compute right eigenvector. call stdlib${ii}$_zlaein( .true., noinit, kr, h, ldh, wk, vr( 1_${ik}$, ks ),work, ldwork, & rwork, eps3, smlnum, iinfo ) if( iinfo>0_${ik}$ ) then info = info + 1_${ik}$ ifailr( ks ) = k else ifailr( ks ) = 0_${ik}$ end if do i = kr + 1, n vr( i, ks ) = czero end do end if ks = ks + 1_${ik}$ end if end do loop_100 return end subroutine stdlib${ii}$_zhsein #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] module subroutine stdlib${ii}$_${ci}$hsein( side, eigsrc, initv, select, n, h, ldh, w, vl,ldvl, vr, ldvr, mm, & !! ZHSEIN: uses inverse iteration to find specified right and/or left !! eigenvectors of a complex upper Hessenberg matrix H. !! The right eigenvector x and the left eigenvector y of the matrix H !! corresponding to an eigenvalue w are defined by: !! H * x = w * x, y**h * H = w * y**h !! where y**h denotes the conjugate transpose of the vector y. m, work, rwork, ifaill,ifailr, info ) ! -- lapack computational 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) :: eigsrc, initv, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldh, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: ifaill(*), ifailr(*) real(${ck}$), intent(out) :: rwork(*) complex(${ck}$), intent(in) :: h(ldh,*) complex(${ck}$), intent(inout) :: vl(ldvl,*), vr(ldvr,*), w(*) complex(${ck}$), intent(out) :: work(*) ! ===================================================================== ! Parameters ! Local Scalars logical(lk) :: bothv, fromqr, leftv, noinit, rightv integer(${ik}$) :: i, iinfo, k, kl, kln, kr, ks, ldwork real(${ck}$) :: eps3, hnorm, smlnum, ulp, unfl complex(${ck}$) :: cdum, wk ! Intrinsic Functions ! Statement Functions real(${ck}$) :: cabs1 ! Statement Function Definitions cabs1( cdum ) = abs( real( cdum,KIND=${ck}$) ) + abs( aimag( cdum ) ) ! Executable Statements ! decode and test the input parameters. bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv fromqr = stdlib_lsame( eigsrc, 'Q' ) noinit = stdlib_lsame( initv, 'N' ) ! set m to the number of columns required to store the selected ! eigenvectors. m = 0_${ik}$ do k = 1, n if( select( k ) )m = m + 1_${ik}$ end do info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.fromqr .and. .not.stdlib_lsame( eigsrc, 'N' ) ) then info = -2_${ik}$ else if( .not.noinit .and. .not.stdlib_lsame( initv, 'U' ) ) then info = -3_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( ldh<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -10_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info = -12_${ik}$ else if( mm<m ) then info = -13_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZHSEIN', -info ) return end if ! quick return if possible. if( n==0 )return ! set machine-dependent constants. unfl = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' ) ulp = stdlib${ii}$_${c2ri(ci)}$lamch( 'PRECISION' ) smlnum = unfl*( n / ulp ) ldwork = n kl = 1_${ik}$ kln = 0_${ik}$ if( fromqr ) then kr = 0_${ik}$ else kr = n end if ks = 1_${ik}$ loop_100: do k = 1, n if( select( k ) ) then ! compute eigenvector(s) corresponding to w(k). if( fromqr ) then ! if affiliation of eigenvalues is known, check whether ! the matrix splits. ! determine kl and kr such that 1 <= kl <= k <= kr <= n ! and h(kl,kl-1) and h(kr+1,kr) are czero (or kl = 1 or ! kr = n). ! then inverse iteration can be performed with the ! submatrix h(kl:n,kl:n) for a left eigenvector, and with ! the submatrix h(1:kr,1:kr) for a right eigenvector. do i = k, kl + 1, -1 if( h( i, i-1 )==czero )go to 30 end do 30 continue kl = i if( k>kr ) then do i = k, n - 1 if( h( i+1, i )==czero )go to 50 end do 50 continue kr = i end if end if if( kl/=kln ) then kln = kl ! compute infinity-norm of submatrix h(kl:kr,kl:kr) if it ! has not ben computed before. hnorm = stdlib${ii}$_${ci}$lanhs( 'I', kr-kl+1, h( kl, kl ), ldh, rwork ) if( stdlib${ii}$_${c2ri(ci)}$isnan( hnorm ) ) then info = -6_${ik}$ return else if( hnorm>zero ) then eps3 = hnorm*ulp else eps3 = smlnum end if end if ! perturb eigenvalue if it is close to any previous ! selected eigenvalues affiliated to the submatrix ! h(kl:kr,kl:kr). close roots are modified by eps3. wk = w( k ) 60 continue do i = k - 1, kl, -1 if( select( i ) .and. cabs1( w( i )-wk )<eps3 ) then wk = wk + eps3 go to 60 end if end do w( k ) = wk if( leftv ) then ! compute left eigenvector. call stdlib${ii}$_${ci}$laein( .false., noinit, n-kl+1, h( kl, kl ), ldh,wk, vl( kl, ks )& , work, ldwork, rwork, eps3,smlnum, iinfo ) if( iinfo>0_${ik}$ ) then info = info + 1_${ik}$ ifaill( ks ) = k else ifaill( ks ) = 0_${ik}$ end if do i = 1, kl - 1 vl( i, ks ) = czero end do end if if( rightv ) then ! compute right eigenvector. call stdlib${ii}$_${ci}$laein( .true., noinit, kr, h, ldh, wk, vr( 1_${ik}$, ks ),work, ldwork, & rwork, eps3, smlnum, iinfo ) if( iinfo>0_${ik}$ ) then info = info + 1_${ik}$ ifailr( ks ) = k else ifailr( ks ) = 0_${ik}$ end if do i = kr + 1, n vr( i, ks ) = czero end do end if ks = ks + 1_${ik}$ end if end do loop_100 return end subroutine stdlib${ii}$_${ci}$hsein #:endif #:endfor pure module subroutine stdlib${ii}$_strevc( side, howmny, select, n, t, ldt, vl, ldvl, vr,ldvr, mm, m, & !! STREVC computes some or all of the right and/or left eigenvectors of !! a real upper quasi-triangular matrix T. !! Matrices of this type are produced by the Schur factorization of !! a real general matrix: A = Q*T*Q**T, as computed by SHSEQR. !! The right eigenvector x and the left eigenvector y of T corresponding !! to an eigenvalue w are defined by: !! T*x = w*x, (y**H)*T = w*(y**H) !! where y**H denotes the conjugate transpose of y. !! The eigenvalues are not input to this routine, but are read directly !! from the diagonal blocks of T. !! This routine returns the matrices X and/or Y of right and left !! eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an !! input matrix. If Q is the orthogonal factor that reduces a matrix !! A to Schur form T, then Q*X and Q*Y are the matrices of right and !! left eigenvectors of A. work, info ) ! -- lapack computational 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) :: howmny, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldt, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(inout) :: select(*) real(sp), intent(in) :: t(ldt,*) real(sp), intent(inout) :: vl(ldvl,*), vr(ldvr,*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: allv, bothv, leftv, over, pair, rightv, somev integer(${ik}$) :: i, ierr, ii, ip, is, j, j1, j2, jnxt, k, ki, n2 real(sp) :: beta, bignum, emax, ovfl, rec, remax, scale, smin, smlnum, ulp, unfl, & vcrit, vmax, wi, wr, xnorm ! Intrinsic Functions ! Local Arrays real(sp) :: x(2_${ik}$,2_${ik}$) ! Executable Statements ! decode and test the input parameters bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv allv = stdlib_lsame( howmny, 'A' ) over = stdlib_lsame( howmny, 'B' ) somev = stdlib_lsame( howmny, 'S' ) info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.allv .and. .not.over .and. .not.somev ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -4_${ik}$ else if( ldt<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -8_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info = -10_${ik}$ else ! set m to the number of columns required to store the selected ! eigenvectors, standardize the array select if necessary, and ! test mm. if( somev ) then m = 0_${ik}$ pair = .false. do j = 1, n if( pair ) then pair = .false. select( j ) = .false. else if( j<n ) then if( t( j+1, j )==zero ) then if( select( j ) )m = m + 1_${ik}$ else pair = .true. if( select( j ) .or. select( j+1 ) ) then select( j ) = .true. m = m + 2_${ik}$ end if end if else if( select( n ) )m = m + 1_${ik}$ end if end if end do else m = n end if if( mm<m ) then info = -11_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'STREVC', -info ) return end if ! quick return if possible. if( n==0 )return ! set the constants to control overflow. unfl = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) ovfl = one / unfl call stdlib${ii}$_slabad( unfl, ovfl ) ulp = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = unfl*( n / ulp ) bignum = ( one-ulp ) / smlnum ! compute 1-norm of each column of strictly upper triangular ! part of t to control overflow in triangular solver. work( 1_${ik}$ ) = zero do j = 2, n work( j ) = zero do i = 1, j - 1 work( j ) = work( j ) + abs( t( i, j ) ) end do end do ! index ip is used to specify the real or complex eigenvalue: ! ip = 0, real eigenvalue, ! 1, first of conjugate complex pair: (wr,wi) ! -1, second of conjugate complex pair: (wr,wi) n2 = 2_${ik}$*n if( rightv ) then ! compute right eigenvectors. ip = 0_${ik}$ is = m loop_140: do ki = n, 1, -1 if( ip==1 )go to 130 if( ki==1 )go to 40 if( t( ki, ki-1 )==zero )go to 40 ip = -1_${ik}$ 40 continue if( somev ) then if( ip==0_${ik}$ ) then if( .not.select( ki ) )go to 130 else if( .not.select( ki-1 ) )go to 130 end if end if ! compute the ki-th eigenvalue (wr,wi). wr = t( ki, ki ) wi = zero if( ip/=0_${ik}$ )wi = sqrt( abs( t( ki, ki-1 ) ) )*sqrt( abs( t( ki-1, ki ) ) ) smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum ) if( ip==0_${ik}$ ) then ! real right eigenvector work( ki+n ) = one ! form right-hand side do k = 1, ki - 1 work( k+n ) = -t( k, ki ) end do ! solve the upper quasi-triangular system: ! (t(1:ki-1,1:ki-1) - wr)*x = scale*work. jnxt = ki - 1_${ik}$ loop_60: do j = ki - 1, 1, -1 if( j>jnxt )cycle loop_60 j1 = j j2 = j jnxt = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnxt = j - 2_${ik}$ end if end if if( j1==j2 ) then ! 1-by-1 diagonal block call stdlib${ii}$_slaln2( .false., 1_${ik}$, 1_${ik}$, smin, one, t( j, j ),ldt, one, one, & work( j+n ), n, wr,zero, x, 2_${ik}$, scale, xnorm, ierr ) ! scale x(1,1) to avoid overflow when updating ! the right-hand side. if( xnorm>one ) then if( work( j )>bignum / xnorm ) then x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ ) / xnorm scale = scale / xnorm end if end if ! scale if necessary if( scale/=one )call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ ) work( j+n ) = x( 1_${ik}$, 1_${ik}$ ) ! update right-hand side call stdlib${ii}$_saxpy( j-1, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ ) else ! 2-by-2 diagonal block call stdlib${ii}$_slaln2( .false., 2_${ik}$, 1_${ik}$, smin, one,t( j-1, j-1 ), ldt, one, & one,work( j-1+n ), n, wr, zero, x, 2_${ik}$,scale, xnorm, ierr ) ! scale x(1,1) and x(2,1) to avoid overflow when ! updating the right-hand side. if( xnorm>one ) then beta = max( work( j-1 ), work( j ) ) if( beta>bignum / xnorm ) then x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ ) / xnorm x( 2_${ik}$, 1_${ik}$ ) = x( 2_${ik}$, 1_${ik}$ ) / xnorm scale = scale / xnorm end if end if ! scale if necessary if( scale/=one )call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ ) work( j-1+n ) = x( 1_${ik}$, 1_${ik}$ ) work( j+n ) = x( 2_${ik}$, 1_${ik}$ ) ! update right-hand side call stdlib${ii}$_saxpy( j-2, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j-1 ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ ) call stdlib${ii}$_saxpy( j-2, -x( 2_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ ) end if end do loop_60 ! copy the vector x or q*x to vr and normalize. if( .not.over ) then call stdlib${ii}$_scopy( ki, work( 1_${ik}$+n ), 1_${ik}$, vr( 1_${ik}$, is ), 1_${ik}$ ) ii = stdlib${ii}$_isamax( ki, vr( 1_${ik}$, is ), 1_${ik}$ ) remax = one / abs( vr( ii, is ) ) call stdlib${ii}$_sscal( ki, remax, vr( 1_${ik}$, is ), 1_${ik}$ ) do k = ki + 1, n vr( k, is ) = zero end do else if( ki>1_${ik}$ )call stdlib${ii}$_sgemv( 'N', n, ki-1, one, vr, ldvr,work( 1_${ik}$+n ), 1_${ik}$, & work( ki+n ),vr( 1_${ik}$, ki ), 1_${ik}$ ) ii = stdlib${ii}$_isamax( n, vr( 1_${ik}$, ki ), 1_${ik}$ ) remax = one / abs( vr( ii, ki ) ) call stdlib${ii}$_sscal( n, remax, vr( 1_${ik}$, ki ), 1_${ik}$ ) end if else ! complex right eigenvector. ! initial solve ! [ (t(ki-1,ki-1) t(ki-1,ki) ) - (wr + i* wi)]*x = 0. ! [ (t(ki,ki-1) t(ki,ki) ) ] if( abs( t( ki-1, ki ) )>=abs( t( ki, ki-1 ) ) ) then work( ki-1+n ) = one work( ki+n2 ) = wi / t( ki-1, ki ) else work( ki-1+n ) = -wi / t( ki, ki-1 ) work( ki+n2 ) = one end if work( ki+n ) = zero work( ki-1+n2 ) = zero ! form right-hand side do k = 1, ki - 2 work( k+n ) = -work( ki-1+n )*t( k, ki-1 ) work( k+n2 ) = -work( ki+n2 )*t( k, ki ) end do ! solve upper quasi-triangular system: ! (t(1:ki-2,1:ki-2) - (wr+i*wi))*x = scale*(work+i*work2) jnxt = ki - 2_${ik}$ loop_90: do j = ki - 2, 1, -1 if( j>jnxt )cycle loop_90 j1 = j j2 = j jnxt = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnxt = j - 2_${ik}$ end if end if if( j1==j2 ) then ! 1-by-1 diagonal block call stdlib${ii}$_slaln2( .false., 1_${ik}$, 2_${ik}$, smin, one, t( j, j ),ldt, one, one, & work( j+n ), n, wr, wi,x, 2_${ik}$, scale, xnorm, ierr ) ! scale x(1,1) and x(1,2) to avoid overflow when ! updating the right-hand side. if( xnorm>one ) then if( work( j )>bignum / xnorm ) then x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ ) / xnorm x( 1_${ik}$, 2_${ik}$ ) = x( 1_${ik}$, 2_${ik}$ ) / xnorm scale = scale / xnorm end if end if ! scale if necessary if( scale/=one ) then call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ ) call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n2 ), 1_${ik}$ ) end if work( j+n ) = x( 1_${ik}$, 1_${ik}$ ) work( j+n2 ) = x( 1_${ik}$, 2_${ik}$ ) ! update the right-hand side call stdlib${ii}$_saxpy( j-1, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ ) call stdlib${ii}$_saxpy( j-1, -x( 1_${ik}$, 2_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n2 ), 1_${ik}$ ) else ! 2-by-2 diagonal block call stdlib${ii}$_slaln2( .false., 2_${ik}$, 2_${ik}$, smin, one,t( j-1, j-1 ), ldt, one, & one,work( j-1+n ), n, wr, wi, x, 2_${ik}$, scale,xnorm, ierr ) ! scale x to avoid overflow when updating ! the right-hand side. if( xnorm>one ) then beta = max( work( j-1 ), work( j ) ) if( beta>bignum / xnorm ) then rec = one / xnorm x( 1_${ik}$, 1_${ik}$ ) = x( 1_${ik}$, 1_${ik}$ )*rec x( 1_${ik}$, 2_${ik}$ ) = x( 1_${ik}$, 2_${ik}$ )*rec x( 2_${ik}$, 1_${ik}$ ) = x( 2_${ik}$, 1_${ik}$ )*rec x( 2_${ik}$, 2_${ik}$ ) = x( 2_${ik}$, 2_${ik}$ )*rec scale = scale*rec end if end if ! scale if necessary if( scale/=one ) then call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n ), 1_${ik}$ ) call stdlib${ii}$_sscal( ki, scale, work( 1_${ik}$+n2 ), 1_${ik}$ ) end if work( j-1+n ) = x( 1_${ik}$, 1_${ik}$ ) work( j+n ) = x( 2_${ik}$, 1_${ik}$ ) work( j-1+n2 ) = x( 1_${ik}$, 2_${ik}$ ) work( j+n2 ) = x( 2_${ik}$, 2_${ik}$ ) ! update the right-hand side call stdlib${ii}$_saxpy( j-2, -x( 1_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j-1 ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ ) call stdlib${ii}$_saxpy( j-2, -x( 2_${ik}$, 1_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n ), 1_${ik}$ ) call stdlib${ii}$_saxpy( j-2, -x( 1_${ik}$, 2_${ik}$ ), t( 1_${ik}$, j-1 ), 1_${ik}$,work( 1_${ik}$+n2 ), 1_${ik}$ ) call stdlib${ii}$_saxpy( j-2, -x( 2_${ik}$, 2_${ik}$ ), t( 1_${ik}$, j ), 1_${ik}$,work( 1_${ik}$+n2 ), 1_${ik}$ ) end if end do loop_90 ! copy the vector x or q*x to vr and normalize. if( .not.over ) then call stdlib${ii}$_scopy( ki, work( 1_${ik}$+n ), 1_${ik}$, vr( 1_${ik}$, is-1 ), 1_${ik}$ ) call stdlib${ii}$_scopy( ki, work( 1_${ik}$+n2 ), 1_${ik}$, vr( 1_${ik}$, is ), 1_${ik}$ ) emax = zero do k = 1, ki emax = max( emax, abs( vr( k, is-1 ) )+abs( vr( k, is ) ) ) end do remax = one / emax call stdlib${ii}$_sscal( ki, remax, vr( 1_${ik}$, is-1 ), 1_${ik}$ ) call stdlib${ii}$_sscal( ki, remax, vr( 1_${ik}$, is ), 1_${ik}$ ) do k = ki + 1, n vr( k, is-1 ) = zero vr( k, is ) = zero end do else if( ki>2_${ik}$ ) then call stdlib${ii}$_sgemv( 'N', n, ki-2, one, vr, ldvr,work( 1_${ik}$+n ), 1_${ik}$, work( ki-& 1_${ik}$+n ),vr( 1_${ik}$, ki-1 ), 1_${ik}$ ) call stdlib${ii}$_sgemv( 'N', n, ki-2, one, vr, ldvr,work( 1_${ik}$+n2 ), 1_${ik}$, work( & ki+n2 ),vr( 1_${ik}$, ki ), 1_${ik}$ ) else call stdlib${ii}$_sscal( n, work( ki-1+n ), vr( 1_${ik}$, ki-1 ), 1_${ik}$ ) call stdlib${ii}$_sscal( n, work( ki+n2 ), vr( 1_${ik}$, ki ), 1_${ik}$ ) end if emax = zero do k = 1, n emax = max( emax, abs( vr( k, ki-1 ) )+abs( vr( k, ki ) ) ) end do remax = one / emax call stdlib${ii}$_sscal( n, remax, vr( 1_${ik}$, ki-1 ), 1_${ik}$ ) call stdlib${ii}$_sscal( n, remax, vr( 1_${ik}$, ki ), 1_${ik}$ ) end if end if is = is - 1_${ik}$ if( ip/=0_${ik}$ )is = is - 1_${ik}$ 130 continue if( ip==1_${ik}$ )ip = 0_${ik}$ if( ip==-1_${ik}$ )ip = 1_${ik}$ end do loop_140 end if if( leftv ) then ! compute left eigenvectors. ip = 0_${ik}$ is = 1_${ik}$ loop_260: do ki = 1, n if( ip==-1 )go to 250 if( ki==n )go to 150 if( t( ki+1, ki )==zero )go to 150 ip = 1_${ik}$ 150 continue if( somev ) then if( .not.select( ki ) )go to 250 end if ! compute the ki-th eigenvalue (wr,wi). wr = t( ki, ki ) wi = zero if( ip/=0_${ik}$ )wi = sqrt( abs( t( ki, ki+1 ) ) )*sqrt( abs( t( ki+1, ki ) ) ) smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum ) if( ip==0_${ik}$ ) then ! real left eigenvector. work( ki+n ) = one ! form right-hand side do k = ki + 1, n work( k+n ) = -t( ki, k ) end do ! solve the quasi-triangular system: ! (t(ki+1:n,ki+1:n) - wr)**t*x = scale*work vmax = one vcrit = bignum jnxt = ki + 1_${ik}$ loop_170: do j = ki + 1, n if( j<jnxt )cycle loop_170 j1 = j j2 = j jnxt = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnxt = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1-by-1 diagonal block ! scale if necessary to avoid overflow when forming ! the right-hand side. if( work( j )>vcrit ) then rec = one / vmax call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ ) vmax = one vcrit = bignum end if work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-1, t( ki+1, j ), 1_${ik}$,work( & ki+1+n ), 1_${ik}$ ) ! solve (t(j,j)-wr)**t*x = work call stdlib${ii}$_slaln2( .false., 1_${ik}$, 1_${ik}$, smin, one, t( j, j ),ldt, one, one, & work( j+n ), n, wr,zero, x, 2_${ik}$, scale, xnorm, ierr ) ! scale if necessary if( scale/=one )call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ ) work( j+n ) = x( 1_${ik}$, 1_${ik}$ ) vmax = max( abs( work( j+n ) ), vmax ) vcrit = bignum / vmax else ! 2-by-2 diagonal block ! scale if necessary to avoid overflow when forming ! the right-hand side. beta = max( work( j ), work( j+1 ) ) if( beta>vcrit ) then rec = one / vmax call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ ) vmax = one vcrit = bignum end if work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-1, t( ki+1, j ), 1_${ik}$,work( & ki+1+n ), 1_${ik}$ ) work( j+1+n ) = work( j+1+n ) -stdlib${ii}$_sdot( j-ki-1, t( ki+1, j+1 ), 1_${ik}$,& work( ki+1+n ), 1_${ik}$ ) ! solve ! [t(j,j)-wr t(j,j+1) ]**t* x = scale*( work1 ) ! [t(j+1,j) t(j+1,j+1)-wr] ( work2 ) call stdlib${ii}$_slaln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j, j ),ldt, one, one, & work( j+n ), n, wr,zero, x, 2_${ik}$, scale, xnorm, ierr ) ! scale if necessary if( scale/=one )call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ ) work( j+n ) = x( 1_${ik}$, 1_${ik}$ ) work( j+1+n ) = x( 2_${ik}$, 1_${ik}$ ) vmax = max( abs( work( j+n ) ),abs( work( j+1+n ) ), vmax ) vcrit = bignum / vmax end if end do loop_170 ! copy the vector x or q*x to vl and normalize. if( .not.over ) then call stdlib${ii}$_scopy( n-ki+1, work( ki+n ), 1_${ik}$, vl( ki, is ), 1_${ik}$ ) ii = stdlib${ii}$_isamax( n-ki+1, vl( ki, is ), 1_${ik}$ ) + ki - 1_${ik}$ remax = one / abs( vl( ii, is ) ) call stdlib${ii}$_sscal( n-ki+1, remax, vl( ki, is ), 1_${ik}$ ) do k = 1, ki - 1 vl( k, is ) = zero end do else if( ki<n )call stdlib${ii}$_sgemv( 'N', n, n-ki, one, vl( 1_${ik}$, ki+1 ), ldvl,work( & ki+1+n ), 1_${ik}$, work( ki+n ),vl( 1_${ik}$, ki ), 1_${ik}$ ) ii = stdlib${ii}$_isamax( n, vl( 1_${ik}$, ki ), 1_${ik}$ ) remax = one / abs( vl( ii, ki ) ) call stdlib${ii}$_sscal( n, remax, vl( 1_${ik}$, ki ), 1_${ik}$ ) end if else ! complex left eigenvector. ! initial solve: ! ((t(ki,ki) t(ki,ki+1) )**t - (wr - i* wi))*x = 0. ! ((t(ki+1,ki) t(ki+1,ki+1)) ) if( abs( t( ki, ki+1 ) )>=abs( t( ki+1, ki ) ) ) then work( ki+n ) = wi / t( ki, ki+1 ) work( ki+1+n2 ) = one else work( ki+n ) = one work( ki+1+n2 ) = -wi / t( ki+1, ki ) end if work( ki+1+n ) = zero work( ki+n2 ) = zero ! form right-hand side do k = ki + 2, n work( k+n ) = -work( ki+n )*t( ki, k ) work( k+n2 ) = -work( ki+1+n2 )*t( ki+1, k ) end do ! solve complex quasi-triangular system: ! ( t(ki+2,n:ki+2,n) - (wr-i*wi) )*x = work1+i*work2 vmax = one vcrit = bignum jnxt = ki + 2_${ik}$ loop_200: do j = ki + 2, n if( j<jnxt )cycle loop_200 j1 = j j2 = j jnxt = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnxt = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1-by-1 diagonal block ! scale if necessary to avoid overflow when ! forming the right-hand side elements. if( work( j )>vcrit ) then rec = one / vmax call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n2 ), 1_${ik}$ ) vmax = one vcrit = bignum end if work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( & ki+2+n ), 1_${ik}$ ) work( j+n2 ) = work( j+n2 ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( & ki+2+n2 ), 1_${ik}$ ) ! solve (t(j,j)-(wr-i*wi))*(x11+i*x12)= wk+i*wk2 call stdlib${ii}$_slaln2( .false., 1_${ik}$, 2_${ik}$, smin, one, t( j, j ),ldt, one, one, & work( j+n ), n, wr,-wi, x, 2_${ik}$, scale, xnorm, ierr ) ! scale if necessary if( scale/=one ) then call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n2 ), 1_${ik}$ ) end if work( j+n ) = x( 1_${ik}$, 1_${ik}$ ) work( j+n2 ) = x( 1_${ik}$, 2_${ik}$ ) vmax = max( abs( work( j+n ) ),abs( work( j+n2 ) ), vmax ) vcrit = bignum / vmax else ! 2-by-2 diagonal block ! scale if necessary to avoid overflow when forming ! the right-hand side elements. beta = max( work( j ), work( j+1 ) ) if( beta>vcrit ) then rec = one / vmax call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-ki+1, rec, work( ki+n2 ), 1_${ik}$ ) vmax = one vcrit = bignum end if work( j+n ) = work( j+n ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( & ki+2+n ), 1_${ik}$ ) work( j+n2 ) = work( j+n2 ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j ), 1_${ik}$,work( & ki+2+n2 ), 1_${ik}$ ) work( j+1+n ) = work( j+1+n ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j+1 ), 1_${ik}$,& work( ki+2+n ), 1_${ik}$ ) work( j+1+n2 ) = work( j+1+n2 ) -stdlib${ii}$_sdot( j-ki-2, t( ki+2, j+1 ), 1_${ik}$,& work( ki+2+n2 ), 1_${ik}$ ) ! solve 2-by-2 complex linear equation ! ([t(j,j) t(j,j+1) ]**t-(wr-i*wi)*i)*x = scale*b ! ([t(j+1,j) t(j+1,j+1)] ) call stdlib${ii}$_slaln2( .true., 2_${ik}$, 2_${ik}$, smin, one, t( j, j ),ldt, one, one, & work( j+n ), n, wr,-wi, x, 2_${ik}$, scale, xnorm, ierr ) ! scale if necessary if( scale/=one ) then call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-ki+1, scale, work( ki+n2 ), 1_${ik}$ ) end if work( j+n ) = x( 1_${ik}$, 1_${ik}$ ) work( j+n2 ) = x( 1_${ik}$, 2_${ik}$ ) work( j+1+n ) = x( 2_${ik}$, 1_${ik}$ ) work( j+1+n2 ) = x( 2_${ik}$, 2_${ik}$ ) vmax = max( abs( x( 1_${ik}$, 1_${ik}$ ) ), abs( x( 1_${ik}$, 2_${ik}$ ) ),abs( x( 2_${ik}$, 1_${ik}$ ) ), abs( x(& 2_${ik}$, 2_${ik}$ ) ), vmax ) vcrit = bignum / vmax end if end do loop_200 ! copy the vector x or q*x to vl and normalize. if( .not.over ) then call stdlib${ii}$_scopy( n-ki+1, work( ki+n ), 1_${ik}$, vl( ki, is ), 1_${ik}$ ) call stdlib${ii}$_scopy( n-ki+1, work( ki+n2 ), 1_${ik}$, vl( ki, is+1 ),1_${ik}$ ) emax = zero do k = ki, n emax = max( emax, abs( vl( k, is ) )+abs( vl( k, is+1 ) ) ) end do remax = one / emax call stdlib${ii}$_sscal( n-ki+1, remax, vl( ki, is ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-ki+1, remax, vl( ki, is+1 ), 1_${ik}$ ) do k = 1, ki - 1 vl( k, is ) = zero vl( k, is+1 ) = zero end do else if( ki<n-1 ) then call stdlib${ii}$_sgemv( 'N', n, n-ki-1, one, vl( 1_${ik}$, ki+2 ),ldvl, work( ki+2+& n ), 1_${ik}$, work( ki+n ),vl( 1_${ik}$, ki ), 1_${ik}$ ) call stdlib${ii}$_sgemv( 'N', n, n-ki-1, one, vl( 1_${ik}$, ki+2 ),ldvl, work( ki+2+& n2 ), 1_${ik}$,work( ki+1+n2 ), vl( 1_${ik}$, ki+1 ), 1_${ik}$ ) else call stdlib${ii}$_sscal( n, work( ki+n ), vl( 1_${ik}$, ki ), 1_${ik}$ ) call stdlib${ii}$_sscal( n, work( ki+1+n2 ), vl( 1_${ik}$, ki+1 ), 1_${ik}$ ) end if emax = zero do k = 1, n emax = max( emax, abs( vl( k, ki ) )+abs( vl( k, ki+1 ) ) ) end do remax = one / emax call stdlib${ii}$_sscal( n, remax, vl( 1_${ik}$, ki ), 1_${ik}$ ) call stdlib${ii}$_sscal( n, remax, vl( 1_${ik}$, ki+1 ), 1_${ik}$ ) end if end if is = is + 1_${ik}$ if( ip/=0_${ik}$ )is = is + 1_${ik}$ 250 continue if( ip==-1_${ik}$ )ip = 0_${ik}$ if( ip==1_${ik}$ )ip = -1_${ik}$ end do loop_260 end if return end subroutine stdlib${ii}$_strevc pure module subroutine stdlib${ii}$_dtrevc( side, howmny, select, n, t, ldt, vl, ldvl, vr,ldvr, mm, m, & !! DTREVC computes some or all of the right and/or left eigenvectors of !! a real upper quasi-triangular matrix T. !! Matrices of this type are produced by the Schur factorization of !! a real general matrix: A = Q*T*Q**T, as computed by DHSEQR. !! The right eigenvector x and the left eigenvector y of T corresponding !! to an eigenvalue w are defined by: !! T*x = w*x, (y**H)*T = w*(y**H) !! where y**H denotes the conjugate transpose of y. !! The eigenvalues are not input to this routine, but are read directly !! from the diagonal blocks of T. !! This routine returns the matrices X and/or Y of right and left !! eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an !! input matrix. If Q is the orthogonal factor that reduces a matrix !! A to Schur form T, then Q*X and Q*Y are the matrices of right and !! left eigenvectors of A. work, info ) ! -- lapack computational 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) :: howmny, side integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: ldt, ldvl, ldvr, mm, n ! Array Arguments logical(lk), intent(inout) :: select(*) real(dp), intent(in) :: t(ldt,*) real(dp), intent(inout) :: vl(ldvl,*), vr(ldvr,*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: allv, bothv, leftv, over, pair, rightv, somev integer(${ik}$) :: i, ierr, ii, ip, is, j, j1, j2, jnxt, k, ki, n2 real(dp) :: beta, bignum, emax, ovfl, rec, remax, scale, smin, smlnum, ulp, unfl, & vcrit, vmax, wi, wr, xnorm ! Intrinsic Functions ! Local Arrays real(dp) :: x(2_${ik}$,2_${ik}$) ! Executable Statements ! decode and test the input parameters bothv = stdlib_lsame( side, 'B' ) rightv = stdlib_lsame( side, 'R' ) .or. bothv leftv = stdlib_lsame( side, 'L' ) .or. bothv allv = stdlib_lsame( howmny, 'A' ) over = stdlib_lsame( howmny, 'B' ) somev = stdlib_lsame( howmny, 'S' ) info = 0_${ik}$ if( .not.rightv .and. .not.leftv ) then info = -1_${ik}$ else if( .not.allv .and. .not.over .and. .not.somev ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -4_${ik}$ else if( ldt<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldvl<1_${ik}$ .or. ( leftv .and. ldvl<n ) ) then info = -8_${ik}$ else if( ldvr<1_${ik}$ .or. ( rightv .and. ldvr<n ) ) then info =