stdlib_lapack_eigv_comp2.fypp Source File


Source Code

#:include "common.fypp" 
submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_comp2
  implicit none


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_stgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alphar, alphai, &
     !! STGSEN reorders the generalized real Schur decomposition of a real
     !! matrix pair (A, B) (in terms of an orthonormal equivalence trans-
     !! formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues
     !! appears in the leading diagonal blocks of the upper quasi-triangular
     !! matrix A and the upper triangular B. The leading columns of Q and
     !! Z form orthonormal bases of the corresponding left and right eigen-
     !! spaces (deflating subspaces). (A, B) must be in generalized real
     !! Schur canonical form (as returned by SGGES), i.e. A is block upper
     !! triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper
     !! triangular.
     !! STGSEN also computes the generalized eigenvalues
     !! w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
     !! of the reordered matrix pair (A, B).
     !! Optionally, STGSEN computes the estimates of reciprocal condition
     !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
     !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
     !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to
     !! the selected cluster and the eigenvalues outside the cluster, resp.,
     !! and norms of "projections" onto left and right eigenspaces w.r.t.
     !! the selected cluster in the (1,1)-block.
               beta, q, ldq, z, ldz, m, pl,pr, dif, work, lwork, iwork, liwork, 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 
           logical(lk), intent(in) :: wantq, wantz
           integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(sp), intent(out) :: pl, pr
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
           real(sp), intent(out) :: alphai(*), alphar(*), beta(*), dif(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: idifjb = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, pair, swap, wantd, wantd1, wantd2, wantp
           integer(${ik}$) :: i, ierr, ijb, k, kase, kk, ks, liwmin, lwmin, mn2, n1, n2
           real(sp) :: dscale, dsum, eps, rdscal, smlnum
           ! Local Arrays 
           integer(${ik}$) :: isave(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then
              info = -14_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'STGSEN', -info )
              return
           end if
           ! get machine constants
           eps = stdlib${ii}$_slamch( 'P' )
           smlnum = stdlib${ii}$_slamch( 'S' ) / eps
           ierr = 0_${ik}$
           wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$
           wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$
           wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$
           wantd = wantd1 .or. wantd2
           ! set m to the dimension of the specified pair of deflating
           ! subspaces.
           m = 0_${ik}$
           pair = .false.
           if( .not.lquery .or. ijob/=0_${ik}$ ) then
           do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 if( k<n ) then
                    if( a( k+1, k )==zero ) then
                       if( select( k ) )m = m + 1_${ik}$
                    else
                       pair = .true.
                       if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$
                    end if
                 else
                    if( select( n ) )m = m + 1_${ik}$
                 end if
              end if
           end do
           end if
           if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 2_${ik}$*m*(n-m) )
              liwmin = max( 1_${ik}$, n+6 )
           else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 4_${ik}$*m*(n-m) )
              liwmin = max( 1_${ik}$, 2_${ik}$*m*(n-m), n+6 )
           else
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16 )
              liwmin = 1_${ik}$
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           if( lwork<lwmin .and. .not.lquery ) then
              info = -22_${ik}$
           else if( liwork<liwmin .and. .not.lquery ) then
              info = -24_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'STGSEN', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible.
           if( m==n .or. m==0_${ik}$ ) then
              if( wantp ) then
                 pl = one
                 pr = one
              end if
              if( wantd ) then
                 dscale = zero
                 dsum = one
                 do i = 1, n
                    call stdlib${ii}$_slassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                    call stdlib${ii}$_slassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                 end do
                 dif( 1_${ik}$ ) = dscale*sqrt( dsum )
                 dif( 2_${ik}$ ) = dif( 1_${ik}$ )
              end if
              go to 60
           end if
           ! collect the selected blocks at the top-left corner of (a, b).
           ks = 0_${ik}$
           pair = .false.
           loop_30: do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 swap = select( k )
                 if( k<n ) then
                    if( a( k+1, k )/=zero ) then
                       pair = .true.
                       swap = swap .or. select( k+1 )
                    end if
                 end if
                 if( swap ) then
                    ks = ks + 1_${ik}$
                    ! swap the k-th block to position ks.
                    ! perform the reordering of diagonal blocks in (a, b)
                    ! by orthogonal transformation matrices and update
                    ! q and z accordingly (if requested):
                    kk = k
                    if( k/=ks )call stdlib${ii}$_stgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,z, ldz,&
                               kk, ks, work, lwork, ierr )
                    if( ierr>0_${ik}$ ) then
                       ! swap is rejected: exit.
                       info = 1_${ik}$
                       if( wantp ) then
                          pl = zero
                          pr = zero
                       end if
                       if( wantd ) then
                          dif( 1_${ik}$ ) = zero
                          dif( 2_${ik}$ ) = zero
                       end if
                       go to 60
                    end if
                    if( pair )ks = ks + 1_${ik}$
                 end if
              end if
           end do loop_30
           if( wantp ) then
              ! solve generalized sylvester equation for r and l
              ! and compute pl and pr.
              n1 = m
              n2 = n - m
              i = n1 + 1_${ik}$
              ijb = 0_${ik}$
              call stdlib${ii}$_slacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 )
              call stdlib${ii}$_slacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 )
              call stdlib${ii}$_stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,&
               i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, &
                         iwork, ierr )
              ! estimate the reciprocal of norms of "projections" onto left
              ! and right eigenspaces.
              rdscal = zero
              dsum = one
              call stdlib${ii}$_slassq( n1*n2, work, 1_${ik}$, rdscal, dsum )
              pl = rdscal*sqrt( dsum )
              if( pl==zero ) then
                 pl = one
              else
                 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
              end if
              rdscal = zero
              dsum = one
              call stdlib${ii}$_slassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum )
              pr = rdscal*sqrt( dsum )
              if( pr==zero ) then
                 pr = one
              else
                 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
              end if
           end if
           if( wantd ) then
              ! compute estimates of difu and difl.
              if( wantd1 ) then
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = idifjb
                 ! frobenius norm-based difu-estimate.
                 call stdlib${ii}$_stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(&
                  i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
                 ! frobenius norm-based difl-estimate.
                 call stdlib${ii}$_stgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),&
                  ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
              else
                 ! compute 1-norm-based estimates of difu and difl using
                 ! reversed communication with stdlib${ii}$_slacn2. in each step a
                 ! generalized sylvester equation or a transposed variant
                 ! is solved.
                 kase = 0_${ik}$
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = 0_${ik}$
                 mn2 = 2_${ik}$*n1*n2
                 ! 1-norm-based estimate of difu.
                 40 continue
                 call stdlib${ii}$_slacn2( mn2, work( mn2+1 ), work, iwork, dif( 1_${ik}$ ),kase, isave )
                           
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation.
                       call stdlib${ii}$_stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_stgsyl( 'T', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 40
                 end if
                 dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ )
                 ! 1-norm-based estimate of difl.
                 50 continue
                 call stdlib${ii}$_slacn2( mn2, work( mn2+1 ), work, iwork, dif( 2_${ik}$ ),kase, isave )
                           
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation.
                       call stdlib${ii}$_stgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_stgsyl( 'T', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 50
                 end if
                 dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ )
              end if
           end if
           60 continue
           ! compute generalized eigenvalues of reordered pair (a, b) and
           ! normalize the generalized schur form.
           pair = .false.
           loop_70: do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 if( k<n ) then
                    if( a( k+1, k )/=zero ) then
                       pair = .true.
                    end if
                 end if
                 if( pair ) then
                   ! compute the eigenvalue(s) at position k.
                    work( 1_${ik}$ ) = a( k, k )
                    work( 2_${ik}$ ) = a( k+1, k )
                    work( 3_${ik}$ ) = a( k, k+1 )
                    work( 4_${ik}$ ) = a( k+1, k+1 )
                    work( 5_${ik}$ ) = b( k, k )
                    work( 6_${ik}$ ) = b( k+1, k )
                    work( 7_${ik}$ ) = b( k, k+1 )
                    work( 8_${ik}$ ) = b( k+1, k+1 )
                    call stdlib${ii}$_slag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta( k ),beta( k+1 ), &
                              alphar( k ), alphar( k+1 ),alphai( k ) )
                    alphai( k+1 ) = -alphai( k )
                 else
                    if( sign( one, b( k, k ) )<zero ) then
                       ! if b(k,k) is negative, make it positive
                       do i = 1, n
                          a( k, i ) = -a( k, i )
                          b( k, i ) = -b( k, i )
                          if( wantq ) q( i, k ) = -q( i, k )
                       end do
                    end if
                    alphar( k ) = a( k, k )
                    alphai( k ) = zero
                    beta( k ) = b( k, k )
                 end if
              end if
           end do loop_70
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_stgsen

     pure module subroutine stdlib${ii}$_dtgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alphar, alphai, &
     !! DTGSEN reorders the generalized real Schur decomposition of a real
     !! matrix pair (A, B) (in terms of an orthonormal equivalence trans-
     !! formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues
     !! appears in the leading diagonal blocks of the upper quasi-triangular
     !! matrix A and the upper triangular B. The leading columns of Q and
     !! Z form orthonormal bases of the corresponding left and right eigen-
     !! spaces (deflating subspaces). (A, B) must be in generalized real
     !! Schur canonical form (as returned by DGGES), i.e. A is block upper
     !! triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper
     !! triangular.
     !! DTGSEN also computes the generalized eigenvalues
     !! w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
     !! of the reordered matrix pair (A, B).
     !! Optionally, DTGSEN computes the estimates of reciprocal condition
     !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
     !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
     !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to
     !! the selected cluster and the eigenvalues outside the cluster, resp.,
     !! and norms of "projections" onto left and right eigenspaces w.r.t.
     !! the selected cluster in the (1,1)-block.
               beta, q, ldq, z, ldz, m, pl,pr, dif, work, lwork, iwork, liwork, 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 
           logical(lk), intent(in) :: wantq, wantz
           integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(dp), intent(out) :: pl, pr
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
           real(dp), intent(out) :: alphai(*), alphar(*), beta(*), dif(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: idifjb = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, pair, swap, wantd, wantd1, wantd2, wantp
           integer(${ik}$) :: i, ierr, ijb, k, kase, kk, ks, liwmin, lwmin, mn2, n1, n2
           real(dp) :: dscale, dsum, eps, rdscal, smlnum
           ! Local Arrays 
           integer(${ik}$) :: isave(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then
              info = -14_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DTGSEN', -info )
              return
           end if
           ! get machine constants
           eps = stdlib${ii}$_dlamch( 'P' )
           smlnum = stdlib${ii}$_dlamch( 'S' ) / eps
           ierr = 0_${ik}$
           wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$
           wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$
           wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$
           wantd = wantd1 .or. wantd2
           ! set m to the dimension of the specified pair of deflating
           ! subspaces.
           m = 0_${ik}$
           pair = .false.
           if( .not.lquery .or. ijob/=0_${ik}$ ) then
           do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 if( k<n ) then
                    if( a( k+1, k )==zero ) then
                       if( select( k ) )m = m + 1_${ik}$
                    else
                       pair = .true.
                       if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$
                    end if
                 else
                    if( select( n ) )m = m + 1_${ik}$
                 end if
              end if
           end do
           end if
           if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 2_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, n+6 )
           else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 4_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+6 )
           else
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16 )
              liwmin = 1_${ik}$
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           if( lwork<lwmin .and. .not.lquery ) then
              info = -22_${ik}$
           else if( liwork<liwmin .and. .not.lquery ) then
              info = -24_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DTGSEN', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible.
           if( m==n .or. m==0_${ik}$ ) then
              if( wantp ) then
                 pl = one
                 pr = one
              end if
              if( wantd ) then
                 dscale = zero
                 dsum = one
                 do i = 1, n
                    call stdlib${ii}$_dlassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                    call stdlib${ii}$_dlassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                 end do
                 dif( 1_${ik}$ ) = dscale*sqrt( dsum )
                 dif( 2_${ik}$ ) = dif( 1_${ik}$ )
              end if
              go to 60
           end if
           ! collect the selected blocks at the top-left corner of (a, b).
           ks = 0_${ik}$
           pair = .false.
           loop_30: do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 swap = select( k )
                 if( k<n ) then
                    if( a( k+1, k )/=zero ) then
                       pair = .true.
                       swap = swap .or. select( k+1 )
                    end if
                 end if
                 if( swap ) then
                    ks = ks + 1_${ik}$
                    ! swap the k-th block to position ks.
                    ! perform the reordering of diagonal blocks in (a, b)
                    ! by orthogonal transformation matrices and update
                    ! q and z accordingly (if requested):
                    kk = k
                    if( k/=ks )call stdlib${ii}$_dtgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,z, ldz,&
                               kk, ks, work, lwork, ierr )
                    if( ierr>0_${ik}$ ) then
                       ! swap is rejected: exit.
                       info = 1_${ik}$
                       if( wantp ) then
                          pl = zero
                          pr = zero
                       end if
                       if( wantd ) then
                          dif( 1_${ik}$ ) = zero
                          dif( 2_${ik}$ ) = zero
                       end if
                       go to 60
                    end if
                    if( pair )ks = ks + 1_${ik}$
                 end if
              end if
           end do loop_30
           if( wantp ) then
              ! solve generalized sylvester equation for r and l
              ! and compute pl and pr.
              n1 = m
              n2 = n - m
              i = n1 + 1_${ik}$
              ijb = 0_${ik}$
              call stdlib${ii}$_dlacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 )
              call stdlib${ii}$_dlacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 )
              call stdlib${ii}$_dtgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,&
               i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, &
                         iwork, ierr )
              ! estimate the reciprocal of norms of "projections" onto left
              ! and right eigenspaces.
              rdscal = zero
              dsum = one
              call stdlib${ii}$_dlassq( n1*n2, work, 1_${ik}$, rdscal, dsum )
              pl = rdscal*sqrt( dsum )
              if( pl==zero ) then
                 pl = one
              else
                 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
              end if
              rdscal = zero
              dsum = one
              call stdlib${ii}$_dlassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum )
              pr = rdscal*sqrt( dsum )
              if( pr==zero ) then
                 pr = one
              else
                 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
              end if
           end if
           if( wantd ) then
              ! compute estimates of difu and difl.
              if( wantd1 ) then
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = idifjb
                 ! frobenius norm-based difu-estimate.
                 call stdlib${ii}$_dtgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(&
                  i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
                 ! frobenius norm-based difl-estimate.
                 call stdlib${ii}$_dtgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),&
                  ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
              else
                 ! compute 1-norm-based estimates of difu and difl using
                 ! reversed communication with stdlib${ii}$_dlacn2. in each step a
                 ! generalized sylvester equation or a transposed variant
                 ! is solved.
                 kase = 0_${ik}$
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = 0_${ik}$
                 mn2 = 2_${ik}$*n1*n2
                 ! 1-norm-based estimate of difu.
                 40 continue
                 call stdlib${ii}$_dlacn2( mn2, work( mn2+1 ), work, iwork, dif( 1_${ik}$ ),kase, isave )
                           
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation.
                       call stdlib${ii}$_dtgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_dtgsyl( 'T', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 40
                 end if
                 dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ )
                 ! 1-norm-based estimate of difl.
                 50 continue
                 call stdlib${ii}$_dlacn2( mn2, work( mn2+1 ), work, iwork, dif( 2_${ik}$ ),kase, isave )
                           
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation.
                       call stdlib${ii}$_dtgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_dtgsyl( 'T', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 50
                 end if
                 dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ )
              end if
           end if
           60 continue
           ! compute generalized eigenvalues of reordered pair (a, b) and
           ! normalize the generalized schur form.
           pair = .false.
           loop_80: do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 if( k<n ) then
                    if( a( k+1, k )/=zero ) then
                       pair = .true.
                    end if
                 end if
                 if( pair ) then
                   ! compute the eigenvalue(s) at position k.
                    work( 1_${ik}$ ) = a( k, k )
                    work( 2_${ik}$ ) = a( k+1, k )
                    work( 3_${ik}$ ) = a( k, k+1 )
                    work( 4_${ik}$ ) = a( k+1, k+1 )
                    work( 5_${ik}$ ) = b( k, k )
                    work( 6_${ik}$ ) = b( k+1, k )
                    work( 7_${ik}$ ) = b( k, k+1 )
                    work( 8_${ik}$ ) = b( k+1, k+1 )
                    call stdlib${ii}$_dlag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta( k ),beta( k+1 ), &
                              alphar( k ), alphar( k+1 ),alphai( k ) )
                    alphai( k+1 ) = -alphai( k )
                 else
                    if( sign( one, b( k, k ) )<zero ) then
                       ! if b(k,k) is negative, make it positive
                       do i = 1, n
                          a( k, i ) = -a( k, i )
                          b( k, i ) = -b( k, i )
                          if( wantq ) q( i, k ) = -q( i, k )
                       end do
                    end if
                    alphar( k ) = a( k, k )
                    alphai( k ) = zero
                    beta( k ) = b( k, k )
                 end if
              end if
           end do loop_80
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_dtgsen

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$tgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alphar, alphai, &
     !! DTGSEN: reorders the generalized real Schur decomposition of a real
     !! matrix pair (A, B) (in terms of an orthonormal equivalence trans-
     !! formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues
     !! appears in the leading diagonal blocks of the upper quasi-triangular
     !! matrix A and the upper triangular B. The leading columns of Q and
     !! Z form orthonormal bases of the corresponding left and right eigen-
     !! spaces (deflating subspaces). (A, B) must be in generalized real
     !! Schur canonical form (as returned by DGGES), i.e. A is block upper
     !! triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper
     !! triangular.
     !! DTGSEN also computes the generalized eigenvalues
     !! w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
     !! of the reordered matrix pair (A, B).
     !! Optionally, DTGSEN computes the estimates of reciprocal condition
     !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
     !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
     !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to
     !! the selected cluster and the eigenvalues outside the cluster, resp.,
     !! and norms of "projections" onto left and right eigenspaces w.r.t.
     !! the selected cluster in the (1,1)-block.
               beta, q, ldq, z, ldz, m, pl,pr, dif, work, lwork, iwork, liwork, 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 
           logical(lk), intent(in) :: wantq, wantz
           integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(${rk}$), intent(out) :: pl, pr
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
           real(${rk}$), intent(out) :: alphai(*), alphar(*), beta(*), dif(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: idifjb = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, pair, swap, wantd, wantd1, wantd2, wantp
           integer(${ik}$) :: i, ierr, ijb, k, kase, kk, ks, liwmin, lwmin, mn2, n1, n2
           real(${rk}$) :: dscale, dsum, eps, rdscal, smlnum
           ! Local Arrays 
           integer(${ik}$) :: isave(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then
              info = -14_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -16_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DTGSEN', -info )
              return
           end if
           ! get machine constants
           eps = stdlib${ii}$_${ri}$lamch( 'P' )
           smlnum = stdlib${ii}$_${ri}$lamch( 'S' ) / eps
           ierr = 0_${ik}$
           wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$
           wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$
           wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$
           wantd = wantd1 .or. wantd2
           ! set m to the dimension of the specified pair of deflating
           ! subspaces.
           m = 0_${ik}$
           pair = .false.
           if( .not.lquery .or. ijob/=0_${ik}$ ) then
           do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 if( k<n ) then
                    if( a( k+1, k )==zero ) then
                       if( select( k ) )m = m + 1_${ik}$
                    else
                       pair = .true.
                       if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$
                    end if
                 else
                    if( select( n ) )m = m + 1_${ik}$
                 end if
              end if
           end do
           end if
           if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 2_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, n+6 )
           else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 4_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+6 )
           else
              lwmin = max( 1_${ik}$, 4_${ik}$*n+16 )
              liwmin = 1_${ik}$
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           if( lwork<lwmin .and. .not.lquery ) then
              info = -22_${ik}$
           else if( liwork<liwmin .and. .not.lquery ) then
              info = -24_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DTGSEN', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible.
           if( m==n .or. m==0_${ik}$ ) then
              if( wantp ) then
                 pl = one
                 pr = one
              end if
              if( wantd ) then
                 dscale = zero
                 dsum = one
                 do i = 1, n
                    call stdlib${ii}$_${ri}$lassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                    call stdlib${ii}$_${ri}$lassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                 end do
                 dif( 1_${ik}$ ) = dscale*sqrt( dsum )
                 dif( 2_${ik}$ ) = dif( 1_${ik}$ )
              end if
              go to 60
           end if
           ! collect the selected blocks at the top-left corner of (a, b).
           ks = 0_${ik}$
           pair = .false.
           loop_30: do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 swap = select( k )
                 if( k<n ) then
                    if( a( k+1, k )/=zero ) then
                       pair = .true.
                       swap = swap .or. select( k+1 )
                    end if
                 end if
                 if( swap ) then
                    ks = ks + 1_${ik}$
                    ! swap the k-th block to position ks.
                    ! perform the reordering of diagonal blocks in (a, b)
                    ! by orthogonal transformation matrices and update
                    ! q and z accordingly (if requested):
                    kk = k
                    if( k/=ks )call stdlib${ii}$_${ri}$tgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,z, ldz,&
                               kk, ks, work, lwork, ierr )
                    if( ierr>0_${ik}$ ) then
                       ! swap is rejected: exit.
                       info = 1_${ik}$
                       if( wantp ) then
                          pl = zero
                          pr = zero
                       end if
                       if( wantd ) then
                          dif( 1_${ik}$ ) = zero
                          dif( 2_${ik}$ ) = zero
                       end if
                       go to 60
                    end if
                    if( pair )ks = ks + 1_${ik}$
                 end if
              end if
           end do loop_30
           if( wantp ) then
              ! solve generalized sylvester equation for r and l
              ! and compute pl and pr.
              n1 = m
              n2 = n - m
              i = n1 + 1_${ik}$
              ijb = 0_${ik}$
              call stdlib${ii}$_${ri}$lacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 )
              call stdlib${ii}$_${ri}$lacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 )
              call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,&
               i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, &
                         iwork, ierr )
              ! estimate the reciprocal of norms of "projections" onto left
              ! and right eigenspaces.
              rdscal = zero
              dsum = one
              call stdlib${ii}$_${ri}$lassq( n1*n2, work, 1_${ik}$, rdscal, dsum )
              pl = rdscal*sqrt( dsum )
              if( pl==zero ) then
                 pl = one
              else
                 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
              end if
              rdscal = zero
              dsum = one
              call stdlib${ii}$_${ri}$lassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum )
              pr = rdscal*sqrt( dsum )
              if( pr==zero ) then
                 pr = one
              else
                 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
              end if
           end if
           if( wantd ) then
              ! compute estimates of difu and difl.
              if( wantd1 ) then
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = idifjb
                 ! frobenius norm-based difu-estimate.
                 call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(&
                  i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
                 ! frobenius norm-based difl-estimate.
                 call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),&
                  ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
              else
                 ! compute 1-norm-based estimates of difu and difl using
                 ! reversed communication with stdlib${ii}$_${ri}$lacn2. in each step a
                 ! generalized sylvester equation or a transposed variant
                 ! is solved.
                 kase = 0_${ik}$
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = 0_${ik}$
                 mn2 = 2_${ik}$*n1*n2
                 ! 1-norm-based estimate of difu.
                 40 continue
                 call stdlib${ii}$_${ri}$lacn2( mn2, work( mn2+1 ), work, iwork, dif( 1_${ik}$ ),kase, isave )
                           
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation.
                       call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_${ri}$tgsyl( 'T', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 40
                 end if
                 dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ )
                 ! 1-norm-based estimate of difl.
                 50 continue
                 call stdlib${ii}$_${ri}$lacn2( mn2, work( mn2+1 ), work, iwork, dif( 2_${ik}$ ),kase, isave )
                           
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation.
                       call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_${ri}$tgsyl( 'T', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 50
                 end if
                 dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ )
              end if
           end if
           60 continue
           ! compute generalized eigenvalues of reordered pair (a, b) and
           ! normalize the generalized schur form.
           pair = .false.
           loop_80: do k = 1, n
              if( pair ) then
                 pair = .false.
              else
                 if( k<n ) then
                    if( a( k+1, k )/=zero ) then
                       pair = .true.
                    end if
                 end if
                 if( pair ) then
                   ! compute the eigenvalue(s) at position k.
                    work( 1_${ik}$ ) = a( k, k )
                    work( 2_${ik}$ ) = a( k+1, k )
                    work( 3_${ik}$ ) = a( k, k+1 )
                    work( 4_${ik}$ ) = a( k+1, k+1 )
                    work( 5_${ik}$ ) = b( k, k )
                    work( 6_${ik}$ ) = b( k+1, k )
                    work( 7_${ik}$ ) = b( k, k+1 )
                    work( 8_${ik}$ ) = b( k+1, k+1 )
                    call stdlib${ii}$_${ri}$lag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta( k ),beta( k+1 ), &
                              alphar( k ), alphar( k+1 ),alphai( k ) )
                    alphai( k+1 ) = -alphai( k )
                 else
                    if( sign( one, b( k, k ) )<zero ) then
                       ! if b(k,k) is negative, make it positive
                       do i = 1, n
                          a( k, i ) = -a( k, i )
                          b( k, i ) = -b( k, i )
                          if( wantq ) q( i, k ) = -q( i, k )
                       end do
                    end if
                    alphar( k ) = a( k, k )
                    alphai( k ) = zero
                    beta( k ) = b( k, k )
                 end if
              end if
           end do loop_80
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_${ri}$tgsen

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_ctgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alpha, beta, q, &
     !! CTGSEN reorders the generalized Schur decomposition of a complex
     !! matrix pair (A, B) (in terms of an unitary equivalence trans-
     !! formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues
     !! appears in the leading diagonal blocks of the pair (A,B). The leading
     !! columns of Q and Z form unitary bases of the corresponding left and
     !! right eigenspaces (deflating subspaces). (A, B) must be in
     !! generalized Schur canonical form, that is, A and B are both upper
     !! triangular.
     !! CTGSEN also computes the generalized eigenvalues
     !! w(j)= ALPHA(j) / BETA(j)
     !! of the reordered matrix pair (A, B).
     !! Optionally, the routine computes estimates of reciprocal condition
     !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
     !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
     !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to
     !! the selected cluster and the eigenvalues outside the cluster, resp.,
     !! and norms of "projections" onto left and right eigenspaces w.r.t.
     !! the selected cluster in the (1,1)-block.
               ldq, z, ldz, m, pl, pr, dif,work, lwork, iwork, liwork, 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 
           logical(lk), intent(in) :: wantq, wantz
           integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(sp), intent(out) :: pl, pr
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(out) :: dif(*)
           complex(sp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
           complex(sp), intent(out) :: alpha(*), beta(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: idifjb = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, swap, wantd, wantd1, wantd2, wantp
           integer(${ik}$) :: i, ierr, ijb, k, kase, ks, liwmin, lwmin, mn2, n1, n2
           real(sp) :: dscale, dsum, rdscal, safmin
           complex(sp) :: temp1, temp2
           ! Local Arrays 
           integer(${ik}$) :: isave(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then
              info = -13_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -15_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CTGSEN', -info )
              return
           end if
           ierr = 0_${ik}$
           wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$
           wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$
           wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$
           wantd = wantd1 .or. wantd2
           ! set m to the dimension of the specified pair of deflating
           ! subspaces.
           m = 0_${ik}$
           if( .not.lquery .or. ijob/=0_${ik}$ ) then
           do k = 1, n
              alpha( k ) = a( k, k )
              beta( k ) = b( k, k )
              if( k<n ) then
                 if( select( k ) )m = m + 1_${ik}$
              else
                 if( select( n ) )m = m + 1_${ik}$
              end if
           end do
           end if
           if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then
              lwmin = max( 1_${ik}$, 2_${ik}$*m*(n-m) )
              liwmin = max( 1_${ik}$, n+2 )
           else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*m*(n-m) )
              liwmin = max( 1_${ik}$, 2_${ik}$*m*(n-m), n+2 )
           else
              lwmin = 1_${ik}$
              liwmin = 1_${ik}$
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           if( lwork<lwmin .and. .not.lquery ) then
              info = -21_${ik}$
           else if( liwork<liwmin .and. .not.lquery ) then
              info = -23_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CTGSEN', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible.
           if( m==n .or. m==0_${ik}$ ) then
              if( wantp ) then
                 pl = one
                 pr = one
              end if
              if( wantd ) then
                 dscale = zero
                 dsum = one
                 do i = 1, n
                    call stdlib${ii}$_classq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                    call stdlib${ii}$_classq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                 end do
                 dif( 1_${ik}$ ) = dscale*sqrt( dsum )
                 dif( 2_${ik}$ ) = dif( 1_${ik}$ )
              end if
              go to 70
           end if
           ! get machine constant
           safmin = stdlib${ii}$_slamch( 'S' )
           ! collect the selected blocks at the top-left corner of (a, b).
           ks = 0_${ik}$
           do k = 1, n
              swap = select( k )
              if( swap ) then
                 ks = ks + 1_${ik}$
                 ! swap the k-th block to position ks. compute unitary q
                 ! and z that will swap adjacent diagonal blocks in (a, b).
                 if( k/=ks )call stdlib${ii}$_ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,ldz, k,&
                            ks, ierr )
                 if( ierr>0_${ik}$ ) then
                    ! swap is rejected: exit.
                    info = 1_${ik}$
                    if( wantp ) then
                       pl = zero
                       pr = zero
                    end if
                    if( wantd ) then
                       dif( 1_${ik}$ ) = zero
                       dif( 2_${ik}$ ) = zero
                    end if
                    go to 70
                 end if
              end if
           end do
           if( wantp ) then
              ! solve generalized sylvester equation for r and l:
                         ! a11 * r - l * a22 = a12
                         ! b11 * r - l * b22 = b12
              n1 = m
              n2 = n - m
              i = n1 + 1_${ik}$
              call stdlib${ii}$_clacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 )
              call stdlib${ii}$_clacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 )
              ijb = 0_${ik}$
              call stdlib${ii}$_ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,&
               i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, &
                         iwork, ierr )
              ! estimate the reciprocal of norms of "projections" onto
              ! left and right eigenspaces
              rdscal = zero
              dsum = one
              call stdlib${ii}$_classq( n1*n2, work, 1_${ik}$, rdscal, dsum )
              pl = rdscal*sqrt( dsum )
              if( pl==zero ) then
                 pl = one
              else
                 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
              end if
              rdscal = zero
              dsum = one
              call stdlib${ii}$_classq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum )
              pr = rdscal*sqrt( dsum )
              if( pr==zero ) then
                 pr = one
              else
                 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
              end if
           end if
           if( wantd ) then
              ! compute estimates difu and difl.
              if( wantd1 ) then
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = idifjb
                 ! frobenius norm-based difu estimate.
                 call stdlib${ii}$_ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(&
                  i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
                 ! frobenius norm-based difl estimate.
                 call stdlib${ii}$_ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),&
                  ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
              else
                 ! compute 1-norm-based estimates of difu and difl using
                 ! reversed communication with stdlib${ii}$_clacn2. in each step a
                 ! generalized sylvester equation or a transposed variant
                 ! is solved.
                 kase = 0_${ik}$
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = 0_${ik}$
                 mn2 = 2_${ik}$*n1*n2
                 ! 1-norm-based estimate of difu.
                 40 continue
                 call stdlib${ii}$_clacn2( mn2, work( mn2+1 ), work, dif( 1_${ik}$ ), kase,isave )
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation
                       call stdlib${ii}$_ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_ctgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 40
                 end if
                 dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ )
                 ! 1-norm-based estimate of difl.
                 50 continue
                 call stdlib${ii}$_clacn2( mn2, work( mn2+1 ), work, dif( 2_${ik}$ ), kase,isave )
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation
                       call stdlib${ii}$_ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_ctgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 50
                 end if
                 dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ )
              end if
           end if
           ! if b(k,k) is complex, make it real and positive (normalization
           ! of the generalized schur form) and store the generalized
           ! eigenvalues of reordered pair (a, b)
           do k = 1, n
              dscale = abs( b( k, k ) )
              if( dscale>safmin ) then
                 temp1 = conjg( b( k, k ) / dscale )
                 temp2 = b( k, k ) / dscale
                 b( k, k ) = dscale
                 call stdlib${ii}$_cscal( n-k, temp1, b( k, k+1 ), ldb )
                 call stdlib${ii}$_cscal( n-k+1, temp1, a( k, k ), lda )
                 if( wantq )call stdlib${ii}$_cscal( n, temp2, q( 1_${ik}$, k ), 1_${ik}$ )
              else
                 b( k, k ) = cmplx( zero, zero,KIND=sp)
              end if
              alpha( k ) = a( k, k )
              beta( k ) = b( k, k )
           end do
           70 continue
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_ctgsen

     pure module subroutine stdlib${ii}$_ztgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alpha, beta, q, &
     !! ZTGSEN reorders the generalized Schur decomposition of a complex
     !! matrix pair (A, B) (in terms of an unitary equivalence trans-
     !! formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues
     !! appears in the leading diagonal blocks of the pair (A,B). The leading
     !! columns of Q and Z form unitary bases of the corresponding left and
     !! right eigenspaces (deflating subspaces). (A, B) must be in
     !! generalized Schur canonical form, that is, A and B are both upper
     !! triangular.
     !! ZTGSEN also computes the generalized eigenvalues
     !! w(j)= ALPHA(j) / BETA(j)
     !! of the reordered matrix pair (A, B).
     !! Optionally, the routine computes estimates of reciprocal condition
     !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
     !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
     !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to
     !! the selected cluster and the eigenvalues outside the cluster, resp.,
     !! and norms of "projections" onto left and right eigenspaces w.r.t.
     !! the selected cluster in the (1,1)-block.
               ldq, z, ldz, m, pl, pr, dif,work, lwork, iwork, liwork, 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 
           logical(lk), intent(in) :: wantq, wantz
           integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(dp), intent(out) :: pl, pr
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(out) :: dif(*)
           complex(dp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
           complex(dp), intent(out) :: alpha(*), beta(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: idifjb = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, swap, wantd, wantd1, wantd2, wantp
           integer(${ik}$) :: i, ierr, ijb, k, kase, ks, liwmin, lwmin, mn2, n1, n2
           real(dp) :: dscale, dsum, rdscal, safmin
           complex(dp) :: temp1, temp2
           ! Local Arrays 
           integer(${ik}$) :: isave(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then
              info = -13_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -15_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZTGSEN', -info )
              return
           end if
           ierr = 0_${ik}$
           wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$
           wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$
           wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$
           wantd = wantd1 .or. wantd2
           ! set m to the dimension of the specified pair of deflating
           ! subspaces.
           m = 0_${ik}$
           if( .not.lquery .or. ijob/=0_${ik}$ ) then
           do k = 1, n
              alpha( k ) = a( k, k )
              beta( k ) = b( k, k )
              if( k<n ) then
                 if( select( k ) )m = m + 1_${ik}$
              else
                 if( select( n ) )m = m + 1_${ik}$
              end if
           end do
           end if
           if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then
              lwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, n+2 )
           else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+2 )
           else
              lwmin = 1_${ik}$
              liwmin = 1_${ik}$
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           if( lwork<lwmin .and. .not.lquery ) then
              info = -21_${ik}$
           else if( liwork<liwmin .and. .not.lquery ) then
              info = -23_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZTGSEN', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible.
           if( m==n .or. m==0_${ik}$ ) then
              if( wantp ) then
                 pl = one
                 pr = one
              end if
              if( wantd ) then
                 dscale = zero
                 dsum = one
                 do i = 1, n
                    call stdlib${ii}$_zlassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                    call stdlib${ii}$_zlassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                 end do
                 dif( 1_${ik}$ ) = dscale*sqrt( dsum )
                 dif( 2_${ik}$ ) = dif( 1_${ik}$ )
              end if
              go to 70
           end if
           ! get machine constant
           safmin = stdlib${ii}$_dlamch( 'S' )
           ! collect the selected blocks at the top-left corner of (a, b).
           ks = 0_${ik}$
           do k = 1, n
              swap = select( k )
              if( swap ) then
                 ks = ks + 1_${ik}$
                 ! swap the k-th block to position ks. compute unitary q
                 ! and z that will swap adjacent diagonal blocks in (a, b).
                 if( k/=ks )call stdlib${ii}$_ztgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,ldz, k,&
                            ks, ierr )
                 if( ierr>0_${ik}$ ) then
                    ! swap is rejected: exit.
                    info = 1_${ik}$
                    if( wantp ) then
                       pl = zero
                       pr = zero
                    end if
                    if( wantd ) then
                       dif( 1_${ik}$ ) = zero
                       dif( 2_${ik}$ ) = zero
                    end if
                    go to 70
                 end if
              end if
           end do
           if( wantp ) then
              ! solve generalized sylvester equation for r and l:
                         ! a11 * r - l * a22 = a12
                         ! b11 * r - l * b22 = b12
              n1 = m
              n2 = n - m
              i = n1 + 1_${ik}$
              call stdlib${ii}$_zlacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 )
              call stdlib${ii}$_zlacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 )
              ijb = 0_${ik}$
              call stdlib${ii}$_ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,&
               i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, &
                         iwork, ierr )
              ! estimate the reciprocal of norms of "projections" onto
              ! left and right eigenspaces
              rdscal = zero
              dsum = one
              call stdlib${ii}$_zlassq( n1*n2, work, 1_${ik}$, rdscal, dsum )
              pl = rdscal*sqrt( dsum )
              if( pl==zero ) then
                 pl = one
              else
                 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
              end if
              rdscal = zero
              dsum = one
              call stdlib${ii}$_zlassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum )
              pr = rdscal*sqrt( dsum )
              if( pr==zero ) then
                 pr = one
              else
                 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
              end if
           end if
           if( wantd ) then
              ! compute estimates difu and difl.
              if( wantd1 ) then
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = idifjb
                 ! frobenius norm-based difu estimate.
                 call stdlib${ii}$_ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(&
                  i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
                 ! frobenius norm-based difl estimate.
                 call stdlib${ii}$_ztgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),&
                  ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
              else
                 ! compute 1-norm-based estimates of difu and difl using
                 ! reversed communication with stdlib${ii}$_zlacn2. in each step a
                 ! generalized sylvester equation or a transposed variant
                 ! is solved.
                 kase = 0_${ik}$
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = 0_${ik}$
                 mn2 = 2_${ik}$*n1*n2
                 ! 1-norm-based estimate of difu.
                 40 continue
                 call stdlib${ii}$_zlacn2( mn2, work( mn2+1 ), work, dif( 1_${ik}$ ), kase,isave )
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation
                       call stdlib${ii}$_ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_ztgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 40
                 end if
                 dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ )
                 ! 1-norm-based estimate of difl.
                 50 continue
                 call stdlib${ii}$_zlacn2( mn2, work( mn2+1 ), work, dif( 2_${ik}$ ), kase,isave )
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation
                       call stdlib${ii}$_ztgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_ztgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 50
                 end if
                 dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ )
              end if
           end if
           ! if b(k,k) is complex, make it real and positive (normalization
           ! of the generalized schur form) and store the generalized
           ! eigenvalues of reordered pair (a, b)
           do k = 1, n
              dscale = abs( b( k, k ) )
              if( dscale>safmin ) then
                 temp1 = conjg( b( k, k ) / dscale )
                 temp2 = b( k, k ) / dscale
                 b( k, k ) = dscale
                 call stdlib${ii}$_zscal( n-k, temp1, b( k, k+1 ), ldb )
                 call stdlib${ii}$_zscal( n-k+1, temp1, a( k, k ), lda )
                 if( wantq )call stdlib${ii}$_zscal( n, temp2, q( 1_${ik}$, k ), 1_${ik}$ )
              else
                 b( k, k ) = cmplx( zero, zero,KIND=dp)
              end if
              alpha( k ) = a( k, k )
              beta( k ) = b( k, k )
           end do
           70 continue
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_ztgsen

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$tgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alpha, beta, q, &
     !! ZTGSEN: reorders the generalized Schur decomposition of a complex
     !! matrix pair (A, B) (in terms of an unitary equivalence trans-
     !! formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues
     !! appears in the leading diagonal blocks of the pair (A,B). The leading
     !! columns of Q and Z form unitary bases of the corresponding left and
     !! right eigenspaces (deflating subspaces). (A, B) must be in
     !! generalized Schur canonical form, that is, A and B are both upper
     !! triangular.
     !! ZTGSEN also computes the generalized eigenvalues
     !! w(j)= ALPHA(j) / BETA(j)
     !! of the reordered matrix pair (A, B).
     !! Optionally, the routine computes estimates of reciprocal condition
     !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11),
     !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s)
     !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to
     !! the selected cluster and the eigenvalues outside the cluster, resp.,
     !! and norms of "projections" onto left and right eigenspaces w.r.t.
     !! the selected cluster in the (1,1)-block.
               ldq, z, ldz, m, pl, pr, dif,work, lwork, iwork, liwork, 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 
           logical(lk), intent(in) :: wantq, wantz
           integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n
           integer(${ik}$), intent(out) :: info, m
           real(${ck}$), intent(out) :: pl, pr
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(${ck}$), intent(out) :: dif(*)
           complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
           complex(${ck}$), intent(out) :: alpha(*), beta(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: idifjb = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, swap, wantd, wantd1, wantd2, wantp
           integer(${ik}$) :: i, ierr, ijb, k, kase, ks, liwmin, lwmin, mn2, n1, n2
           real(${ck}$) :: dscale, dsum, rdscal, safmin
           complex(${ck}$) :: temp1, temp2
           ! Local Arrays 
           integer(${ik}$) :: isave(3_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ )
           if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -5_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -9_${ik}$
           else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then
              info = -13_${ik}$
           else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then
              info = -15_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZTGSEN', -info )
              return
           end if
           ierr = 0_${ik}$
           wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$
           wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$
           wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$
           wantd = wantd1 .or. wantd2
           ! set m to the dimension of the specified pair of deflating
           ! subspaces.
           m = 0_${ik}$
           if( .not.lquery .or. ijob/=0_${ik}$ ) then
           do k = 1, n
              alpha( k ) = a( k, k )
              beta( k ) = b( k, k )
              if( k<n ) then
                 if( select( k ) )m = m + 1_${ik}$
              else
                 if( select( n ) )m = m + 1_${ik}$
              end if
           end do
           end if
           if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then
              lwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, n+2 )
           else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then
              lwmin = max( 1_${ik}$, 4_${ik}$*m*( n-m ) )
              liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+2 )
           else
              lwmin = 1_${ik}$
              liwmin = 1_${ik}$
           end if
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           if( lwork<lwmin .and. .not.lquery ) then
              info = -21_${ik}$
           else if( liwork<liwmin .and. .not.lquery ) then
              info = -23_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZTGSEN', -info )
              return
           else if( lquery ) then
              return
           end if
           ! quick return if possible.
           if( m==n .or. m==0_${ik}$ ) then
              if( wantp ) then
                 pl = one
                 pr = one
              end if
              if( wantd ) then
                 dscale = zero
                 dsum = one
                 do i = 1, n
                    call stdlib${ii}$_${ci}$lassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                    call stdlib${ii}$_${ci}$lassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum )
                 end do
                 dif( 1_${ik}$ ) = dscale*sqrt( dsum )
                 dif( 2_${ik}$ ) = dif( 1_${ik}$ )
              end if
              go to 70
           end if
           ! get machine constant
           safmin = stdlib${ii}$_${c2ri(ci)}$lamch( 'S' )
           ! collect the selected blocks at the top-left corner of (a, b).
           ks = 0_${ik}$
           do k = 1, n
              swap = select( k )
              if( swap ) then
                 ks = ks + 1_${ik}$
                 ! swap the k-th block to position ks. compute unitary q
                 ! and z that will swap adjacent diagonal blocks in (a, b).
                 if( k/=ks )call stdlib${ii}$_${ci}$tgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,ldz, k,&
                            ks, ierr )
                 if( ierr>0_${ik}$ ) then
                    ! swap is rejected: exit.
                    info = 1_${ik}$
                    if( wantp ) then
                       pl = zero
                       pr = zero
                    end if
                    if( wantd ) then
                       dif( 1_${ik}$ ) = zero
                       dif( 2_${ik}$ ) = zero
                    end if
                    go to 70
                 end if
              end if
           end do
           if( wantp ) then
              ! solve generalized sylvester equation for r and l:
                         ! a11 * r - l * a22 = a12
                         ! b11 * r - l * b22 = b12
              n1 = m
              n2 = n - m
              i = n1 + 1_${ik}$
              call stdlib${ii}$_${ci}$lacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 )
              call stdlib${ii}$_${ci}$lacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 )
              ijb = 0_${ik}$
              call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,&
               i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, &
                         iwork, ierr )
              ! estimate the reciprocal of norms of "projections" onto
              ! left and right eigenspaces
              rdscal = zero
              dsum = one
              call stdlib${ii}$_${ci}$lassq( n1*n2, work, 1_${ik}$, rdscal, dsum )
              pl = rdscal*sqrt( dsum )
              if( pl==zero ) then
                 pl = one
              else
                 pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) )
              end if
              rdscal = zero
              dsum = one
              call stdlib${ii}$_${ci}$lassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum )
              pr = rdscal*sqrt( dsum )
              if( pr==zero ) then
                 pr = one
              else
                 pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) )
              end if
           end if
           if( wantd ) then
              ! compute estimates difu and difl.
              if( wantd1 ) then
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = idifjb
                 ! frobenius norm-based difu estimate.
                 call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(&
                  i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
                 ! frobenius norm-based difl estimate.
                 call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),&
                  ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-&
                            2_${ik}$*n1*n2, iwork, ierr )
              else
                 ! compute 1-norm-based estimates of difu and difl using
                 ! reversed communication with stdlib${ii}$_${ci}$lacn2. in each step a
                 ! generalized sylvester equation or a transposed variant
                 ! is solved.
                 kase = 0_${ik}$
                 n1 = m
                 n2 = n - m
                 i = n1 + 1_${ik}$
                 ijb = 0_${ik}$
                 mn2 = 2_${ik}$*n1*n2
                 ! 1-norm-based estimate of difu.
                 40 continue
                 call stdlib${ii}$_${ci}$lacn2( mn2, work( mn2+1 ), work, dif( 1_${ik}$ ), kase,isave )
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation
                       call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_${ci}$tgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 40
                 end if
                 dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ )
                 ! 1-norm-based estimate of difl.
                 50 continue
                 call stdlib${ii}$_${ci}$lacn2( mn2, work( mn2+1 ), work, dif( 2_${ik}$ ), kase,isave )
                 if( kase/=0_${ik}$ ) then
                    if( kase==1_${ik}$ ) then
                       ! solve generalized sylvester equation
                       call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( &
                       i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    else
                       ! solve the transposed variant.
                       call stdlib${ii}$_${ci}$tgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b, &
                       ldb, b( i, i ), ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )&
                                 , lwork-2*n1*n2, iwork,ierr )
                    end if
                    go to 50
                 end if
                 dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ )
              end if
           end if
           ! if b(k,k) is complex, make it real and positive (normalization
           ! of the generalized schur form) and store the generalized
           ! eigenvalues of reordered pair (a, b)
           do k = 1, n
              dscale = abs( b( k, k ) )
              if( dscale>safmin ) then
                 temp1 = conjg( b( k, k ) / dscale )
                 temp2 = b( k, k ) / dscale
                 b( k, k ) = dscale
                 call stdlib${ii}$_${ci}$scal( n-k, temp1, b( k, k+1 ), ldb )
                 call stdlib${ii}$_${ci}$scal( n-k+1, temp1, a( k, k ), lda )
                 if( wantq )call stdlib${ii}$_${ci}$scal( n, temp2, q( 1_${ik}$, k ), 1_${ik}$ )
              else
                 b( k, k ) = cmplx( zero, zero,KIND=${ck}$)
              end if
              alpha( k ) = a( k, k )
              beta( k ) = b( k, k )
           end do
           70 continue
           work( 1_${ik}$ ) = lwmin
           iwork( 1_${ik}$ ) = liwmin
           return
     end subroutine stdlib${ii}$_${ci}$tgsen

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_stgsna( job, howmny, select, n, a, lda, b, ldb, vl,ldvl, vr, ldvr, s, &
     !! STGSNA estimates reciprocal condition numbers for specified
     !! eigenvalues and/or eigenvectors of a matrix pair (A, B) in
     !! generalized real Schur canonical form (or of any matrix pair
     !! (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
     !! Z**T denotes the transpose of Z.
     !! (A, B) must be in generalized real Schur form (as returned by SGGES),
     !! i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
     !! blocks. B is upper triangular.
               dif, mm, m, work, lwork,iwork, 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, job
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: lda, ldb, ldvl, ldvr, lwork, mm, n
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(in) :: a(lda,*), b(ldb,*), vl(ldvl,*), vr(ldvr,*)
           real(sp), intent(out) :: dif(*), s(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: difdri = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, pair, somcon, wantbh, wantdf, wants
           integer(${ik}$) :: i, ierr, ifst, ilst, iz, k, ks, lwmin, n1, n2
           real(sp) :: alphai, alphar, alprqt, beta, c1, c2, cond, eps, lnrm, rnrm, root1, root2, &
                     scale, smlnum, tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv, uhbvi
           ! Local Arrays 
           real(sp) :: dummy(1_${ik}$), dummy1(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           wantbh = stdlib_lsame( job, 'B' )
           wants = stdlib_lsame( job, 'E' ) .or. wantbh
           wantdf = stdlib_lsame( job, 'V' ) .or. wantbh
           somcon = stdlib_lsame( howmny, 'S' )
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ )
           if( .not.wants .and. .not.wantdf ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( howmny, 'A' ) .and. .not.somcon ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           else if( wants .and. ldvl<n ) then
              info = -10_${ik}$
           else if( wants .and. ldvr<n ) then
              info = -12_${ik}$
           else
              ! set m to the number of eigenpairs for which condition numbers
              ! are required, and test mm.
              if( somcon ) then
                 m = 0_${ik}$
                 pair = .false.
                 do k = 1, n
                    if( pair ) then
                       pair = .false.
                    else
                       if( k<n ) then
                          if( a( k+1, k )==zero ) then
                             if( select( k ) )m = m + 1_${ik}$
                          else
                             pair = .true.
                             if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$
                          end if
                       else
                          if( select( n ) )m = m + 1_${ik}$
                       end if
                    end if
                 end do
              else
                 m = n
              end if
              if( n==0_${ik}$ ) then
                 lwmin = 1_${ik}$
              else if( stdlib_lsame( job, 'V' ) .or. stdlib_lsame( job, 'B' ) ) then
                 lwmin = 2_${ik}$*n*( n + 2_${ik}$ ) + 16_${ik}$
              else
                 lwmin = n
              end if
              work( 1_${ik}$ ) = lwmin
              if( mm<m ) then
                 info = -15_${ik}$
              else if( lwork<lwmin .and. .not.lquery ) then
                 info = -18_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'STGSNA', -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' ) / eps
           ks = 0_${ik}$
           pair = .false.
           loop_20: do k = 1, n
              ! determine whether a(k,k) begins a 1-by-1 or 2-by-2 block.
              if( pair ) then
                 pair = .false.
                 cycle loop_20
              else
                 if( k<n )pair = a( k+1, k )/=zero
              end if
              ! determine whether condition numbers are required for the k-th
              ! eigenpair.
              if( somcon ) then
                 if( pair ) then
                    if( .not.select( k ) .and. .not.select( k+1 ) )cycle loop_20
                 else
                    if( .not.select( k ) )cycle loop_20
                 end if
              end if
              ks = ks + 1_${ik}$
              if( wants ) then
                 ! compute the reciprocal condition number of the k-th
                 ! eigenvalue.
                 if( pair ) then
                    ! complex eigenvalue pair.
                    rnrm = stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, vr( &
                              1_${ik}$, ks+1 ), 1_${ik}$ ) )
                    lnrm = stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, vl( &
                              1_${ik}$, ks+1 ), 1_${ik}$ ) )
                    call stdlib${ii}$_sgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    tmprr = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    tmpri = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    call stdlib${ii}$_sgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ )
                              
                    tmpii = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    tmpir = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    uhav = tmprr + tmpii
                    uhavi = tmpir - tmpri
                    call stdlib${ii}$_sgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    tmprr = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    tmpri = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    call stdlib${ii}$_sgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ )
                              
                    tmpii = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    tmpir = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    uhbv = tmprr + tmpii
                    uhbvi = tmpir - tmpri
                    uhav = stdlib${ii}$_slapy2( uhav, uhavi )
                    uhbv = stdlib${ii}$_slapy2( uhbv, uhbvi )
                    cond = stdlib${ii}$_slapy2( uhav, uhbv )
                    s( ks ) = cond / ( rnrm*lnrm )
                    s( ks+1 ) = s( ks )
                 else
                    ! real eigenvalue.
                    rnrm = stdlib${ii}$_snrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ )
                    lnrm = stdlib${ii}$_snrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    call stdlib${ii}$_sgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    uhav = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    call stdlib${ii}$_sgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    uhbv = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    cond = stdlib${ii}$_slapy2( uhav, uhbv )
                    if( cond==zero ) then
                       s( ks ) = -one
                    else
                       s( ks ) = cond / ( rnrm*lnrm )
                    end if
                 end if
              end if
              if( wantdf ) then
                 if( n==1_${ik}$ ) then
                    dif( ks ) = stdlib${ii}$_slapy2( a( 1_${ik}$, 1_${ik}$ ), b( 1_${ik}$, 1_${ik}$ ) )
                    cycle loop_20
                 end if
                 ! estimate the reciprocal condition number of the k-th
                 ! eigenvectors.
                 if( pair ) then
                    ! copy the  2-by 2 pencil beginning at (a(k,k), b(k, k)).
                    ! compute the eigenvalue(s) at position k.
                    work( 1_${ik}$ ) = a( k, k )
                    work( 2_${ik}$ ) = a( k+1, k )
                    work( 3_${ik}$ ) = a( k, k+1 )
                    work( 4_${ik}$ ) = a( k+1, k+1 )
                    work( 5_${ik}$ ) = b( k, k )
                    work( 6_${ik}$ ) = b( k+1, k )
                    work( 7_${ik}$ ) = b( k, k+1 )
                    work( 8_${ik}$ ) = b( k+1, k+1 )
                    call stdlib${ii}$_slag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta,dummy1( 1_${ik}$ ), &
                              alphar, dummy( 1_${ik}$ ), alphai )
                    alprqt = one
                    c1 = two*( alphar*alphar+alphai*alphai+beta*beta )
                    c2 = four*beta*beta*alphai*alphai
                    root1 = c1 + sqrt( c1*c1-4.0_sp*c2 )
                    root2 = c2 / root1
                    root1 = root1 / two
                    cond = min( sqrt( root1 ), sqrt( root2 ) )
                 end if
                 ! copy the matrix (a, b) to the array work and swap the
                 ! diagonal block beginning at a(k,k) to the (1,1) position.
                 call stdlib${ii}$_slacpy( 'FULL', n, n, a, lda, work, n )
                 call stdlib${ii}$_slacpy( 'FULL', n, n, b, ldb, work( n*n+1 ), n )
                 ifst = k
                 ilst = 1_${ik}$
                 call stdlib${ii}$_stgexc( .false., .false., n, work, n, work( n*n+1 ), n,dummy, 1_${ik}$, &
                           dummy1, 1_${ik}$, ifst, ilst,work( n*n*2_${ik}$+1 ), lwork-2*n*n, ierr )
                 if( ierr>0_${ik}$ ) then
                    ! ill-conditioned problem - swap rejected.
                    dif( ks ) = zero
                 else
                    ! reordering successful, solve generalized sylvester
                    ! equation for r and l,
                               ! a22 * r - l * a11 = a12
                               ! b22 * r - l * b11 = b12,
                    ! and compute estimate of difl((a11,b11), (a22, b22)).
                    n1 = 1_${ik}$
                    if( work( 2_${ik}$ )/=zero )n1 = 2_${ik}$
                    n2 = n - n1
                    if( n2==0_${ik}$ ) then
                       dif( ks ) = cond
                    else
                       i = n*n + 1_${ik}$
                       iz = 2_${ik}$*n*n + 1_${ik}$
                       call stdlib${ii}$_stgsyl( 'N', difdri, n2, n1, work( n*n1+n1+1 ),n, work, n, &
                       work( n1+1 ), n,work( n*n1+n1+i ), n, work( i ), n,work( n1+i ), n, scale, &
                                 dif( ks ),work( iz+1 ), lwork-2*n*n, iwork, ierr )
                       if( pair )dif( ks ) = min( max( one, alprqt )*dif( ks ),cond )
                    end if
                 end if
                 if( pair )dif( ks+1 ) = dif( ks )
              end if
              if( pair )ks = ks + 1_${ik}$
           end do loop_20
           work( 1_${ik}$ ) = lwmin
           return
     end subroutine stdlib${ii}$_stgsna

     pure module subroutine stdlib${ii}$_dtgsna( job, howmny, select, n, a, lda, b, ldb, vl,ldvl, vr, ldvr, s, &
     !! DTGSNA estimates reciprocal condition numbers for specified
     !! eigenvalues and/or eigenvectors of a matrix pair (A, B) in
     !! generalized real Schur canonical form (or of any matrix pair
     !! (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
     !! Z**T denotes the transpose of Z.
     !! (A, B) must be in generalized real Schur form (as returned by DGGES),
     !! i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
     !! blocks. B is upper triangular.
               dif, mm, m, work, lwork,iwork, 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, job
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: lda, ldb, ldvl, ldvr, lwork, mm, n
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(in) :: a(lda,*), b(ldb,*), vl(ldvl,*), vr(ldvr,*)
           real(dp), intent(out) :: dif(*), s(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: difdri = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, pair, somcon, wantbh, wantdf, wants
           integer(${ik}$) :: i, ierr, ifst, ilst, iz, k, ks, lwmin, n1, n2
           real(dp) :: alphai, alphar, alprqt, beta, c1, c2, cond, eps, lnrm, rnrm, root1, root2, &
                     scale, smlnum, tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv, uhbvi
           ! Local Arrays 
           real(dp) :: dummy(1_${ik}$), dummy1(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           wantbh = stdlib_lsame( job, 'B' )
           wants = stdlib_lsame( job, 'E' ) .or. wantbh
           wantdf = stdlib_lsame( job, 'V' ) .or. wantbh
           somcon = stdlib_lsame( howmny, 'S' )
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ )
           if( .not.wants .and. .not.wantdf ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( howmny, 'A' ) .and. .not.somcon ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           else if( wants .and. ldvl<n ) then
              info = -10_${ik}$
           else if( wants .and. ldvr<n ) then
              info = -12_${ik}$
           else
              ! set m to the number of eigenpairs for which condition numbers
              ! are required, and test mm.
              if( somcon ) then
                 m = 0_${ik}$
                 pair = .false.
                 do k = 1, n
                    if( pair ) then
                       pair = .false.
                    else
                       if( k<n ) then
                          if( a( k+1, k )==zero ) then
                             if( select( k ) )m = m + 1_${ik}$
                          else
                             pair = .true.
                             if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$
                          end if
                       else
                          if( select( n ) )m = m + 1_${ik}$
                       end if
                    end if
                 end do
              else
                 m = n
              end if
              if( n==0_${ik}$ ) then
                 lwmin = 1_${ik}$
              else if( stdlib_lsame( job, 'V' ) .or. stdlib_lsame( job, 'B' ) ) then
                 lwmin = 2_${ik}$*n*( n + 2_${ik}$ ) + 16_${ik}$
              else
                 lwmin = n
              end if
              work( 1_${ik}$ ) = lwmin
              if( mm<m ) then
                 info = -15_${ik}$
              else if( lwork<lwmin .and. .not.lquery ) then
                 info = -18_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DTGSNA', -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' ) / eps
           ks = 0_${ik}$
           pair = .false.
           loop_20: do k = 1, n
              ! determine whether a(k,k) begins a 1-by-1 or 2-by-2 block.
              if( pair ) then
                 pair = .false.
                 cycle loop_20
              else
                 if( k<n )pair = a( k+1, k )/=zero
              end if
              ! determine whether condition numbers are required for the k-th
              ! eigenpair.
              if( somcon ) then
                 if( pair ) then
                    if( .not.select( k ) .and. .not.select( k+1 ) )cycle loop_20
                 else
                    if( .not.select( k ) )cycle loop_20
                 end if
              end if
              ks = ks + 1_${ik}$
              if( wants ) then
                 ! compute the reciprocal condition number of the k-th
                 ! eigenvalue.
                 if( pair ) then
                    ! complex eigenvalue pair.
                    rnrm = stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, vr( &
                              1_${ik}$, ks+1 ), 1_${ik}$ ) )
                    lnrm = stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, vl( &
                              1_${ik}$, ks+1 ), 1_${ik}$ ) )
                    call stdlib${ii}$_dgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    tmprr = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    tmpri = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    call stdlib${ii}$_dgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ )
                              
                    tmpii = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    tmpir = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    uhav = tmprr + tmpii
                    uhavi = tmpir - tmpri
                    call stdlib${ii}$_dgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    tmprr = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    tmpri = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    call stdlib${ii}$_dgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ )
                              
                    tmpii = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    tmpir = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    uhbv = tmprr + tmpii
                    uhbvi = tmpir - tmpri
                    uhav = stdlib${ii}$_dlapy2( uhav, uhavi )
                    uhbv = stdlib${ii}$_dlapy2( uhbv, uhbvi )
                    cond = stdlib${ii}$_dlapy2( uhav, uhbv )
                    s( ks ) = cond / ( rnrm*lnrm )
                    s( ks+1 ) = s( ks )
                 else
                    ! real eigenvalue.
                    rnrm = stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ )
                    lnrm = stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    call stdlib${ii}$_dgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    uhav = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    call stdlib${ii}$_dgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    uhbv = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    cond = stdlib${ii}$_dlapy2( uhav, uhbv )
                    if( cond==zero ) then
                       s( ks ) = -one
                    else
                       s( ks ) = cond / ( rnrm*lnrm )
                    end if
                 end if
              end if
              if( wantdf ) then
                 if( n==1_${ik}$ ) then
                    dif( ks ) = stdlib${ii}$_dlapy2( a( 1_${ik}$, 1_${ik}$ ), b( 1_${ik}$, 1_${ik}$ ) )
                    cycle loop_20
                 end if
                 ! estimate the reciprocal condition number of the k-th
                 ! eigenvectors.
                 if( pair ) then
                    ! copy the  2-by 2 pencil beginning at (a(k,k), b(k, k)).
                    ! compute the eigenvalue(s) at position k.
                    work( 1_${ik}$ ) = a( k, k )
                    work( 2_${ik}$ ) = a( k+1, k )
                    work( 3_${ik}$ ) = a( k, k+1 )
                    work( 4_${ik}$ ) = a( k+1, k+1 )
                    work( 5_${ik}$ ) = b( k, k )
                    work( 6_${ik}$ ) = b( k+1, k )
                    work( 7_${ik}$ ) = b( k, k+1 )
                    work( 8_${ik}$ ) = b( k+1, k+1 )
                    call stdlib${ii}$_dlag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta,dummy1( 1_${ik}$ ), &
                              alphar, dummy( 1_${ik}$ ), alphai )
                    alprqt = one
                    c1 = two*( alphar*alphar+alphai*alphai+beta*beta )
                    c2 = four*beta*beta*alphai*alphai
                    root1 = c1 + sqrt( c1*c1-4.0_dp*c2 )
                    root2 = c2 / root1
                    root1 = root1 / two
                    cond = min( sqrt( root1 ), sqrt( root2 ) )
                 end if
                 ! copy the matrix (a, b) to the array work and swap the
                 ! diagonal block beginning at a(k,k) to the (1,1) position.
                 call stdlib${ii}$_dlacpy( 'FULL', n, n, a, lda, work, n )
                 call stdlib${ii}$_dlacpy( 'FULL', n, n, b, ldb, work( n*n+1 ), n )
                 ifst = k
                 ilst = 1_${ik}$
                 call stdlib${ii}$_dtgexc( .false., .false., n, work, n, work( n*n+1 ), n,dummy, 1_${ik}$, &
                           dummy1, 1_${ik}$, ifst, ilst,work( n*n*2_${ik}$+1 ), lwork-2*n*n, ierr )
                 if( ierr>0_${ik}$ ) then
                    ! ill-conditioned problem - swap rejected.
                    dif( ks ) = zero
                 else
                    ! reordering successful, solve generalized sylvester
                    ! equation for r and l,
                               ! a22 * r - l * a11 = a12
                               ! b22 * r - l * b11 = b12,
                    ! and compute estimate of difl((a11,b11), (a22, b22)).
                    n1 = 1_${ik}$
                    if( work( 2_${ik}$ )/=zero )n1 = 2_${ik}$
                    n2 = n - n1
                    if( n2==0_${ik}$ ) then
                       dif( ks ) = cond
                    else
                       i = n*n + 1_${ik}$
                       iz = 2_${ik}$*n*n + 1_${ik}$
                       call stdlib${ii}$_dtgsyl( 'N', difdri, n2, n1, work( n*n1+n1+1 ),n, work, n, &
                       work( n1+1 ), n,work( n*n1+n1+i ), n, work( i ), n,work( n1+i ), n, scale, &
                                 dif( ks ),work( iz+1 ), lwork-2*n*n, iwork, ierr )
                       if( pair )dif( ks ) = min( max( one, alprqt )*dif( ks ),cond )
                    end if
                 end if
                 if( pair )dif( ks+1 ) = dif( ks )
              end if
              if( pair )ks = ks + 1_${ik}$
           end do loop_20
           work( 1_${ik}$ ) = lwmin
           return
     end subroutine stdlib${ii}$_dtgsna

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$tgsna( job, howmny, select, n, a, lda, b, ldb, vl,ldvl, vr, ldvr, s, &
     !! DTGSNA: estimates reciprocal condition numbers for specified
     !! eigenvalues and/or eigenvectors of a matrix pair (A, B) in
     !! generalized real Schur canonical form (or of any matrix pair
     !! (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where
     !! Z**T denotes the transpose of Z.
     !! (A, B) must be in generalized real Schur form (as returned by DGGES),
     !! i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
     !! blocks. B is upper triangular.
               dif, mm, m, work, lwork,iwork, 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) :: howmny, job
           integer(${ik}$), intent(out) :: info, m
           integer(${ik}$), intent(in) :: lda, ldb, ldvl, ldvr, lwork, mm, n
           ! Array Arguments 
           logical(lk), intent(in) :: select(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(in) :: a(lda,*), b(ldb,*), vl(ldvl,*), vr(ldvr,*)
           real(${rk}$), intent(out) :: dif(*), s(*), work(*)
        ! =====================================================================
           ! Parameters 
           integer(${ik}$), parameter :: difdri = 3_${ik}$
           
           
           ! Local Scalars 
           logical(lk) :: lquery, pair, somcon, wantbh, wantdf, wants
           integer(${ik}$) :: i, ierr, ifst, ilst, iz, k, ks, lwmin, n1, n2
           real(${rk}$) :: alphai, alphar, alprqt, beta, c1, c2, cond, eps, lnrm, rnrm, root1, root2, &
                     scale, smlnum, tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv, uhbvi
           ! Local Arrays 
           real(${rk}$) :: dummy(1_${ik}$), dummy1(1_${ik}$)
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters
           wantbh = stdlib_lsame( job, 'B' )
           wants = stdlib_lsame( job, 'E' ) .or. wantbh
           wantdf = stdlib_lsame( job, 'V' ) .or. wantbh
           somcon = stdlib_lsame( howmny, 'S' )
           info = 0_${ik}$
           lquery = ( lwork==-1_${ik}$ )
           if( .not.wants .and. .not.wantdf ) then
              info = -1_${ik}$
           else if( .not.stdlib_lsame( howmny, 'A' ) .and. .not.somcon ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -4_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -8_${ik}$
           else if( wants .and. ldvl<n ) then
              info = -10_${ik}$
           else if( wants .and. ldvr<n ) then
              info = -12_${ik}$
           else
              ! set m to the number of eigenpairs for which condition numbers
              ! are required, and test mm.
              if( somcon ) then
                 m = 0_${ik}$
                 pair = .false.
                 do k = 1, n
                    if( pair ) then
                       pair = .false.
                    else
                       if( k<n ) then
                          if( a( k+1, k )==zero ) then
                             if( select( k ) )m = m + 1_${ik}$
                          else
                             pair = .true.
                             if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$
                          end if
                       else
                          if( select( n ) )m = m + 1_${ik}$
                       end if
                    end if
                 end do
              else
                 m = n
              end if
              if( n==0_${ik}$ ) then
                 lwmin = 1_${ik}$
              else if( stdlib_lsame( job, 'V' ) .or. stdlib_lsame( job, 'B' ) ) then
                 lwmin = 2_${ik}$*n*( n + 2_${ik}$ ) + 16_${ik}$
              else
                 lwmin = n
              end if
              work( 1_${ik}$ ) = lwmin
              if( mm<m ) then
                 info = -15_${ik}$
              else if( lwork<lwmin .and. .not.lquery ) then
                 info = -18_${ik}$
              end if
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DTGSNA', -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' ) / eps
           ks = 0_${ik}$
           pair = .false.
           loop_20: do k = 1, n
              ! determine whether a(k,k) begins a 1-by-1 or 2-by-2 block.
              if( pair ) then
                 pair = .false.
                 cycle loop_20
              else
                 if( k<n )pair = a( k+1, k )/=zero
              end if
              ! determine whether condition numbers are required for the k-th
              ! eigenpair.
              if( somcon ) then
                 if( pair ) then
                    if( .not.select( k ) .and. .not.select( k+1 ) )cycle loop_20
                 else
                    if( .not.select( k ) )cycle loop_20
                 end if
              end if
              ks = ks + 1_${ik}$
              if( wants ) then
                 ! compute the reciprocal condition number of the k-th
                 ! eigenvalue.
                 if( pair ) then
                    ! complex eigenvalue pair.
                    rnrm = stdlib${ii}$_${ri}$lapy2( stdlib${ii}$_${ri}$nrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_${ri}$nrm2( n, vr( &
                              1_${ik}$, ks+1 ), 1_${ik}$ ) )
                    lnrm = stdlib${ii}$_${ri}$lapy2( stdlib${ii}$_${ri}$nrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_${ri}$nrm2( n, vl( &
                              1_${ik}$, ks+1 ), 1_${ik}$ ) )
                    call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    tmprr = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    tmpri = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ )
                              
                    tmpii = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    tmpir = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    uhav = tmprr + tmpii
                    uhavi = tmpir - tmpri
                    call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    tmprr = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    tmpri = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ )
                              
                    tmpii = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ )
                    tmpir = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    uhbv = tmprr + tmpii
                    uhbvi = tmpir - tmpri
                    uhav = stdlib${ii}$_${ri}$lapy2( uhav, uhavi )
                    uhbv = stdlib${ii}$_${ri}$lapy2( uhbv, uhbvi )
                    cond = stdlib${ii}$_${ri}$lapy2( uhav, uhbv )
                    s( ks ) = cond / ( rnrm*lnrm )
                    s( ks+1 ) = s( ks )
                 else
                    ! real eigenvalue.
                    rnrm = stdlib${ii}$_${ri}$nrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ )
                    lnrm = stdlib${ii}$_${ri}$nrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    uhav = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ )
                              
                    uhbv = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ )
                    cond = stdlib${ii}$_${ri}$lapy2( uhav, uhbv )
                    if( cond==zero ) then
                       s( ks ) = -one
                    else
                       s( ks ) = cond / ( rnrm*lnrm )
                    end if
                 end if
              end if
              if( wantdf ) then
                 if( n==1_${ik}$ ) then
                    dif( ks ) = stdlib${ii}$_${ri}$lapy2( a( 1_${ik}$, 1_${ik}$ ), b( 1_${ik}$, 1_${ik}$ ) )
                    cycle loop_20
                 end if
                 ! estimate the reciprocal condition number of the k-th
                 ! eigenvectors.
                 if( pair ) then
                    ! copy the  2-by 2 pencil beginning at (a(k,k), b(k, k)).
                    ! compute the eigenvalue(s) at position k.
                    work( 1_${ik}$ ) = a( k, k )
                    work( 2_${ik}$ ) = a( k+1, k )
                    work( 3_${ik}$ ) = a( k, k