stdlib_lapack_eigv_comp.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_sggbal( job, n, a, lda, b, ldb, ilo, ihi, lscale,rscale, work, info )
     !! SGGBAL balances a pair of general real matrices (A,B).  This
     !! involves, first, permuting A and B by similarity transformations to
     !! isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
     !! elements on the diagonal; and second, applying a diagonal similarity
     !! transformation to rows and columns ILO to IHI to make the rows
     !! and columns as close in norm as possible. Both steps are optional.
     !! Balancing may reduce the 1-norm of the matrices, and improve the
     !! accuracy of the computed eigenvalues and/or eigenvectors in the
     !! generalized eigenvalue problem A*x = lambda*B*x.
               
        ! -- 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) :: job
           integer(${ik}$), intent(out) :: ihi, ilo, info
           integer(${ik}$), intent(in) :: lda, ldb, n
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*), b(ldb,*)
           real(sp), intent(out) :: lscale(*), rscale(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: sclfac = ten
           
           
           ! Local Scalars 
           integer(${ik}$) :: i, icab, iflow, ip1, ir, irab, it, j, jc, jp1, k, kount, l, lcab, lm1, &
                     lrab, lsfmax, lsfmin, m, nr, nrp2
           real(sp) :: alpha, basl, beta, cab, cmax, coef, coef2, coef5, cor, ew, ewc, gamma, &
                     pgamma, rab, sfmax, sfmin, sum, t, ta, tb, tc
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'N' ) .and. .not.stdlib_lsame( job, 'P' ) &
                     .and..not.stdlib_lsame( job, 'S' ) .and. .not.stdlib_lsame( job, 'B' ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SGGBAL', -info )
              return
           end if
           ! quick return if possible
           if( n==0_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              return
           end if
           if( n==1_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              lscale( 1_${ik}$ ) = one
              rscale( 1_${ik}$ ) = one
              return
           end if
           if( stdlib_lsame( job, 'N' ) ) then
              ilo = 1_${ik}$
              ihi = n
              do i = 1, n
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           k = 1_${ik}$
           l = n
           if( stdlib_lsame( job, 'S' ) )go to 190
           go to 30
           ! permute the matrices a and b to isolate the eigenvalues.
           ! find row with one nonzero in columns 1 through l
           20 continue
           l = lm1
           if( l/=1 )go to 30
           rscale( 1_${ik}$ ) = one
           lscale( 1_${ik}$ ) = one
           go to 190
           30 continue
           lm1 = l - 1_${ik}$
           loop_80: do i = l, 1, -1
              do j = 1, lm1
                 jp1 = j + 1_${ik}$
                 if( a( i, j )/=zero .or. b( i, j )/=zero )go to 50
              end do
              j = l
              go to 70
              50 continue
              do j = jp1, l
                 if( a( i, j )/=zero .or. b( i, j )/=zero )cycle loop_80
              end do
              j = jp1 - 1_${ik}$
              70 continue
              m = l
              iflow = 1_${ik}$
              go to 160
           end do loop_80
           go to 100
           ! find column with one nonzero in rows k through n
           90 continue
           k = k + 1_${ik}$
           100 continue
           loop_150: do j = k, l
              do i = k, lm1
                 ip1 = i + 1_${ik}$
                 if( a( i, j )/=zero .or. b( i, j )/=zero )go to 120
              end do
              i = l
              go to 140
              120 continue
              do i = ip1, l
                 if( a( i, j )/=zero .or. b( i, j )/=zero )cycle loop_150
              end do
              i = ip1 - 1_${ik}$
              140 continue
              m = k
              iflow = 2_${ik}$
              go to 160
           end do loop_150
           go to 190
           ! permute rows m and i
           160 continue
           lscale( m ) = i
           if( i==m )go to 170
           call stdlib${ii}$_sswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
           call stdlib${ii}$_sswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
           ! permute columns m and j
           170 continue
           rscale( m ) = j
           if( j==m )go to 180
           call stdlib${ii}$_sswap( l, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, m ), 1_${ik}$ )
           call stdlib${ii}$_sswap( l, b( 1_${ik}$, j ), 1_${ik}$, b( 1_${ik}$, m ), 1_${ik}$ )
           180 continue
           go to ( 20, 90 )iflow
           190 continue
           ilo = k
           ihi = l
           if( stdlib_lsame( job, 'P' ) ) then
              do i = ilo, ihi
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           if( ilo==ihi )return
           ! balance the submatrix in rows ilo to ihi.
           nr = ihi - ilo + 1_${ik}$
           do i = ilo, ihi
              rscale( i ) = zero
              lscale( i ) = zero
              work( i ) = zero
              work( i+n ) = zero
              work( i+2*n ) = zero
              work( i+3*n ) = zero
              work( i+4*n ) = zero
              work( i+5*n ) = zero
           end do
           ! compute right side vector in resulting linear equations
           basl = log10( sclfac )
           do i = ilo, ihi
              do j = ilo, ihi
                 tb = b( i, j )
                 ta = a( i, j )
                 if( ta==zero )go to 210
                 ta = log10( abs( ta ) ) / basl
                 210 continue
                 if( tb==zero )go to 220
                 tb = log10( abs( tb ) ) / basl
                 220 continue
                 work( i+4*n ) = work( i+4*n ) - ta - tb
                 work( j+5*n ) = work( j+5*n ) - ta - tb
              end do
           end do
           coef = one / real( 2_${ik}$*nr,KIND=sp)
           coef2 = coef*coef
           coef5 = half*coef2
           nrp2 = nr + 2_${ik}$
           beta = zero
           it = 1_${ik}$
           ! start generalized conjugate gradient iteration
           250 continue
           gamma = stdlib${ii}$_sdot( nr, work( ilo+4*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ ) +stdlib${ii}$_sdot( nr, &
                     work( ilo+5*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           ew = zero
           ewc = zero
           do i = ilo, ihi
              ew = ew + work( i+4*n )
              ewc = ewc + work( i+5*n )
           end do
           gamma = coef*gamma - coef2*( ew**2_${ik}$+ewc**2_${ik}$ ) - coef5*( ew-ewc )**2_${ik}$
           if( gamma==zero )go to 350
           if( it/=1_${ik}$ )beta = gamma / pgamma
           t = coef5*( ewc-three*ew )
           tc = coef5*( ew-three*ewc )
           call stdlib${ii}$_sscal( nr, beta, work( ilo ), 1_${ik}$ )
           call stdlib${ii}$_sscal( nr, beta, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_saxpy( nr, coef, work( ilo+4*n ), 1_${ik}$, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_saxpy( nr, coef, work( ilo+5*n ), 1_${ik}$, work( ilo ), 1_${ik}$ )
           do i = ilo, ihi
              work( i ) = work( i ) + tc
              work( i+n ) = work( i+n ) + t
           end do
           ! apply matrix to vector
           do i = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_290: do j = ilo, ihi
                 if( a( i, j )==zero )go to 280
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
                 280 continue
                 if( b( i, j )==zero )cycle loop_290
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
              end do loop_290
              work( i+2*n ) = real( kount,KIND=sp)*work( i+n ) + sum
           end do
           do j = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_320: do i = ilo, ihi
                 if( a( i, j )==zero )go to 310
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
                 310 continue
                 if( b( i, j )==zero )cycle loop_320
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
              end do loop_320
              work( j+3*n ) = real( kount,KIND=sp)*work( j ) + sum
           end do
           sum = stdlib${ii}$_sdot( nr, work( ilo+n ), 1_${ik}$, work( ilo+2*n ), 1_${ik}$ ) +stdlib${ii}$_sdot( nr, work( &
                     ilo ), 1_${ik}$, work( ilo+3*n ), 1_${ik}$ )
           alpha = gamma / sum
           ! determine correction to current iteration
           cmax = zero
           do i = ilo, ihi
              cor = alpha*work( i+n )
              if( abs( cor )>cmax )cmax = abs( cor )
              lscale( i ) = lscale( i ) + cor
              cor = alpha*work( i )
              if( abs( cor )>cmax )cmax = abs( cor )
              rscale( i ) = rscale( i ) + cor
           end do
           if( cmax<half )go to 350
           call stdlib${ii}$_saxpy( nr, -alpha, work( ilo+2*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ )
           call stdlib${ii}$_saxpy( nr, -alpha, work( ilo+3*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           pgamma = gamma
           it = it + 1_${ik}$
           if( it<=nrp2 )go to 250
           ! end generalized conjugate gradient iteration
           350 continue
           sfmin = stdlib${ii}$_slamch( 'S' )
           sfmax = one / sfmin
           lsfmin = int( log10( sfmin ) / basl+one,KIND=${ik}$)
           lsfmax = int( log10( sfmax ) / basl,KIND=${ik}$)
           do i = ilo, ihi
              irab = stdlib${ii}$_isamax( n-ilo+1, a( i, ilo ), lda )
              rab = abs( a( i, irab+ilo-1 ) )
              irab = stdlib${ii}$_isamax( n-ilo+1, b( i, ilo ), ldb )
              rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
              lrab = int( log10( rab+sfmin ) / basl+one,KIND=${ik}$)
              ir = lscale( i ) + sign( half, lscale( i ) )
              ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
              lscale( i ) = sclfac**ir
              icab = stdlib${ii}$_isamax( ihi, a( 1_${ik}$, i ), 1_${ik}$ )
              cab = abs( a( icab, i ) )
              icab = stdlib${ii}$_isamax( ihi, b( 1_${ik}$, i ), 1_${ik}$ )
              cab = max( cab, abs( b( icab, i ) ) )
              lcab = int( log10( cab+sfmin ) / basl+one,KIND=${ik}$)
              jc = rscale( i ) + sign( half, rscale( i ) )
              jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
              rscale( i ) = sclfac**jc
           end do
           ! row scaling of matrices a and b
           do i = ilo, ihi
              call stdlib${ii}$_sscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
              call stdlib${ii}$_sscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
           end do
           ! column scaling of matrices a and b
           do j = ilo, ihi
              call stdlib${ii}$_sscal( ihi, rscale( j ), a( 1_${ik}$, j ), 1_${ik}$ )
              call stdlib${ii}$_sscal( ihi, rscale( j ), b( 1_${ik}$, j ), 1_${ik}$ )
           end do
           return
     end subroutine stdlib${ii}$_sggbal

     pure module subroutine stdlib${ii}$_dggbal( job, n, a, lda, b, ldb, ilo, ihi, lscale,rscale, work, info )
     !! DGGBAL balances a pair of general real matrices (A,B).  This
     !! involves, first, permuting A and B by similarity transformations to
     !! isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
     !! elements on the diagonal; and second, applying a diagonal similarity
     !! transformation to rows and columns ILO to IHI to make the rows
     !! and columns as close in norm as possible. Both steps are optional.
     !! Balancing may reduce the 1-norm of the matrices, and improve the
     !! accuracy of the computed eigenvalues and/or eigenvectors in the
     !! generalized eigenvalue problem A*x = lambda*B*x.
               
        ! -- 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) :: job
           integer(${ik}$), intent(out) :: ihi, ilo, info
           integer(${ik}$), intent(in) :: lda, ldb, n
           ! Array Arguments 
           real(dp), intent(inout) :: a(lda,*), b(ldb,*)
           real(dp), intent(out) :: lscale(*), rscale(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: sclfac = ten
           
           
           ! Local Scalars 
           integer(${ik}$) :: i, icab, iflow, ip1, ir, irab, it, j, jc, jp1, k, kount, l, lcab, lm1, &
                     lrab, lsfmax, lsfmin, m, nr, nrp2
           real(dp) :: alpha, basl, beta, cab, cmax, coef, coef2, coef5, cor, ew, ewc, gamma, &
                     pgamma, rab, sfmax, sfmin, sum, t, ta, tb, tc
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'N' ) .and. .not.stdlib_lsame( job, 'P' ) &
                     .and..not.stdlib_lsame( job, 'S' ) .and. .not.stdlib_lsame( job, 'B' ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGGBAL', -info )
              return
           end if
           ! quick return if possible
           if( n==0_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              return
           end if
           if( n==1_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              lscale( 1_${ik}$ ) = one
              rscale( 1_${ik}$ ) = one
              return
           end if
           if( stdlib_lsame( job, 'N' ) ) then
              ilo = 1_${ik}$
              ihi = n
              do i = 1, n
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           k = 1_${ik}$
           l = n
           if( stdlib_lsame( job, 'S' ) )go to 190
           go to 30
           ! permute the matrices a and b to isolate the eigenvalues.
           ! find row with one nonzero in columns 1 through l
           20 continue
           l = lm1
           if( l/=1 )go to 30
           rscale( 1_${ik}$ ) = one
           lscale( 1_${ik}$ ) = one
           go to 190
           30 continue
           lm1 = l - 1_${ik}$
           loop_80: do i = l, 1, -1
              do j = 1, lm1
                 jp1 = j + 1_${ik}$
                 if( a( i, j )/=zero .or. b( i, j )/=zero )go to 50
              end do
              j = l
              go to 70
              50 continue
              do j = jp1, l
                 if( a( i, j )/=zero .or. b( i, j )/=zero )cycle loop_80
              end do
              j = jp1 - 1_${ik}$
              70 continue
              m = l
              iflow = 1_${ik}$
              go to 160
           end do loop_80
           go to 100
           ! find column with one nonzero in rows k through n
           90 continue
           k = k + 1_${ik}$
           100 continue
           loop_150: do j = k, l
              do i = k, lm1
                 ip1 = i + 1_${ik}$
                 if( a( i, j )/=zero .or. b( i, j )/=zero )go to 120
              end do
              i = l
              go to 140
              120 continue
              do i = ip1, l
                 if( a( i, j )/=zero .or. b( i, j )/=zero )cycle loop_150
              end do
              i = ip1 - 1_${ik}$
              140 continue
              m = k
              iflow = 2_${ik}$
              go to 160
           end do loop_150
           go to 190
           ! permute rows m and i
           160 continue
           lscale( m ) = i
           if( i==m )go to 170
           call stdlib${ii}$_dswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
           call stdlib${ii}$_dswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
           ! permute columns m and j
           170 continue
           rscale( m ) = j
           if( j==m )go to 180
           call stdlib${ii}$_dswap( l, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, m ), 1_${ik}$ )
           call stdlib${ii}$_dswap( l, b( 1_${ik}$, j ), 1_${ik}$, b( 1_${ik}$, m ), 1_${ik}$ )
           180 continue
           go to ( 20, 90 )iflow
           190 continue
           ilo = k
           ihi = l
           if( stdlib_lsame( job, 'P' ) ) then
              do i = ilo, ihi
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           if( ilo==ihi )return
           ! balance the submatrix in rows ilo to ihi.
           nr = ihi - ilo + 1_${ik}$
           do i = ilo, ihi
              rscale( i ) = zero
              lscale( i ) = zero
              work( i ) = zero
              work( i+n ) = zero
              work( i+2*n ) = zero
              work( i+3*n ) = zero
              work( i+4*n ) = zero
              work( i+5*n ) = zero
           end do
           ! compute right side vector in resulting linear equations
           basl = log10( sclfac )
           do i = ilo, ihi
              do j = ilo, ihi
                 tb = b( i, j )
                 ta = a( i, j )
                 if( ta==zero )go to 210
                 ta = log10( abs( ta ) ) / basl
                 210 continue
                 if( tb==zero )go to 220
                 tb = log10( abs( tb ) ) / basl
                 220 continue
                 work( i+4*n ) = work( i+4*n ) - ta - tb
                 work( j+5*n ) = work( j+5*n ) - ta - tb
              end do
           end do
           coef = one / real( 2_${ik}$*nr,KIND=dp)
           coef2 = coef*coef
           coef5 = half*coef2
           nrp2 = nr + 2_${ik}$
           beta = zero
           it = 1_${ik}$
           ! start generalized conjugate gradient iteration
           250 continue
           gamma = stdlib${ii}$_ddot( nr, work( ilo+4*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ ) +stdlib${ii}$_ddot( nr, &
                     work( ilo+5*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           ew = zero
           ewc = zero
           do i = ilo, ihi
              ew = ew + work( i+4*n )
              ewc = ewc + work( i+5*n )
           end do
           gamma = coef*gamma - coef2*( ew**2_${ik}$+ewc**2_${ik}$ ) - coef5*( ew-ewc )**2_${ik}$
           if( gamma==zero )go to 350
           if( it/=1_${ik}$ )beta = gamma / pgamma
           t = coef5*( ewc-three*ew )
           tc = coef5*( ew-three*ewc )
           call stdlib${ii}$_dscal( nr, beta, work( ilo ), 1_${ik}$ )
           call stdlib${ii}$_dscal( nr, beta, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_daxpy( nr, coef, work( ilo+4*n ), 1_${ik}$, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_daxpy( nr, coef, work( ilo+5*n ), 1_${ik}$, work( ilo ), 1_${ik}$ )
           do i = ilo, ihi
              work( i ) = work( i ) + tc
              work( i+n ) = work( i+n ) + t
           end do
           ! apply matrix to vector
           do i = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_290: do j = ilo, ihi
                 if( a( i, j )==zero )go to 280
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
                 280 continue
                 if( b( i, j )==zero )cycle loop_290
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
              end do loop_290
              work( i+2*n ) = real( kount,KIND=dp)*work( i+n ) + sum
           end do
           do j = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_320: do i = ilo, ihi
                 if( a( i, j )==zero )go to 310
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
                 310 continue
                 if( b( i, j )==zero )cycle loop_320
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
              end do loop_320
              work( j+3*n ) = real( kount,KIND=dp)*work( j ) + sum
           end do
           sum = stdlib${ii}$_ddot( nr, work( ilo+n ), 1_${ik}$, work( ilo+2*n ), 1_${ik}$ ) +stdlib${ii}$_ddot( nr, work( &
                     ilo ), 1_${ik}$, work( ilo+3*n ), 1_${ik}$ )
           alpha = gamma / sum
           ! determine correction to current iteration
           cmax = zero
           do i = ilo, ihi
              cor = alpha*work( i+n )
              if( abs( cor )>cmax )cmax = abs( cor )
              lscale( i ) = lscale( i ) + cor
              cor = alpha*work( i )
              if( abs( cor )>cmax )cmax = abs( cor )
              rscale( i ) = rscale( i ) + cor
           end do
           if( cmax<half )go to 350
           call stdlib${ii}$_daxpy( nr, -alpha, work( ilo+2*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ )
           call stdlib${ii}$_daxpy( nr, -alpha, work( ilo+3*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           pgamma = gamma
           it = it + 1_${ik}$
           if( it<=nrp2 )go to 250
           ! end generalized conjugate gradient iteration
           350 continue
           sfmin = stdlib${ii}$_dlamch( 'S' )
           sfmax = one / sfmin
           lsfmin = int( log10( sfmin ) / basl+one,KIND=${ik}$)
           lsfmax = int( log10( sfmax ) / basl,KIND=${ik}$)
           do i = ilo, ihi
              irab = stdlib${ii}$_idamax( n-ilo+1, a( i, ilo ), lda )
              rab = abs( a( i, irab+ilo-1 ) )
              irab = stdlib${ii}$_idamax( n-ilo+1, b( i, ilo ), ldb )
              rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
              lrab = int( log10( rab+sfmin ) / basl+one,KIND=${ik}$)
              ir = int(lscale( i ) + sign( half, lscale( i ) ),KIND=${ik}$)
              ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
              lscale( i ) = sclfac**ir
              icab = stdlib${ii}$_idamax( ihi, a( 1_${ik}$, i ), 1_${ik}$ )
              cab = abs( a( icab, i ) )
              icab = stdlib${ii}$_idamax( ihi, b( 1_${ik}$, i ), 1_${ik}$ )
              cab = max( cab, abs( b( icab, i ) ) )
              lcab = int( log10( cab+sfmin ) / basl+one,KIND=${ik}$)
              jc = int(rscale( i ) + sign( half, rscale( i ) ),KIND=${ik}$)
              jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
              rscale( i ) = sclfac**jc
           end do
           ! row scaling of matrices a and b
           do i = ilo, ihi
              call stdlib${ii}$_dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
              call stdlib${ii}$_dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
           end do
           ! column scaling of matrices a and b
           do j = ilo, ihi
              call stdlib${ii}$_dscal( ihi, rscale( j ), a( 1_${ik}$, j ), 1_${ik}$ )
              call stdlib${ii}$_dscal( ihi, rscale( j ), b( 1_${ik}$, j ), 1_${ik}$ )
           end do
           return
     end subroutine stdlib${ii}$_dggbal

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$ggbal( job, n, a, lda, b, ldb, ilo, ihi, lscale,rscale, work, info )
     !! DGGBAL: balances a pair of general real matrices (A,B).  This
     !! involves, first, permuting A and B by similarity transformations to
     !! isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
     !! elements on the diagonal; and second, applying a diagonal similarity
     !! transformation to rows and columns ILO to IHI to make the rows
     !! and columns as close in norm as possible. Both steps are optional.
     !! Balancing may reduce the 1-norm of the matrices, and improve the
     !! accuracy of the computed eigenvalues and/or eigenvectors in the
     !! generalized eigenvalue problem A*x = lambda*B*x.
               
        ! -- 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) :: job
           integer(${ik}$), intent(out) :: ihi, ilo, info
           integer(${ik}$), intent(in) :: lda, ldb, n
           ! Array Arguments 
           real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*)
           real(${rk}$), intent(out) :: lscale(*), rscale(*), work(*)
        ! =====================================================================
           ! Parameters 
           real(${rk}$), parameter :: sclfac = ten
           
           
           ! Local Scalars 
           integer(${ik}$) :: i, icab, iflow, ip1, ir, irab, it, j, jc, jp1, k, kount, l, lcab, lm1, &
                     lrab, lsfmax, lsfmin, m, nr, nrp2
           real(${rk}$) :: alpha, basl, beta, cab, cmax, coef, coef2, coef5, cor, ew, ewc, gamma, &
                     pgamma, rab, sfmax, sfmin, sum, t, ta, tb, tc
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'N' ) .and. .not.stdlib_lsame( job, 'P' ) &
                     .and..not.stdlib_lsame( job, 'S' ) .and. .not.stdlib_lsame( job, 'B' ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGGBAL', -info )
              return
           end if
           ! quick return if possible
           if( n==0_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              return
           end if
           if( n==1_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              lscale( 1_${ik}$ ) = one
              rscale( 1_${ik}$ ) = one
              return
           end if
           if( stdlib_lsame( job, 'N' ) ) then
              ilo = 1_${ik}$
              ihi = n
              do i = 1, n
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           k = 1_${ik}$
           l = n
           if( stdlib_lsame( job, 'S' ) )go to 190
           go to 30
           ! permute the matrices a and b to isolate the eigenvalues.
           ! find row with one nonzero in columns 1 through l
           20 continue
           l = lm1
           if( l/=1 )go to 30
           rscale( 1_${ik}$ ) = one
           lscale( 1_${ik}$ ) = one
           go to 190
           30 continue
           lm1 = l - 1_${ik}$
           loop_80: do i = l, 1, -1
              do j = 1, lm1
                 jp1 = j + 1_${ik}$
                 if( a( i, j )/=zero .or. b( i, j )/=zero )go to 50
              end do
              j = l
              go to 70
              50 continue
              do j = jp1, l
                 if( a( i, j )/=zero .or. b( i, j )/=zero )cycle loop_80
              end do
              j = jp1 - 1_${ik}$
              70 continue
              m = l
              iflow = 1_${ik}$
              go to 160
           end do loop_80
           go to 100
           ! find column with one nonzero in rows k through n
           90 continue
           k = k + 1_${ik}$
           100 continue
           loop_150: do j = k, l
              do i = k, lm1
                 ip1 = i + 1_${ik}$
                 if( a( i, j )/=zero .or. b( i, j )/=zero )go to 120
              end do
              i = l
              go to 140
              120 continue
              do i = ip1, l
                 if( a( i, j )/=zero .or. b( i, j )/=zero )cycle loop_150
              end do
              i = ip1 - 1_${ik}$
              140 continue
              m = k
              iflow = 2_${ik}$
              go to 160
           end do loop_150
           go to 190
           ! permute rows m and i
           160 continue
           lscale( m ) = i
           if( i==m )go to 170
           call stdlib${ii}$_${ri}$swap( n-k+1, a( i, k ), lda, a( m, k ), lda )
           call stdlib${ii}$_${ri}$swap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
           ! permute columns m and j
           170 continue
           rscale( m ) = j
           if( j==m )go to 180
           call stdlib${ii}$_${ri}$swap( l, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, m ), 1_${ik}$ )
           call stdlib${ii}$_${ri}$swap( l, b( 1_${ik}$, j ), 1_${ik}$, b( 1_${ik}$, m ), 1_${ik}$ )
           180 continue
           go to ( 20, 90 )iflow
           190 continue
           ilo = k
           ihi = l
           if( stdlib_lsame( job, 'P' ) ) then
              do i = ilo, ihi
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           if( ilo==ihi )return
           ! balance the submatrix in rows ilo to ihi.
           nr = ihi - ilo + 1_${ik}$
           do i = ilo, ihi
              rscale( i ) = zero
              lscale( i ) = zero
              work( i ) = zero
              work( i+n ) = zero
              work( i+2*n ) = zero
              work( i+3*n ) = zero
              work( i+4*n ) = zero
              work( i+5*n ) = zero
           end do
           ! compute right side vector in resulting linear equations
           basl = log10( sclfac )
           do i = ilo, ihi
              do j = ilo, ihi
                 tb = b( i, j )
                 ta = a( i, j )
                 if( ta==zero )go to 210
                 ta = log10( abs( ta ) ) / basl
                 210 continue
                 if( tb==zero )go to 220
                 tb = log10( abs( tb ) ) / basl
                 220 continue
                 work( i+4*n ) = work( i+4*n ) - ta - tb
                 work( j+5*n ) = work( j+5*n ) - ta - tb
              end do
           end do
           coef = one / real( 2_${ik}$*nr,KIND=${rk}$)
           coef2 = coef*coef
           coef5 = half*coef2
           nrp2 = nr + 2_${ik}$
           beta = zero
           it = 1_${ik}$
           ! start generalized conjugate gradient iteration
           250 continue
           gamma = stdlib${ii}$_${ri}$dot( nr, work( ilo+4*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ ) +stdlib${ii}$_${ri}$dot( nr, &
                     work( ilo+5*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           ew = zero
           ewc = zero
           do i = ilo, ihi
              ew = ew + work( i+4*n )
              ewc = ewc + work( i+5*n )
           end do
           gamma = coef*gamma - coef2*( ew**2_${ik}$+ewc**2_${ik}$ ) - coef5*( ew-ewc )**2_${ik}$
           if( gamma==zero )go to 350
           if( it/=1_${ik}$ )beta = gamma / pgamma
           t = coef5*( ewc-three*ew )
           tc = coef5*( ew-three*ewc )
           call stdlib${ii}$_${ri}$scal( nr, beta, work( ilo ), 1_${ik}$ )
           call stdlib${ii}$_${ri}$scal( nr, beta, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_${ri}$axpy( nr, coef, work( ilo+4*n ), 1_${ik}$, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_${ri}$axpy( nr, coef, work( ilo+5*n ), 1_${ik}$, work( ilo ), 1_${ik}$ )
           do i = ilo, ihi
              work( i ) = work( i ) + tc
              work( i+n ) = work( i+n ) + t
           end do
           ! apply matrix to vector
           do i = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_290: do j = ilo, ihi
                 if( a( i, j )==zero )go to 280
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
                 280 continue
                 if( b( i, j )==zero )cycle loop_290
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
              end do loop_290
              work( i+2*n ) = real( kount,KIND=${rk}$)*work( i+n ) + sum
           end do
           do j = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_320: do i = ilo, ihi
                 if( a( i, j )==zero )go to 310
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
                 310 continue
                 if( b( i, j )==zero )cycle loop_320
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
              end do loop_320
              work( j+3*n ) = real( kount,KIND=${rk}$)*work( j ) + sum
           end do
           sum = stdlib${ii}$_${ri}$dot( nr, work( ilo+n ), 1_${ik}$, work( ilo+2*n ), 1_${ik}$ ) +stdlib${ii}$_${ri}$dot( nr, work( &
                     ilo ), 1_${ik}$, work( ilo+3*n ), 1_${ik}$ )
           alpha = gamma / sum
           ! determine correction to current iteration
           cmax = zero
           do i = ilo, ihi
              cor = alpha*work( i+n )
              if( abs( cor )>cmax )cmax = abs( cor )
              lscale( i ) = lscale( i ) + cor
              cor = alpha*work( i )
              if( abs( cor )>cmax )cmax = abs( cor )
              rscale( i ) = rscale( i ) + cor
           end do
           if( cmax<half )go to 350
           call stdlib${ii}$_${ri}$axpy( nr, -alpha, work( ilo+2*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ )
           call stdlib${ii}$_${ri}$axpy( nr, -alpha, work( ilo+3*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           pgamma = gamma
           it = it + 1_${ik}$
           if( it<=nrp2 )go to 250
           ! end generalized conjugate gradient iteration
           350 continue
           sfmin = stdlib${ii}$_${ri}$lamch( 'S' )
           sfmax = one / sfmin
           lsfmin = int( log10( sfmin ) / basl+one,KIND=${ik}$)
           lsfmax = int( log10( sfmax ) / basl,KIND=${ik}$)
           do i = ilo, ihi
              irab = stdlib${ii}$_i${ri}$amax( n-ilo+1, a( i, ilo ), lda )
              rab = abs( a( i, irab+ilo-1 ) )
              irab = stdlib${ii}$_i${ri}$amax( n-ilo+1, b( i, ilo ), ldb )
              rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
              lrab = int( log10( rab+sfmin ) / basl+one,KIND=${ik}$)
              ir = int(lscale( i ) + sign( half, lscale( i ) ),KIND=${ik}$)
              ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
              lscale( i ) = sclfac**ir
              icab = stdlib${ii}$_i${ri}$amax( ihi, a( 1_${ik}$, i ), 1_${ik}$ )
              cab = abs( a( icab, i ) )
              icab = stdlib${ii}$_i${ri}$amax( ihi, b( 1_${ik}$, i ), 1_${ik}$ )
              cab = max( cab, abs( b( icab, i ) ) )
              lcab = int( log10( cab+sfmin ) / basl+one,KIND=${ik}$)
              jc = int(rscale( i ) + sign( half, rscale( i ) ),KIND=${ik}$)
              jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
              rscale( i ) = sclfac**jc
           end do
           ! row scaling of matrices a and b
           do i = ilo, ihi
              call stdlib${ii}$_${ri}$scal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
              call stdlib${ii}$_${ri}$scal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
           end do
           ! column scaling of matrices a and b
           do j = ilo, ihi
              call stdlib${ii}$_${ri}$scal( ihi, rscale( j ), a( 1_${ik}$, j ), 1_${ik}$ )
              call stdlib${ii}$_${ri}$scal( ihi, rscale( j ), b( 1_${ik}$, j ), 1_${ik}$ )
           end do
           return
     end subroutine stdlib${ii}$_${ri}$ggbal

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_cggbal( job, n, a, lda, b, ldb, ilo, ihi, lscale,rscale, work, info )
     !! CGGBAL balances a pair of general complex matrices (A,B).  This
     !! involves, first, permuting A and B by similarity transformations to
     !! isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
     !! elements on the diagonal; and second, applying a diagonal similarity
     !! transformation to rows and columns ILO to IHI to make the rows
     !! and columns as close in norm as possible. Both steps are optional.
     !! Balancing may reduce the 1-norm of the matrices, and improve the
     !! accuracy of the computed eigenvalues and/or eigenvectors in the
     !! generalized eigenvalue problem A*x = lambda*B*x.
               
        ! -- 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) :: job
           integer(${ik}$), intent(out) :: ihi, ilo, info
           integer(${ik}$), intent(in) :: lda, ldb, n
           ! Array Arguments 
           real(sp), intent(out) :: lscale(*), rscale(*), work(*)
           complex(sp), intent(inout) :: a(lda,*), b(ldb,*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: sclfac = ten
           
           
           
           ! Local Scalars 
           integer(${ik}$) :: i, icab, iflow, ip1, ir, irab, it, j, jc, jp1, k, kount, l, lcab, lm1, &
                     lrab, lsfmax, lsfmin, m, nr, nrp2
           real(sp) :: alpha, basl, beta, cab, cmax, coef, coef2, coef5, cor, ew, ewc, gamma, &
                     pgamma, rab, sfmax, sfmin, sum, t, ta, tb, tc
           complex(sp) :: cdum
           ! Intrinsic Functions 
           ! Statement Functions 
           real(sp) :: cabs1
           ! Statement Function Definitions 
           cabs1( cdum ) = abs( real( cdum,KIND=sp) ) + abs( aimag( cdum ) )
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'N' ) .and. .not.stdlib_lsame( job, 'P' ) &
                     .and..not.stdlib_lsame( job, 'S' ) .and. .not.stdlib_lsame( job, 'B' ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CGGBAL', -info )
              return
           end if
           ! quick return if possible
           if( n==0_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              return
           end if
           if( n==1_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              lscale( 1_${ik}$ ) = one
              rscale( 1_${ik}$ ) = one
              return
           end if
           if( stdlib_lsame( job, 'N' ) ) then
              ilo = 1_${ik}$
              ihi = n
              do i = 1, n
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           k = 1_${ik}$
           l = n
           if( stdlib_lsame( job, 'S' ) )go to 190
           go to 30
           ! permute the matrices a and b to isolate the eigenvalues.
           ! find row with one nonzero in columns 1 through l
           20 continue
           l = lm1
           if( l/=1 )go to 30
           rscale( 1_${ik}$ ) = one
           lscale( 1_${ik}$ ) = one
           go to 190
           30 continue
           lm1 = l - 1_${ik}$
           loop_80: do i = l, 1, -1
              do j = 1, lm1
                 jp1 = j + 1_${ik}$
                 if( a( i, j )/=czero .or. b( i, j )/=czero )go to 50
              end do
              j = l
              go to 70
              50 continue
              do j = jp1, l
                 if( a( i, j )/=czero .or. b( i, j )/=czero )cycle loop_80
              end do
              j = jp1 - 1_${ik}$
              70 continue
              m = l
              iflow = 1_${ik}$
              go to 160
           end do loop_80
           go to 100
           ! find column with one nonzero in rows k through n
           90 continue
           k = k + 1_${ik}$
           100 continue
           loop_150: do j = k, l
              do i = k, lm1
                 ip1 = i + 1_${ik}$
                 if( a( i, j )/=czero .or. b( i, j )/=czero )go to 120
              end do
              i = l
              go to 140
              120 continue
              do i = ip1, l
                 if( a( i, j )/=czero .or. b( i, j )/=czero )cycle loop_150
              end do
              i = ip1 - 1_${ik}$
              140 continue
              m = k
              iflow = 2_${ik}$
              go to 160
           end do loop_150
           go to 190
           ! permute rows m and i
           160 continue
           lscale( m ) = i
           if( i==m )go to 170
           call stdlib${ii}$_cswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
           call stdlib${ii}$_cswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
           ! permute columns m and j
           170 continue
           rscale( m ) = j
           if( j==m )go to 180
           call stdlib${ii}$_cswap( l, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, m ), 1_${ik}$ )
           call stdlib${ii}$_cswap( l, b( 1_${ik}$, j ), 1_${ik}$, b( 1_${ik}$, m ), 1_${ik}$ )
           180 continue
           go to ( 20, 90 )iflow
           190 continue
           ilo = k
           ihi = l
           if( stdlib_lsame( job, 'P' ) ) then
              do i = ilo, ihi
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           if( ilo==ihi )return
           ! balance the submatrix in rows ilo to ihi.
           nr = ihi - ilo + 1_${ik}$
           do i = ilo, ihi
              rscale( i ) = zero
              lscale( i ) = zero
              work( i ) = zero
              work( i+n ) = zero
              work( i+2*n ) = zero
              work( i+3*n ) = zero
              work( i+4*n ) = zero
              work( i+5*n ) = zero
           end do
           ! compute right side vector in resulting linear equations
           basl = log10( sclfac )
           do i = ilo, ihi
              do j = ilo, ihi
                 if( a( i, j )==czero ) then
                    ta = zero
                    go to 210
                 end if
                 ta = log10( cabs1( a( i, j ) ) ) / basl
                 210 continue
                 if( b( i, j )==czero ) then
                    tb = zero
                    go to 220
                 end if
                 tb = log10( cabs1( b( i, j ) ) ) / basl
                 220 continue
                 work( i+4*n ) = work( i+4*n ) - ta - tb
                 work( j+5*n ) = work( j+5*n ) - ta - tb
              end do
           end do
           coef = one / real( 2_${ik}$*nr,KIND=sp)
           coef2 = coef*coef
           coef5 = half*coef2
           nrp2 = nr + 2_${ik}$
           beta = zero
           it = 1_${ik}$
           ! start generalized conjugate gradient iteration
           250 continue
           gamma = stdlib${ii}$_sdot( nr, work( ilo+4*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ ) +stdlib${ii}$_sdot( nr, &
                     work( ilo+5*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           ew = zero
           ewc = zero
           do i = ilo, ihi
              ew = ew + work( i+4*n )
              ewc = ewc + work( i+5*n )
           end do
           gamma = coef*gamma - coef2*( ew**2_${ik}$+ewc**2_${ik}$ ) - coef5*( ew-ewc )**2_${ik}$
           if( gamma==zero )go to 350
           if( it/=1_${ik}$ )beta = gamma / pgamma
           t = coef5*( ewc-three*ew )
           tc = coef5*( ew-three*ewc )
           call stdlib${ii}$_sscal( nr, beta, work( ilo ), 1_${ik}$ )
           call stdlib${ii}$_sscal( nr, beta, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_saxpy( nr, coef, work( ilo+4*n ), 1_${ik}$, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_saxpy( nr, coef, work( ilo+5*n ), 1_${ik}$, work( ilo ), 1_${ik}$ )
           do i = ilo, ihi
              work( i ) = work( i ) + tc
              work( i+n ) = work( i+n ) + t
           end do
           ! apply matrix to vector
           do i = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_290: do j = ilo, ihi
                 if( a( i, j )==czero )go to 280
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
                 280 continue
                 if( b( i, j )==czero )cycle loop_290
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
              end do loop_290
              work( i+2*n ) = real( kount,KIND=sp)*work( i+n ) + sum
           end do
           do j = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_320: do i = ilo, ihi
                 if( a( i, j )==czero )go to 310
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
                 310 continue
                 if( b( i, j )==czero )cycle loop_320
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
              end do loop_320
              work( j+3*n ) = real( kount,KIND=sp)*work( j ) + sum
           end do
           sum = stdlib${ii}$_sdot( nr, work( ilo+n ), 1_${ik}$, work( ilo+2*n ), 1_${ik}$ ) +stdlib${ii}$_sdot( nr, work( &
                     ilo ), 1_${ik}$, work( ilo+3*n ), 1_${ik}$ )
           alpha = gamma / sum
           ! determine correction to current iteration
           cmax = zero
           do i = ilo, ihi
              cor = alpha*work( i+n )
              if( abs( cor )>cmax )cmax = abs( cor )
              lscale( i ) = lscale( i ) + cor
              cor = alpha*work( i )
              if( abs( cor )>cmax )cmax = abs( cor )
              rscale( i ) = rscale( i ) + cor
           end do
           if( cmax<half )go to 350
           call stdlib${ii}$_saxpy( nr, -alpha, work( ilo+2*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ )
           call stdlib${ii}$_saxpy( nr, -alpha, work( ilo+3*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           pgamma = gamma
           it = it + 1_${ik}$
           if( it<=nrp2 )go to 250
           ! end generalized conjugate gradient iteration
           350 continue
           sfmin = stdlib${ii}$_slamch( 'S' )
           sfmax = one / sfmin
           lsfmin = int( log10( sfmin ) / basl+one,KIND=${ik}$)
           lsfmax = int( log10( sfmax ) / basl,KIND=${ik}$)
           do i = ilo, ihi
              irab = stdlib${ii}$_icamax( n-ilo+1, a( i, ilo ), lda )
              rab = abs( a( i, irab+ilo-1 ) )
              irab = stdlib${ii}$_icamax( n-ilo+1, b( i, ilo ), ldb )
              rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
              lrab = int( log10( rab+sfmin ) / basl+one,KIND=${ik}$)
              ir = lscale( i ) + sign( half, lscale( i ) )
              ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
              lscale( i ) = sclfac**ir
              icab = stdlib${ii}$_icamax( ihi, a( 1_${ik}$, i ), 1_${ik}$ )
              cab = abs( a( icab, i ) )
              icab = stdlib${ii}$_icamax( ihi, b( 1_${ik}$, i ), 1_${ik}$ )
              cab = max( cab, abs( b( icab, i ) ) )
              lcab = int( log10( cab+sfmin ) / basl+one,KIND=${ik}$)
              jc = rscale( i ) + sign( half, rscale( i ) )
              jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
              rscale( i ) = sclfac**jc
           end do
           ! row scaling of matrices a and b
           do i = ilo, ihi
              call stdlib${ii}$_csscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
              call stdlib${ii}$_csscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
           end do
           ! column scaling of matrices a and b
           do j = ilo, ihi
              call stdlib${ii}$_csscal( ihi, rscale( j ), a( 1_${ik}$, j ), 1_${ik}$ )
              call stdlib${ii}$_csscal( ihi, rscale( j ), b( 1_${ik}$, j ), 1_${ik}$ )
           end do
           return
     end subroutine stdlib${ii}$_cggbal

     pure module subroutine stdlib${ii}$_zggbal( job, n, a, lda, b, ldb, ilo, ihi, lscale,rscale, work, info )
     !! ZGGBAL balances a pair of general complex matrices (A,B).  This
     !! involves, first, permuting A and B by similarity transformations to
     !! isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
     !! elements on the diagonal; and second, applying a diagonal similarity
     !! transformation to rows and columns ILO to IHI to make the rows
     !! and columns as close in norm as possible. Both steps are optional.
     !! Balancing may reduce the 1-norm of the matrices, and improve the
     !! accuracy of the computed eigenvalues and/or eigenvectors in the
     !! generalized eigenvalue problem A*x = lambda*B*x.
               
        ! -- 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) :: job
           integer(${ik}$), intent(out) :: ihi, ilo, info
           integer(${ik}$), intent(in) :: lda, ldb, n
           ! Array Arguments 
           real(dp), intent(out) :: lscale(*), rscale(*), work(*)
           complex(dp), intent(inout) :: a(lda,*), b(ldb,*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: sclfac = ten
           
           
           
           ! Local Scalars 
           integer(${ik}$) :: i, icab, iflow, ip1, ir, irab, it, j, jc, jp1, k, kount, l, lcab, lm1, &
                     lrab, lsfmax, lsfmin, m, nr, nrp2
           real(dp) :: alpha, basl, beta, cab, cmax, coef, coef2, coef5, cor, ew, ewc, gamma, &
                     pgamma, rab, sfmax, sfmin, sum, t, ta, tb, tc
           complex(dp) :: cdum
           ! Intrinsic Functions 
           ! Statement Functions 
           real(dp) :: cabs1
           ! Statement Function Definitions 
           cabs1( cdum ) = abs( real( cdum,KIND=dp) ) + abs( aimag( cdum ) )
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'N' ) .and. .not.stdlib_lsame( job, 'P' ) &
                     .and..not.stdlib_lsame( job, 'S' ) .and. .not.stdlib_lsame( job, 'B' ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGGBAL', -info )
              return
           end if
           ! quick return if possible
           if( n==0_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              return
           end if
           if( n==1_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              lscale( 1_${ik}$ ) = one
              rscale( 1_${ik}$ ) = one
              return
           end if
           if( stdlib_lsame( job, 'N' ) ) then
              ilo = 1_${ik}$
              ihi = n
              do i = 1, n
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           k = 1_${ik}$
           l = n
           if( stdlib_lsame( job, 'S' ) )go to 190
           go to 30
           ! permute the matrices a and b to isolate the eigenvalues.
           ! find row with one nonzero in columns 1 through l
           20 continue
           l = lm1
           if( l/=1 )go to 30
           rscale( 1_${ik}$ ) = 1_${ik}$
           lscale( 1_${ik}$ ) = 1_${ik}$
           go to 190
           30 continue
           lm1 = l - 1_${ik}$
           loop_80: do i = l, 1, -1
              do j = 1, lm1
                 jp1 = j + 1_${ik}$
                 if( a( i, j )/=czero .or. b( i, j )/=czero )go to 50
              end do
              j = l
              go to 70
              50 continue
              do j = jp1, l
                 if( a( i, j )/=czero .or. b( i, j )/=czero )cycle loop_80
              end do
              j = jp1 - 1_${ik}$
              70 continue
              m = l
              iflow = 1_${ik}$
              go to 160
           end do loop_80
           go to 100
           ! find column with one nonzero in rows k through n
           90 continue
           k = k + 1_${ik}$
           100 continue
           loop_150: do j = k, l
              do i = k, lm1
                 ip1 = i + 1_${ik}$
                 if( a( i, j )/=czero .or. b( i, j )/=czero )go to 120
              end do
              i = l
              go to 140
              120 continue
              do i = ip1, l
                 if( a( i, j )/=czero .or. b( i, j )/=czero )cycle loop_150
              end do
              i = ip1 - 1_${ik}$
              140 continue
              m = k
              iflow = 2_${ik}$
              go to 160
           end do loop_150
           go to 190
           ! permute rows m and i
           160 continue
           lscale( m ) = i
           if( i==m )go to 170
           call stdlib${ii}$_zswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
           call stdlib${ii}$_zswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
           ! permute columns m and j
           170 continue
           rscale( m ) = j
           if( j==m )go to 180
           call stdlib${ii}$_zswap( l, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, m ), 1_${ik}$ )
           call stdlib${ii}$_zswap( l, b( 1_${ik}$, j ), 1_${ik}$, b( 1_${ik}$, m ), 1_${ik}$ )
           180 continue
           go to ( 20, 90 )iflow
           190 continue
           ilo = k
           ihi = l
           if( stdlib_lsame( job, 'P' ) ) then
              do i = ilo, ihi
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           if( ilo==ihi )return
           ! balance the submatrix in rows ilo to ihi.
           nr = ihi - ilo + 1_${ik}$
           do i = ilo, ihi
              rscale( i ) = zero
              lscale( i ) = zero
              work( i ) = zero
              work( i+n ) = zero
              work( i+2*n ) = zero
              work( i+3*n ) = zero
              work( i+4*n ) = zero
              work( i+5*n ) = zero
           end do
           ! compute right side vector in resulting linear equations
           basl = log10( sclfac )
           do i = ilo, ihi
              do j = ilo, ihi
                 if( a( i, j )==czero ) then
                    ta = zero
                    go to 210
                 end if
                 ta = log10( cabs1( a( i, j ) ) ) / basl
                 210 continue
                 if( b( i, j )==czero ) then
                    tb = zero
                    go to 220
                 end if
                 tb = log10( cabs1( b( i, j ) ) ) / basl
                 220 continue
                 work( i+4*n ) = work( i+4*n ) - ta - tb
                 work( j+5*n ) = work( j+5*n ) - ta - tb
              end do
           end do
           coef = one / real( 2_${ik}$*nr,KIND=dp)
           coef2 = coef*coef
           coef5 = half*coef2
           nrp2 = nr + 2_${ik}$
           beta = zero
           it = 1_${ik}$
           ! start generalized conjugate gradient iteration
           250 continue
           gamma = stdlib${ii}$_ddot( nr, work( ilo+4*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ ) +stdlib${ii}$_ddot( nr, &
                     work( ilo+5*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           ew = zero
           ewc = zero
           do i = ilo, ihi
              ew = ew + work( i+4*n )
              ewc = ewc + work( i+5*n )
           end do
           gamma = coef*gamma - coef2*( ew**2_${ik}$+ewc**2_${ik}$ ) - coef5*( ew-ewc )**2_${ik}$
           if( gamma==zero )go to 350
           if( it/=1_${ik}$ )beta = gamma / pgamma
           t = coef5*( ewc-three*ew )
           tc = coef5*( ew-three*ewc )
           call stdlib${ii}$_dscal( nr, beta, work( ilo ), 1_${ik}$ )
           call stdlib${ii}$_dscal( nr, beta, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_daxpy( nr, coef, work( ilo+4*n ), 1_${ik}$, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_daxpy( nr, coef, work( ilo+5*n ), 1_${ik}$, work( ilo ), 1_${ik}$ )
           do i = ilo, ihi
              work( i ) = work( i ) + tc
              work( i+n ) = work( i+n ) + t
           end do
           ! apply matrix to vector
           do i = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_290: do j = ilo, ihi
                 if( a( i, j )==czero )go to 280
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
                 280 continue
                 if( b( i, j )==czero )cycle loop_290
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
              end do loop_290
              work( i+2*n ) = real( kount,KIND=dp)*work( i+n ) + sum
           end do
           do j = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_320: do i = ilo, ihi
                 if( a( i, j )==czero )go to 310
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
                 310 continue
                 if( b( i, j )==czero )cycle loop_320
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
              end do loop_320
              work( j+3*n ) = real( kount,KIND=dp)*work( j ) + sum
           end do
           sum = stdlib${ii}$_ddot( nr, work( ilo+n ), 1_${ik}$, work( ilo+2*n ), 1_${ik}$ ) +stdlib${ii}$_ddot( nr, work( &
                     ilo ), 1_${ik}$, work( ilo+3*n ), 1_${ik}$ )
           alpha = gamma / sum
           ! determine correction to current iteration
           cmax = zero
           do i = ilo, ihi
              cor = alpha*work( i+n )
              if( abs( cor )>cmax )cmax = abs( cor )
              lscale( i ) = lscale( i ) + cor
              cor = alpha*work( i )
              if( abs( cor )>cmax )cmax = abs( cor )
              rscale( i ) = rscale( i ) + cor
           end do
           if( cmax<half )go to 350
           call stdlib${ii}$_daxpy( nr, -alpha, work( ilo+2*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ )
           call stdlib${ii}$_daxpy( nr, -alpha, work( ilo+3*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           pgamma = gamma
           it = it + 1_${ik}$
           if( it<=nrp2 )go to 250
           ! end generalized conjugate gradient iteration
           350 continue
           sfmin = stdlib${ii}$_dlamch( 'S' )
           sfmax = one / sfmin
           lsfmin = int( log10( sfmin ) / basl+one,KIND=${ik}$)
           lsfmax = int( log10( sfmax ) / basl,KIND=${ik}$)
           do i = ilo, ihi
              irab = stdlib${ii}$_izamax( n-ilo+1, a( i, ilo ), lda )
              rab = abs( a( i, irab+ilo-1 ) )
              irab = stdlib${ii}$_izamax( n-ilo+1, b( i, ilo ), ldb )
              rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
              lrab = int( log10( rab+sfmin ) / basl+one,KIND=${ik}$)
              ir = int(lscale( i ) + sign( half, lscale( i ) ),KIND=${ik}$)
              ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
              lscale( i ) = sclfac**ir
              icab = stdlib${ii}$_izamax( ihi, a( 1_${ik}$, i ), 1_${ik}$ )
              cab = abs( a( icab, i ) )
              icab = stdlib${ii}$_izamax( ihi, b( 1_${ik}$, i ), 1_${ik}$ )
              cab = max( cab, abs( b( icab, i ) ) )
              lcab = int( log10( cab+sfmin ) / basl+one,KIND=${ik}$)
              jc = int(rscale( i ) + sign( half, rscale( i ) ),KIND=${ik}$)
              jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
              rscale( i ) = sclfac**jc
           end do
           ! row scaling of matrices a and b
           do i = ilo, ihi
              call stdlib${ii}$_zdscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
              call stdlib${ii}$_zdscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
           end do
           ! column scaling of matrices a and b
           do j = ilo, ihi
              call stdlib${ii}$_zdscal( ihi, rscale( j ), a( 1_${ik}$, j ), 1_${ik}$ )
              call stdlib${ii}$_zdscal( ihi, rscale( j ), b( 1_${ik}$, j ), 1_${ik}$ )
           end do
           return
     end subroutine stdlib${ii}$_zggbal

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$ggbal( job, n, a, lda, b, ldb, ilo, ihi, lscale,rscale, work, info )
     !! ZGGBAL: balances a pair of general complex matrices (A,B).  This
     !! involves, first, permuting A and B by similarity transformations to
     !! isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
     !! elements on the diagonal; and second, applying a diagonal similarity
     !! transformation to rows and columns ILO to IHI to make the rows
     !! and columns as close in norm as possible. Both steps are optional.
     !! Balancing may reduce the 1-norm of the matrices, and improve the
     !! accuracy of the computed eigenvalues and/or eigenvectors in the
     !! generalized eigenvalue problem A*x = lambda*B*x.
               
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: job
           integer(${ik}$), intent(out) :: ihi, ilo, info
           integer(${ik}$), intent(in) :: lda, ldb, n
           ! Array Arguments 
           real(${ck}$), intent(out) :: lscale(*), rscale(*), work(*)
           complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*)
        ! =====================================================================
           ! Parameters 
           real(${ck}$), parameter :: sclfac = ten
           
           
           
           ! Local Scalars 
           integer(${ik}$) :: i, icab, iflow, ip1, ir, irab, it, j, jc, jp1, k, kount, l, lcab, lm1, &
                     lrab, lsfmax, lsfmin, m, nr, nrp2
           real(${ck}$) :: alpha, basl, beta, cab, cmax, coef, coef2, coef5, cor, ew, ewc, gamma, &
                     pgamma, rab, sfmax, sfmin, sum, t, ta, tb, tc
           complex(${ck}$) :: cdum
           ! Intrinsic Functions 
           ! Statement Functions 
           real(${ck}$) :: cabs1
           ! Statement Function Definitions 
           cabs1( cdum ) = abs( real( cdum,KIND=${ck}$) ) + abs( aimag( cdum ) )
           ! Executable Statements 
           ! test the input parameters
           info = 0_${ik}$
           if( .not.stdlib_lsame( job, 'N' ) .and. .not.stdlib_lsame( job, 'P' ) &
                     .and..not.stdlib_lsame( job, 'S' ) .and. .not.stdlib_lsame( job, 'B' ) ) then
              info = -1_${ik}$
           else if( n<0_${ik}$ ) then
              info = -2_${ik}$
           else if( lda<max( 1_${ik}$, n ) ) then
              info = -4_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGGBAL', -info )
              return
           end if
           ! quick return if possible
           if( n==0_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              return
           end if
           if( n==1_${ik}$ ) then
              ilo = 1_${ik}$
              ihi = n
              lscale( 1_${ik}$ ) = one
              rscale( 1_${ik}$ ) = one
              return
           end if
           if( stdlib_lsame( job, 'N' ) ) then
              ilo = 1_${ik}$
              ihi = n
              do i = 1, n
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           k = 1_${ik}$
           l = n
           if( stdlib_lsame( job, 'S' ) )go to 190
           go to 30
           ! permute the matrices a and b to isolate the eigenvalues.
           ! find row with one nonzero in columns 1 through l
           20 continue
           l = lm1
           if( l/=1 )go to 30
           rscale( 1_${ik}$ ) = 1_${ik}$
           lscale( 1_${ik}$ ) = 1_${ik}$
           go to 190
           30 continue
           lm1 = l - 1_${ik}$
           loop_80: do i = l, 1, -1
              do j = 1, lm1
                 jp1 = j + 1_${ik}$
                 if( a( i, j )/=czero .or. b( i, j )/=czero )go to 50
              end do
              j = l
              go to 70
              50 continue
              do j = jp1, l
                 if( a( i, j )/=czero .or. b( i, j )/=czero )cycle loop_80
              end do
              j = jp1 - 1_${ik}$
              70 continue
              m = l
              iflow = 1_${ik}$
              go to 160
           end do loop_80
           go to 100
           ! find column with one nonzero in rows k through n
           90 continue
           k = k + 1_${ik}$
           100 continue
           loop_150: do j = k, l
              do i = k, lm1
                 ip1 = i + 1_${ik}$
                 if( a( i, j )/=czero .or. b( i, j )/=czero )go to 120
              end do
              i = l
              go to 140
              120 continue
              do i = ip1, l
                 if( a( i, j )/=czero .or. b( i, j )/=czero )cycle loop_150
              end do
              i = ip1 - 1_${ik}$
              140 continue
              m = k
              iflow = 2_${ik}$
              go to 160
           end do loop_150
           go to 190
           ! permute rows m and i
           160 continue
           lscale( m ) = i
           if( i==m )go to 170
           call stdlib${ii}$_${ci}$swap( n-k+1, a( i, k ), lda, a( m, k ), lda )
           call stdlib${ii}$_${ci}$swap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
           ! permute columns m and j
           170 continue
           rscale( m ) = j
           if( j==m )go to 180
           call stdlib${ii}$_${ci}$swap( l, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, m ), 1_${ik}$ )
           call stdlib${ii}$_${ci}$swap( l, b( 1_${ik}$, j ), 1_${ik}$, b( 1_${ik}$, m ), 1_${ik}$ )
           180 continue
           go to ( 20, 90 )iflow
           190 continue
           ilo = k
           ihi = l
           if( stdlib_lsame( job, 'P' ) ) then
              do i = ilo, ihi
                 lscale( i ) = one
                 rscale( i ) = one
              end do
              return
           end if
           if( ilo==ihi )return
           ! balance the submatrix in rows ilo to ihi.
           nr = ihi - ilo + 1_${ik}$
           do i = ilo, ihi
              rscale( i ) = zero
              lscale( i ) = zero
              work( i ) = zero
              work( i+n ) = zero
              work( i+2*n ) = zero
              work( i+3*n ) = zero
              work( i+4*n ) = zero
              work( i+5*n ) = zero
           end do
           ! compute right side vector in resulting linear equations
           basl = log10( sclfac )
           do i = ilo, ihi
              do j = ilo, ihi
                 if( a( i, j )==czero ) then
                    ta = zero
                    go to 210
                 end if
                 ta = log10( cabs1( a( i, j ) ) ) / basl
                 210 continue
                 if( b( i, j )==czero ) then
                    tb = zero
                    go to 220
                 end if
                 tb = log10( cabs1( b( i, j ) ) ) / basl
                 220 continue
                 work( i+4*n ) = work( i+4*n ) - ta - tb
                 work( j+5*n ) = work( j+5*n ) - ta - tb
              end do
           end do
           coef = one / real( 2_${ik}$*nr,KIND=${ck}$)
           coef2 = coef*coef
           coef5 = half*coef2
           nrp2 = nr + 2_${ik}$
           beta = zero
           it = 1_${ik}$
           ! start generalized conjugate gradient iteration
           250 continue
           gamma = stdlib${ii}$_${c2ri(ci)}$dot( nr, work( ilo+4*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ ) +stdlib${ii}$_${c2ri(ci)}$dot( nr, &
                     work( ilo+5*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           ew = zero
           ewc = zero
           do i = ilo, ihi
              ew = ew + work( i+4*n )
              ewc = ewc + work( i+5*n )
           end do
           gamma = coef*gamma - coef2*( ew**2_${ik}$+ewc**2_${ik}$ ) - coef5*( ew-ewc )**2_${ik}$
           if( gamma==zero )go to 350
           if( it/=1_${ik}$ )beta = gamma / pgamma
           t = coef5*( ewc-three*ew )
           tc = coef5*( ew-three*ewc )
           call stdlib${ii}$_${c2ri(ci)}$scal( nr, beta, work( ilo ), 1_${ik}$ )
           call stdlib${ii}$_${c2ri(ci)}$scal( nr, beta, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_${c2ri(ci)}$axpy( nr, coef, work( ilo+4*n ), 1_${ik}$, work( ilo+n ), 1_${ik}$ )
           call stdlib${ii}$_${c2ri(ci)}$axpy( nr, coef, work( ilo+5*n ), 1_${ik}$, work( ilo ), 1_${ik}$ )
           do i = ilo, ihi
              work( i ) = work( i ) + tc
              work( i+n ) = work( i+n ) + t
           end do
           ! apply matrix to vector
           do i = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_290: do j = ilo, ihi
                 if( a( i, j )==czero )go to 280
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
                 280 continue
                 if( b( i, j )==czero )cycle loop_290
                 kount = kount + 1_${ik}$
                 sum = sum + work( j )
              end do loop_290
              work( i+2*n ) = real( kount,KIND=${ck}$)*work( i+n ) + sum
           end do
           do j = ilo, ihi
              kount = 0_${ik}$
              sum = zero
              loop_320: do i = ilo, ihi
                 if( a( i, j )==czero )go to 310
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
                 310 continue
                 if( b( i, j )==czero )cycle loop_320
                 kount = kount + 1_${ik}$
                 sum = sum + work( i+n )
              end do loop_320
              work( j+3*n ) = real( kount,KIND=${ck}$)*work( j ) + sum
           end do
           sum = stdlib${ii}$_${c2ri(ci)}$dot( nr, work( ilo+n ), 1_${ik}$, work( ilo+2*n ), 1_${ik}$ ) +stdlib${ii}$_${c2ri(ci)}$dot( nr, work( &
                     ilo ), 1_${ik}$, work( ilo+3*n ), 1_${ik}$ )
           alpha = gamma / sum
           ! determine correction to current iteration
           cmax = zero
           do i = ilo, ihi
              cor = alpha*work( i+n )
              if( abs( cor )>cmax )cmax = abs( cor )
              lscale( i ) = lscale( i ) + cor
              cor = alpha*work( i )
              if( abs( cor )>cmax )cmax = abs( cor )
              rscale( i ) = rscale( i ) + cor
           end do
           if( cmax<half )go to 350
           call stdlib${ii}$_${c2ri(ci)}$axpy( nr, -alpha, work( ilo+2*n ), 1_${ik}$, work( ilo+4*n ), 1_${ik}$ )
           call stdlib${ii}$_${c2ri(ci)}$axpy( nr, -alpha, work( ilo+3*n ), 1_${ik}$, work( ilo+5*n ), 1_${ik}$ )
           pgamma = gamma
           it = it + 1_${ik}$
           if( it<=nrp2 )go to 250
           ! end generalized conjugate gradient iteration
           350 continue
           sfmin = stdlib${ii}$_${c2ri(ci)}$lamch( 'S' )
           sfmax = one / sfmin
           lsfmin = int( log10( sfmin ) / basl+one,KIND=${ik}$)
           lsfmax = int( log10( sfmax ) / basl,KIND=${ik}$)
           do i = ilo, ihi
              irab = stdlib${ii}$_i${ci}$amax( n-ilo+1, a( i, ilo ), lda )
              rab = abs( a( i, irab+ilo-1 ) )
              irab = stdlib${ii}$_i${ci}$amax( n-ilo+1, b( i, ilo ), ldb )
              rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
              lrab = int( log10( rab+sfmin ) / basl+one,KIND=${ik}$)
              ir = int(lscale( i ) + sign( half, lscale( i ) ),KIND=${ik}$)
              ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
              lscale( i ) = sclfac**ir
              icab = stdlib${ii}$_i${ci}$amax( ihi, a( 1_${ik}$, i ), 1_${ik}$ )
              cab = abs( a( icab, i ) )
              icab = stdlib${ii}$_i${ci}$amax( ihi, b( 1_${ik}$, i ), 1_${ik}$ )
              cab = max( cab, abs( b( icab, i ) ) )
              lcab = int( log10( cab+sfmin ) / basl+one,KIND=${ik}$)
              jc = int(rscale( i ) + sign( half, rscale( i ) ),KIND=${ik}$)
              jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
              rscale( i ) = sclfac**jc
           end do
           ! row scaling of matrices a and b
           do i = ilo, ihi
              call stdlib${ii}$_${ci}$dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
              call stdlib${ii}$_${ci}$dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
           end do
           ! column scaling of matrices a and b
           do j = ilo, ihi
              call stdlib${ii}$_${ci}$dscal( ihi, rscale( j ), a( 1_${ik}$, j ), 1_${ik}$ )
              call stdlib${ii}$_${ci}$dscal( ihi, rscale( j ), b( 1_${ik}$, j ), 1_${ik}$ )
           end do
           return
     end subroutine stdlib${ii}$_${ci}$ggbal

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sgghrd( compq, compz, n, ilo, ihi, a, lda, b, ldb, q,ldq, z, ldz, &
     !! SGGHRD reduces a pair of real matrices (A,B) to generalized upper
     !! Hessenberg form using orthogonal transformations, where A is a
     !! general matrix and B is upper triangular.  The form of the
     !! generalized eigenvalue problem is
     !! A*x = lambda*B*x,
     !! and B is typically made upper triangular by computing its QR
     !! factorization and moving the orthogonal matrix Q to the left side
     !! of the equation.
     !! This subroutine simultaneously reduces A to a Hessenberg matrix H:
     !! Q**T*A*Z = H
     !! and transforms B to another upper triangular matrix T:
     !! Q**T*B*Z = T
     !! in order to reduce the problem to its standard form
     !! H*y = lambda*T*y
     !! where y = Z**T*x.
     !! The orthogonal matrices Q and Z are determined as products of Givens
     !! rotations.  They may either be formed explicitly, or they may be
     !! postmultiplied into input matrices Q1 and Z1, so that
     !! Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T
     !! Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T
     !! If Q1 is the orthogonal matrix from the QR factorization of B in the
     !! original equation A*x = lambda*B*x, then SGGHRD reduces the original
     !! problem to generalized Hessenberg form.
               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) :: compq, compz
           integer(${ik}$), intent(in) :: ihi, ilo, lda, ldb, ldq, ldz, n
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: ilq, ilz
           integer(${ik}$) :: icompq, icompz, jcol, jrow
           real(sp) :: c, s, temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode compq
           if( stdlib_lsame( compq, 'N' ) ) then
              ilq = .false.
              icompq = 1_${ik}$
           else if( stdlib_lsame( compq, 'V' ) ) then
              ilq = .true.
              icompq = 2_${ik}$
           else if( stdlib_lsame( compq, 'I' ) ) then
              ilq = .true.
              icompq = 3_${ik}$
           else
              icompq = 0_${ik}$
           end if
           ! decode compz
           if( stdlib_lsame( compz, 'N' ) ) then
              ilz = .false.
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              ilz = .true.
              icompz = 2_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              ilz = .true.
              icompz = 3_${ik}$
           else
              icompz = 0_${ik}$
           end if
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( icompz<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ ) then
              info = -4_${ik}$
           else if( ihi>n .or. ihi<ilo-1 ) 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( ( ilq .and. ldq<n ) .or. ldq<1_${ik}$ ) then
              info = -11_${ik}$
           else if( ( ilz .and. ldz<n ) .or. ldz<1_${ik}$ ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SGGHRD', -info )
              return
           end if
           ! initialize q and z if desired.
           if( icompq==3_${ik}$ )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, q, ldq )
           if( icompz==3_${ik}$ )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, z, ldz )
           ! quick return if possible
           if( n<=1 )return
           ! zero out lower triangle of b
           do jcol = 1, n - 1
              do jrow = jcol + 1, n
                 b( jrow, jcol ) = zero
              end do
           end do
           ! reduce a and b
           do jcol = ilo, ihi - 2
              do jrow = ihi, jcol + 2, -1
                 ! step 1: rotate rows jrow-1, jrow to kill a(jrow,jcol)
                 temp = a( jrow-1, jcol )
                 call stdlib${ii}$_slartg( temp, a( jrow, jcol ), c, s,a( jrow-1, jcol ) )
                 a( jrow, jcol ) = zero
                 call stdlib${ii}$_srot( n-jcol, a( jrow-1, jcol+1 ), lda,a( jrow, jcol+1 ), lda, c, s )
                           
                 call stdlib${ii}$_srot( n+2-jrow, b( jrow-1, jrow-1 ), ldb,b( jrow, jrow-1 ), ldb, c, &
                           s )
                 if( ilq )call stdlib${ii}$_srot( n, q( 1_${ik}$, jrow-1 ), 1_${ik}$, q( 1_${ik}$, jrow ), 1_${ik}$, c, s )
                 ! step 2: rotate columns jrow, jrow-1 to kill b(jrow,jrow-1)
                 temp = b( jrow, jrow )
                 call stdlib${ii}$_slartg( temp, b( jrow, jrow-1 ), c, s,b( jrow, jrow ) )
                 b( jrow, jrow-1 ) = zero
                 call stdlib${ii}$_srot( ihi, a( 1_${ik}$, jrow ), 1_${ik}$, a( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
                 call stdlib${ii}$_srot( jrow-1, b( 1_${ik}$, jrow ), 1_${ik}$, b( 1_${ik}$, jrow-1 ), 1_${ik}$, c,s )
                 if( ilz )call stdlib${ii}$_srot( n, z( 1_${ik}$, jrow ), 1_${ik}$, z( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
              end do
           end do
           return
     end subroutine stdlib${ii}$_sgghrd

     pure module subroutine stdlib${ii}$_dgghrd( compq, compz, n, ilo, ihi, a, lda, b, ldb, q,ldq, z, ldz, &
     !! DGGHRD reduces a pair of real matrices (A,B) to generalized upper
     !! Hessenberg form using orthogonal transformations, where A is a
     !! general matrix and B is upper triangular.  The form of the
     !! generalized eigenvalue problem is
     !! A*x = lambda*B*x,
     !! and B is typically made upper triangular by computing its QR
     !! factorization and moving the orthogonal matrix Q to the left side
     !! of the equation.
     !! This subroutine simultaneously reduces A to a Hessenberg matrix H:
     !! Q**T*A*Z = H
     !! and transforms B to another upper triangular matrix T:
     !! Q**T*B*Z = T
     !! in order to reduce the problem to its standard form
     !! H*y = lambda*T*y
     !! where y = Z**T*x.
     !! The orthogonal matrices Q and Z are determined as products of Givens
     !! rotations.  They may either be formed explicitly, or they may be
     !! postmultiplied into input matrices Q1 and Z1, so that
     !! Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T
     !! Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T
     !! If Q1 is the orthogonal matrix from the QR factorization of B in the
     !! original equation A*x = lambda*B*x, then DGGHRD reduces the original
     !! problem to generalized Hessenberg form.
               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) :: compq, compz
           integer(${ik}$), intent(in) :: ihi, ilo, lda, ldb, ldq, ldz, n
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           real(dp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: ilq, ilz
           integer(${ik}$) :: icompq, icompz, jcol, jrow
           real(dp) :: c, s, temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode compq
           if( stdlib_lsame( compq, 'N' ) ) then
              ilq = .false.
              icompq = 1_${ik}$
           else if( stdlib_lsame( compq, 'V' ) ) then
              ilq = .true.
              icompq = 2_${ik}$
           else if( stdlib_lsame( compq, 'I' ) ) then
              ilq = .true.
              icompq = 3_${ik}$
           else
              icompq = 0_${ik}$
           end if
           ! decode compz
           if( stdlib_lsame( compz, 'N' ) ) then
              ilz = .false.
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              ilz = .true.
              icompz = 2_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              ilz = .true.
              icompz = 3_${ik}$
           else
              icompz = 0_${ik}$
           end if
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( icompz<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ ) then
              info = -4_${ik}$
           else if( ihi>n .or. ihi<ilo-1 ) 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( ( ilq .and. ldq<n ) .or. ldq<1_${ik}$ ) then
              info = -11_${ik}$
           else if( ( ilz .and. ldz<n ) .or. ldz<1_${ik}$ ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGGHRD', -info )
              return
           end if
           ! initialize q and z if desired.
           if( icompq==3_${ik}$ )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, q, ldq )
           if( icompz==3_${ik}$ )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, z, ldz )
           ! quick return if possible
           if( n<=1 )return
           ! zero out lower triangle of b
           do jcol = 1, n - 1
              do jrow = jcol + 1, n
                 b( jrow, jcol ) = zero
              end do
           end do
           ! reduce a and b
           do jcol = ilo, ihi - 2
              do jrow = ihi, jcol + 2, -1
                 ! step 1: rotate rows jrow-1, jrow to kill a(jrow,jcol)
                 temp = a( jrow-1, jcol )
                 call stdlib${ii}$_dlartg( temp, a( jrow, jcol ), c, s,a( jrow-1, jcol ) )
                 a( jrow, jcol ) = zero
                 call stdlib${ii}$_drot( n-jcol, a( jrow-1, jcol+1 ), lda,a( jrow, jcol+1 ), lda, c, s )
                           
                 call stdlib${ii}$_drot( n+2-jrow, b( jrow-1, jrow-1 ), ldb,b( jrow, jrow-1 ), ldb, c, &
                           s )
                 if( ilq )call stdlib${ii}$_drot( n, q( 1_${ik}$, jrow-1 ), 1_${ik}$, q( 1_${ik}$, jrow ), 1_${ik}$, c, s )
                 ! step 2: rotate columns jrow, jrow-1 to kill b(jrow,jrow-1)
                 temp = b( jrow, jrow )
                 call stdlib${ii}$_dlartg( temp, b( jrow, jrow-1 ), c, s,b( jrow, jrow ) )
                 b( jrow, jrow-1 ) = zero
                 call stdlib${ii}$_drot( ihi, a( 1_${ik}$, jrow ), 1_${ik}$, a( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
                 call stdlib${ii}$_drot( jrow-1, b( 1_${ik}$, jrow ), 1_${ik}$, b( 1_${ik}$, jrow-1 ), 1_${ik}$, c,s )
                 if( ilz )call stdlib${ii}$_drot( n, z( 1_${ik}$, jrow ), 1_${ik}$, z( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
              end do
           end do
           return
     end subroutine stdlib${ii}$_dgghrd

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$gghrd( compq, compz, n, ilo, ihi, a, lda, b, ldb, q,ldq, z, ldz, &
     !! DGGHRD: reduces a pair of real matrices (A,B) to generalized upper
     !! Hessenberg form using orthogonal transformations, where A is a
     !! general matrix and B is upper triangular.  The form of the
     !! generalized eigenvalue problem is
     !! A*x = lambda*B*x,
     !! and B is typically made upper triangular by computing its QR
     !! factorization and moving the orthogonal matrix Q to the left side
     !! of the equation.
     !! This subroutine simultaneously reduces A to a Hessenberg matrix H:
     !! Q**T*A*Z = H
     !! and transforms B to another upper triangular matrix T:
     !! Q**T*B*Z = T
     !! in order to reduce the problem to its standard form
     !! H*y = lambda*T*y
     !! where y = Z**T*x.
     !! The orthogonal matrices Q and Z are determined as products of Givens
     !! rotations.  They may either be formed explicitly, or they may be
     !! postmultiplied into input matrices Q1 and Z1, so that
     !! Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T
     !! Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T
     !! If Q1 is the orthogonal matrix from the QR factorization of B in the
     !! original equation A*x = lambda*B*x, then DGGHRD reduces the original
     !! problem to generalized Hessenberg form.
               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) :: compq, compz
           integer(${ik}$), intent(in) :: ihi, ilo, lda, ldb, ldq, ldz, n
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: ilq, ilz
           integer(${ik}$) :: icompq, icompz, jcol, jrow
           real(${rk}$) :: c, s, temp
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode compq
           if( stdlib_lsame( compq, 'N' ) ) then
              ilq = .false.
              icompq = 1_${ik}$
           else if( stdlib_lsame( compq, 'V' ) ) then
              ilq = .true.
              icompq = 2_${ik}$
           else if( stdlib_lsame( compq, 'I' ) ) then
              ilq = .true.
              icompq = 3_${ik}$
           else
              icompq = 0_${ik}$
           end if
           ! decode compz
           if( stdlib_lsame( compz, 'N' ) ) then
              ilz = .false.
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              ilz = .true.
              icompz = 2_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              ilz = .true.
              icompz = 3_${ik}$
           else
              icompz = 0_${ik}$
           end if
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( icompz<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ ) then
              info = -4_${ik}$
           else if( ihi>n .or. ihi<ilo-1 ) 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( ( ilq .and. ldq<n ) .or. ldq<1_${ik}$ ) then
              info = -11_${ik}$
           else if( ( ilz .and. ldz<n ) .or. ldz<1_${ik}$ ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DGGHRD', -info )
              return
           end if
           ! initialize q and z if desired.
           if( icompq==3_${ik}$ )call stdlib${ii}$_${ri}$laset( 'FULL', n, n, zero, one, q, ldq )
           if( icompz==3_${ik}$ )call stdlib${ii}$_${ri}$laset( 'FULL', n, n, zero, one, z, ldz )
           ! quick return if possible
           if( n<=1 )return
           ! zero out lower triangle of b
           do jcol = 1, n - 1
              do jrow = jcol + 1, n
                 b( jrow, jcol ) = zero
              end do
           end do
           ! reduce a and b
           do jcol = ilo, ihi - 2
              do jrow = ihi, jcol + 2, -1
                 ! step 1: rotate rows jrow-1, jrow to kill a(jrow,jcol)
                 temp = a( jrow-1, jcol )
                 call stdlib${ii}$_${ri}$lartg( temp, a( jrow, jcol ), c, s,a( jrow-1, jcol ) )
                 a( jrow, jcol ) = zero
                 call stdlib${ii}$_${ri}$rot( n-jcol, a( jrow-1, jcol+1 ), lda,a( jrow, jcol+1 ), lda, c, s )
                           
                 call stdlib${ii}$_${ri}$rot( n+2-jrow, b( jrow-1, jrow-1 ), ldb,b( jrow, jrow-1 ), ldb, c, &
                           s )
                 if( ilq )call stdlib${ii}$_${ri}$rot( n, q( 1_${ik}$, jrow-1 ), 1_${ik}$, q( 1_${ik}$, jrow ), 1_${ik}$, c, s )
                 ! step 2: rotate columns jrow, jrow-1 to kill b(jrow,jrow-1)
                 temp = b( jrow, jrow )
                 call stdlib${ii}$_${ri}$lartg( temp, b( jrow, jrow-1 ), c, s,b( jrow, jrow ) )
                 b( jrow, jrow-1 ) = zero
                 call stdlib${ii}$_${ri}$rot( ihi, a( 1_${ik}$, jrow ), 1_${ik}$, a( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
                 call stdlib${ii}$_${ri}$rot( jrow-1, b( 1_${ik}$, jrow ), 1_${ik}$, b( 1_${ik}$, jrow-1 ), 1_${ik}$, c,s )
                 if( ilz )call stdlib${ii}$_${ri}$rot( n, z( 1_${ik}$, jrow ), 1_${ik}$, z( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
              end do
           end do
           return
     end subroutine stdlib${ii}$_${ri}$gghrd

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_cgghrd( compq, compz, n, ilo, ihi, a, lda, b, ldb, q,ldq, z, ldz, &
     !! CGGHRD reduces a pair of complex matrices (A,B) to generalized upper
     !! Hessenberg form using unitary transformations, where A is a
     !! general matrix and B is upper triangular.  The form of the generalized
     !! eigenvalue problem is
     !! A*x = lambda*B*x,
     !! and B is typically made upper triangular by computing its QR
     !! factorization and moving the unitary matrix Q to the left side
     !! of the equation.
     !! This subroutine simultaneously reduces A to a Hessenberg matrix H:
     !! Q**H*A*Z = H
     !! and transforms B to another upper triangular matrix T:
     !! Q**H*B*Z = T
     !! in order to reduce the problem to its standard form
     !! H*y = lambda*T*y
     !! where y = Z**H*x.
     !! The unitary matrices Q and Z are determined as products of Givens
     !! rotations.  They may either be formed explicitly, or they may be
     !! postmultiplied into input matrices Q1 and Z1, so that
     !! Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
     !! Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
     !! If Q1 is the unitary matrix from the QR factorization of B in the
     !! original equation A*x = lambda*B*x, then CGGHRD reduces the original
     !! problem to generalized Hessenberg form.
               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) :: compq, compz
           integer(${ik}$), intent(in) :: ihi, ilo, lda, ldb, ldq, ldz, n
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           complex(sp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: ilq, ilz
           integer(${ik}$) :: icompq, icompz, jcol, jrow
           real(sp) :: c
           complex(sp) :: ctemp, s
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode compq
           if( stdlib_lsame( compq, 'N' ) ) then
              ilq = .false.
              icompq = 1_${ik}$
           else if( stdlib_lsame( compq, 'V' ) ) then
              ilq = .true.
              icompq = 2_${ik}$
           else if( stdlib_lsame( compq, 'I' ) ) then
              ilq = .true.
              icompq = 3_${ik}$
           else
              icompq = 0_${ik}$
           end if
           ! decode compz
           if( stdlib_lsame( compz, 'N' ) ) then
              ilz = .false.
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              ilz = .true.
              icompz = 2_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              ilz = .true.
              icompz = 3_${ik}$
           else
              icompz = 0_${ik}$
           end if
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( icompz<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ ) then
              info = -4_${ik}$
           else if( ihi>n .or. ihi<ilo-1 ) 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( ( ilq .and. ldq<n ) .or. ldq<1_${ik}$ ) then
              info = -11_${ik}$
           else if( ( ilz .and. ldz<n ) .or. ldz<1_${ik}$ ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CGGHRD', -info )
              return
           end if
           ! initialize q and z if desired.
           if( icompq==3_${ik}$ )call stdlib${ii}$_claset( 'FULL', n, n, czero, cone, q, ldq )
           if( icompz==3_${ik}$ )call stdlib${ii}$_claset( 'FULL', n, n, czero, cone, z, ldz )
           ! quick return if possible
           if( n<=1 )return
           ! zero out lower triangle of b
           do jcol = 1, n - 1
              do jrow = jcol + 1, n
                 b( jrow, jcol ) = czero
              end do
           end do
           ! reduce a and b
           do jcol = ilo, ihi - 2
              do jrow = ihi, jcol + 2, -1
                 ! step 1: rotate rows jrow-1, jrow to kill a(jrow,jcol)
                 ctemp = a( jrow-1, jcol )
                 call stdlib${ii}$_clartg( ctemp, a( jrow, jcol ), c, s,a( jrow-1, jcol ) )
                 a( jrow, jcol ) = czero
                 call stdlib${ii}$_crot( n-jcol, a( jrow-1, jcol+1 ), lda,a( jrow, jcol+1 ), lda, c, s )
                           
                 call stdlib${ii}$_crot( n+2-jrow, b( jrow-1, jrow-1 ), ldb,b( jrow, jrow-1 ), ldb, c, &
                           s )
                 if( ilq )call stdlib${ii}$_crot( n, q( 1_${ik}$, jrow-1 ), 1_${ik}$, q( 1_${ik}$, jrow ), 1_${ik}$, c,conjg( s ) )
                           
                 ! step 2: rotate columns jrow, jrow-1 to kill b(jrow,jrow-1)
                 ctemp = b( jrow, jrow )
                 call stdlib${ii}$_clartg( ctemp, b( jrow, jrow-1 ), c, s,b( jrow, jrow ) )
                 b( jrow, jrow-1 ) = czero
                 call stdlib${ii}$_crot( ihi, a( 1_${ik}$, jrow ), 1_${ik}$, a( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
                 call stdlib${ii}$_crot( jrow-1, b( 1_${ik}$, jrow ), 1_${ik}$, b( 1_${ik}$, jrow-1 ), 1_${ik}$, c,s )
                 if( ilz )call stdlib${ii}$_crot( n, z( 1_${ik}$, jrow ), 1_${ik}$, z( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
              end do
           end do
           return
     end subroutine stdlib${ii}$_cgghrd

     pure module subroutine stdlib${ii}$_zgghrd( compq, compz, n, ilo, ihi, a, lda, b, ldb, q,ldq, z, ldz, &
     !! ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper
     !! Hessenberg form using unitary transformations, where A is a
     !! general matrix and B is upper triangular.  The form of the
     !! generalized eigenvalue problem is
     !! A*x = lambda*B*x,
     !! and B is typically made upper triangular by computing its QR
     !! factorization and moving the unitary matrix Q to the left side
     !! of the equation.
     !! This subroutine simultaneously reduces A to a Hessenberg matrix H:
     !! Q**H*A*Z = H
     !! and transforms B to another upper triangular matrix T:
     !! Q**H*B*Z = T
     !! in order to reduce the problem to its standard form
     !! H*y = lambda*T*y
     !! where y = Z**H*x.
     !! The unitary matrices Q and Z are determined as products of Givens
     !! rotations.  They may either be formed explicitly, or they may be
     !! postmultiplied into input matrices Q1 and Z1, so that
     !! Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
     !! Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
     !! If Q1 is the unitary matrix from the QR factorization of B in the
     !! original equation A*x = lambda*B*x, then ZGGHRD reduces the original
     !! problem to generalized Hessenberg form.
               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) :: compq, compz
           integer(${ik}$), intent(in) :: ihi, ilo, lda, ldb, ldq, ldz, n
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           complex(dp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: ilq, ilz
           integer(${ik}$) :: icompq, icompz, jcol, jrow
           real(dp) :: c
           complex(dp) :: ctemp, s
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode compq
           if( stdlib_lsame( compq, 'N' ) ) then
              ilq = .false.
              icompq = 1_${ik}$
           else if( stdlib_lsame( compq, 'V' ) ) then
              ilq = .true.
              icompq = 2_${ik}$
           else if( stdlib_lsame( compq, 'I' ) ) then
              ilq = .true.
              icompq = 3_${ik}$
           else
              icompq = 0_${ik}$
           end if
           ! decode compz
           if( stdlib_lsame( compz, 'N' ) ) then
              ilz = .false.
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              ilz = .true.
              icompz = 2_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              ilz = .true.
              icompz = 3_${ik}$
           else
              icompz = 0_${ik}$
           end if
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( icompz<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ ) then
              info = -4_${ik}$
           else if( ihi>n .or. ihi<ilo-1 ) 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( ( ilq .and. ldq<n ) .or. ldq<1_${ik}$ ) then
              info = -11_${ik}$
           else if( ( ilz .and. ldz<n ) .or. ldz<1_${ik}$ ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGGHRD', -info )
              return
           end if
           ! initialize q and z if desired.
           if( icompq==3_${ik}$ )call stdlib${ii}$_zlaset( 'FULL', n, n, czero, cone, q, ldq )
           if( icompz==3_${ik}$ )call stdlib${ii}$_zlaset( 'FULL', n, n, czero, cone, z, ldz )
           ! quick return if possible
           if( n<=1 )return
           ! zero out lower triangle of b
           do jcol = 1, n - 1
              do jrow = jcol + 1, n
                 b( jrow, jcol ) = czero
              end do
           end do
           ! reduce a and b
           do jcol = ilo, ihi - 2
              do jrow = ihi, jcol + 2, -1
                 ! step 1: rotate rows jrow-1, jrow to kill a(jrow,jcol)
                 ctemp = a( jrow-1, jcol )
                 call stdlib${ii}$_zlartg( ctemp, a( jrow, jcol ), c, s,a( jrow-1, jcol ) )
                 a( jrow, jcol ) = czero
                 call stdlib${ii}$_zrot( n-jcol, a( jrow-1, jcol+1 ), lda,a( jrow, jcol+1 ), lda, c, s )
                           
                 call stdlib${ii}$_zrot( n+2-jrow, b( jrow-1, jrow-1 ), ldb,b( jrow, jrow-1 ), ldb, c, &
                           s )
                 if( ilq )call stdlib${ii}$_zrot( n, q( 1_${ik}$, jrow-1 ), 1_${ik}$, q( 1_${ik}$, jrow ), 1_${ik}$, c,conjg( s ) )
                           
                 ! step 2: rotate columns jrow, jrow-1 to kill b(jrow,jrow-1)
                 ctemp = b( jrow, jrow )
                 call stdlib${ii}$_zlartg( ctemp, b( jrow, jrow-1 ), c, s,b( jrow, jrow ) )
                 b( jrow, jrow-1 ) = czero
                 call stdlib${ii}$_zrot( ihi, a( 1_${ik}$, jrow ), 1_${ik}$, a( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
                 call stdlib${ii}$_zrot( jrow-1, b( 1_${ik}$, jrow ), 1_${ik}$, b( 1_${ik}$, jrow-1 ), 1_${ik}$, c,s )
                 if( ilz )call stdlib${ii}$_zrot( n, z( 1_${ik}$, jrow ), 1_${ik}$, z( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
              end do
           end do
           return
     end subroutine stdlib${ii}$_zgghrd

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$gghrd( compq, compz, n, ilo, ihi, a, lda, b, ldb, q,ldq, z, ldz, &
     !! ZGGHRD: reduces a pair of complex matrices (A,B) to generalized upper
     !! Hessenberg form using unitary transformations, where A is a
     !! general matrix and B is upper triangular.  The form of the
     !! generalized eigenvalue problem is
     !! A*x = lambda*B*x,
     !! and B is typically made upper triangular by computing its QR
     !! factorization and moving the unitary matrix Q to the left side
     !! of the equation.
     !! This subroutine simultaneously reduces A to a Hessenberg matrix H:
     !! Q**H*A*Z = H
     !! and transforms B to another upper triangular matrix T:
     !! Q**H*B*Z = T
     !! in order to reduce the problem to its standard form
     !! H*y = lambda*T*y
     !! where y = Z**H*x.
     !! The unitary matrices Q and Z are determined as products of Givens
     !! rotations.  They may either be formed explicitly, or they may be
     !! postmultiplied into input matrices Q1 and Z1, so that
     !! Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
     !! Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
     !! If Q1 is the unitary matrix from the QR factorization of B in the
     !! original equation A*x = lambda*B*x, then ZGGHRD reduces the original
     !! problem to generalized Hessenberg form.
               info )
        ! -- lapack computational routine --
        ! -- lapack is a software package provided by univ. of tennessee,    --
        ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
           use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
           ! Scalar Arguments 
           character, intent(in) :: compq, compz
           integer(${ik}$), intent(in) :: ihi, ilo, lda, ldb, ldq, ldz, n
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: ilq, ilz
           integer(${ik}$) :: icompq, icompz, jcol, jrow
           real(${ck}$) :: c
           complex(${ck}$) :: ctemp, s
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode compq
           if( stdlib_lsame( compq, 'N' ) ) then
              ilq = .false.
              icompq = 1_${ik}$
           else if( stdlib_lsame( compq, 'V' ) ) then
              ilq = .true.
              icompq = 2_${ik}$
           else if( stdlib_lsame( compq, 'I' ) ) then
              ilq = .true.
              icompq = 3_${ik}$
           else
              icompq = 0_${ik}$
           end if
           ! decode compz
           if( stdlib_lsame( compz, 'N' ) ) then
              ilz = .false.
              icompz = 1_${ik}$
           else if( stdlib_lsame( compz, 'V' ) ) then
              ilz = .true.
              icompz = 2_${ik}$
           else if( stdlib_lsame( compz, 'I' ) ) then
              ilz = .true.
              icompz = 3_${ik}$
           else
              icompz = 0_${ik}$
           end if
           ! test the input parameters.
           info = 0_${ik}$
           if( icompq<=0_${ik}$ ) then
              info = -1_${ik}$
           else if( icompz<=0_${ik}$ ) then
              info = -2_${ik}$
           else if( n<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ilo<1_${ik}$ ) then
              info = -4_${ik}$
           else if( ihi>n .or. ihi<ilo-1 ) 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( ( ilq .and. ldq<n ) .or. ldq<1_${ik}$ ) then
              info = -11_${ik}$
           else if( ( ilz .and. ldz<n ) .or. ldz<1_${ik}$ ) then
              info = -13_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZGGHRD', -info )
              return
           end if
           ! initialize q and z if desired.
           if( icompq==3_${ik}$ )call stdlib${ii}$_${ci}$laset( 'FULL', n, n, czero, cone, q, ldq )
           if( icompz==3_${ik}$ )call stdlib${ii}$_${ci}$laset( 'FULL', n, n, czero, cone, z, ldz )
           ! quick return if possible
           if( n<=1 )return
           ! zero out lower triangle of b
           do jcol = 1, n - 1
              do jrow = jcol + 1, n
                 b( jrow, jcol ) = czero
              end do
           end do
           ! reduce a and b
           do jcol = ilo, ihi - 2
              do jrow = ihi, jcol + 2, -1
                 ! step 1: rotate rows jrow-1, jrow to kill a(jrow,jcol)
                 ctemp = a( jrow-1, jcol )
                 call stdlib${ii}$_${ci}$lartg( ctemp, a( jrow, jcol ), c, s,a( jrow-1, jcol ) )
                 a( jrow, jcol ) = czero
                 call stdlib${ii}$_${ci}$rot( n-jcol, a( jrow-1, jcol+1 ), lda,a( jrow, jcol+1 ), lda, c, s )
                           
                 call stdlib${ii}$_${ci}$rot( n+2-jrow, b( jrow-1, jrow-1 ), ldb,b( jrow, jrow-1 ), ldb, c, &
                           s )
                 if( ilq )call stdlib${ii}$_${ci}$rot( n, q( 1_${ik}$, jrow-1 ), 1_${ik}$, q( 1_${ik}$, jrow ), 1_${ik}$, c,conjg( s ) )
                           
                 ! step 2: rotate columns jrow, jrow-1 to kill b(jrow,jrow-1)
                 ctemp = b( jrow, jrow )
                 call stdlib${ii}$_${ci}$lartg( ctemp, b( jrow, jrow-1 ), c, s,b( jrow, jrow ) )
                 b( jrow, jrow-1 ) = czero
                 call stdlib${ii}$_${ci}$rot( ihi, a( 1_${ik}$, jrow ), 1_${ik}$, a( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
                 call stdlib${ii}$_${ci}$rot( jrow-1, b( 1_${ik}$, jrow ), 1_${ik}$, b( 1_${ik}$, jrow-1 ), 1_${ik}$, c,s )
                 if( ilz )call stdlib${ii}$_${ci}$rot( n, z( 1_${ik}$, jrow ), 1_${ik}$, z( 1_${ik}$, jrow-1 ), 1_${ik}$, c, s )
              end do
           end do
           return
     end subroutine stdlib${ii}$_${ci}$gghrd

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_sgghd3( compq, compz, n, ilo, ihi, a, lda, b, ldb, q,ldq, z, ldz, &
     !! SGGHD3 reduces a pair of real matrices (A,B) to generalized upper
     !! Hessenberg form using orthogonal transformations, where A is a
     !! general matrix and B is upper triangular.  The form of the
     !! generalized eigenvalue problem is
     !! A*x = lambda*B*x,
     !! and B is typically made upper triangular by computing its QR
     !! factorization and moving the orthogonal matrix Q to the left side
     !! of the equation.
     !! This subroutine simultaneously reduces A to a Hessenberg matrix H:
     !! Q**T*A*Z = H
     !! and transforms B to another upper triangular matrix T:
     !! Q**T*B*Z = T
     !! in order to reduce the problem to its standard form
     !! H*y = lambda*T*y
     !! where y = Z**T*x.
     !! The orthogonal matrices Q and Z are determined as products of Givens
     !! rotations.  They may either be formed explicitly, or they may be
     !! postmultiplied into input matrices Q1 and Z1, so that
     !! Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T
     !! Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T
     !! If Q1 is the orthogonal matrix from the QR factorization of B in the
     !! original equation A*x = lambda*B*x, then SGGHD3 reduces the original
     !! problem to generalized Hessenberg form.
     !! This is a blocked variant of SGGHRD, using matrix-matrix
     !! multiplications for parts of the computation to enhance performance.
               work, lwork, 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) :: compq, compz
           integer(${ik}$), intent(in) :: ihi, ilo, lda, ldb, ldq, ldz, n, lwork
           integer(${ik}$), intent(out) :: info
           ! Array Arguments 
           real(sp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: blk22, initq, initz, lquery, wantq, wantz
           character :: compq2, compz2
           integer(${ik}$) :: cola, i, ierr, j, j0, jcol, jj, jrow, k, kacc22, len, lwkopt, n2nb, nb,&
                      nblst, nbmin, nh, nnb, nx, ppw, ppwo, pw, top, topq
           real(sp) :: c, c1, c2, s, s1, s2, temp, temp1, temp2, temp3
           ! Intrinsic Functions 
           ! Executable Statements 
           ! decode and test the input parameters.
           info = 0_${ik}$
           nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SGGHD3',