#: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',