#:include "common.fypp" submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_gen3 implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES module subroutine stdlib${ii}$_slaqtr( ltran, lreal, n, t, ldt, b, w, scale, x, work,info ) !! SLAQTR solves the real quasi-triangular system !! op(T)*p = scale*c, if LREAL = .TRUE. !! or the complex quasi-triangular systems !! op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE. !! in real arithmetic, where T is upper quasi-triangular. !! If LREAL = .FALSE., then the first diagonal block of T must be !! 1 by 1, B is the specially structured matrix !! B = [ b(1) b(2) ... b(n) ] !! [ w ] !! [ w ] !! [ . ] !! [ w ] !! op(A) = A or A**T, A**T denotes the transpose of !! matrix A. !! On input, X = [ c ]. On output, X = [ p ]. !! [ d ] [ q ] !! This subroutine is designed for the condition number estimation !! in routine STRSNA. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: lreal, ltran integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldt, n real(sp), intent(out) :: scale real(sp), intent(in) :: w ! Array Arguments real(sp), intent(in) :: b(*), t(ldt,*) real(sp), intent(out) :: work(*) real(sp), intent(inout) :: x(*) ! ===================================================================== ! Local Scalars logical(lk) :: notran integer(${ik}$) :: i, ierr, j, j1, j2, jnext, k, n1, n2 real(sp) :: bignum, eps, rec, scaloc, si, smin, sminw, smlnum, sr, tjj, tmp, xj, xmax, & xnorm, z ! Local Arrays real(sp) :: d(2_${ik}$,2_${ik}$), v(2_${ik}$,2_${ik}$) ! Intrinsic Functions ! Executable Statements ! do not test the input parameters for errors notran = .not.ltran info = 0_${ik}$ ! quick return if possible if( n==0 )return ! set constants to control overflow eps = stdlib${ii}$_slamch( 'P' ) smlnum = stdlib${ii}$_slamch( 'S' ) / eps bignum = one / smlnum xnorm = stdlib${ii}$_slange( 'M', n, n, t, ldt, d ) if( .not.lreal )xnorm = max( xnorm, abs( w ), stdlib${ii}$_slange( 'M', n, 1_${ik}$, b, n, d ) ) smin = max( smlnum, eps*xnorm ) ! compute 1-norm of each column of strictly upper triangular ! part of t to control overflow in triangular solver. work( 1_${ik}$ ) = zero do j = 2, n work( j ) = stdlib${ii}$_sasum( j-1, t( 1_${ik}$, j ), 1_${ik}$ ) end do if( .not.lreal ) then do i = 2, n work( i ) = work( i ) + abs( b( i ) ) end do end if n2 = 2_${ik}$*n n1 = n if( .not.lreal )n1 = n2 k = stdlib${ii}$_isamax( n1, x, 1_${ik}$ ) xmax = abs( x( k ) ) scale = one if( xmax>bignum ) then scale = bignum / xmax call stdlib${ii}$_sscal( n1, scale, x, 1_${ik}$ ) xmax = bignum end if if( lreal ) then if( notran ) then ! solve t*p = scale*c jnext = n loop_30: do j = n, 1, -1 if( j>jnext )cycle loop_30 j1 = j j2 = j jnext = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnext = j - 2_${ik}$ end if end if if( j1==j2 ) then ! meet 1 by 1 diagonal block ! scale to avoid overflow when computing ! x(j) = b(j)/t(j,j) xj = abs( x( j1 ) ) tjj = abs( t( j1, j1 ) ) tmp = t( j1, j1 ) if( tjj<smin ) then tmp = smin tjj = smin info = 1_${ik}$ end if if( xj==zero )cycle loop_30 if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) / tmp xj = abs( x( j1 ) ) ! scale x if necessary to avoid overflow when adding a ! multiple of column j1 of t. if( xj>one ) then rec = one / xj if( work( j1 )>( bignum-xmax )*rec ) then call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ ) scale = scale*rec end if end if if( j1>1_${ik}$ ) then call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) k = stdlib${ii}$_isamax( j1-1, x, 1_${ik}$ ) xmax = abs( x( k ) ) end if else ! meet 2 by 2 diagonal block ! call 2 by 2 linear system solve, to take ! care of possible overflow by scaling factor. d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) call stdlib${ii}$_slaln2( .false., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d,& 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_sscal( n, scaloc, x, 1_${ik}$ ) scale = scale*scaloc end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) ! scale v(1,1) (= x(j1)) and/or v(2,1) (=x(j2)) ! to avoid overflow in updating right-hand side. xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) ), abs( v( 2_${ik}$, 1_${ik}$ ) ) ) if( xj>one ) then rec = one / xj if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ ) scale = scale*rec end if end if ! update right-hand side if( j1>1_${ik}$ ) then call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_saxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ ) k = stdlib${ii}$_isamax( j1-1, x, 1_${ik}$ ) xmax = abs( x( k ) ) end if end if end do loop_30 else ! solve t**t*p = scale*c jnext = 1_${ik}$ loop_40: do j = 1, n if( j<jnext )cycle loop_40 j1 = j j2 = j jnext = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnext = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = abs( x( j1 ) ) if( xmax>one ) then rec = one / xmax if( work( j1 )>( bignum-xj )*rec ) then call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) xj = abs( x( j1 ) ) tjj = abs( t( j1, j1 ) ) tmp = t( j1, j1 ) if( tjj<smin ) then tmp = smin tjj = smin info = 1_${ik}$ end if if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) / tmp xmax = max( xmax, abs( x( j1 ) ) ) else ! 2 by 2 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side elements by inner product. xj = max( abs( x( j1 ) ), abs( x( j2 ) ) ) if( xmax>one ) then rec = one / xmax if( max( work( j2 ), work( j1 ) )>( bignum-xj )*rec ) then call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ ) call stdlib${ii}$_slaln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d, & 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_sscal( n, scaloc, x, 1_${ik}$ ) scale = scale*scaloc end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax ) end if end do loop_40 end if else sminw = max( eps*abs( w ), smin ) if( notran ) then ! solve (t + ib)*(p+iq) = c+id jnext = n loop_70: do j = n, 1, -1 if( j>jnext )cycle loop_70 j1 = j j2 = j jnext = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnext = j - 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in division z = w if( j1==1_${ik}$ )z = b( 1_${ik}$ ) xj = abs( x( j1 ) ) + abs( x( n+j1 ) ) tjj = abs( t( j1, j1 ) ) + abs( z ) tmp = t( j1, j1 ) if( tjj<sminw ) then tmp = sminw tjj = sminw info = 1_${ik}$ end if if( xj==zero )cycle loop_70 if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if call stdlib${ii}$_sladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si ) x( j1 ) = sr x( n+j1 ) = si xj = abs( x( j1 ) ) + abs( x( n+j1 ) ) ! scale x if necessary to avoid overflow when adding a ! multiple of column j1 of t. if( xj>one ) then rec = one / xj if( work( j1 )>( bignum-xmax )*rec ) then call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec end if end if if( j1>1_${ik}$ ) then call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_saxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) xmax = zero do k = 1, j1 - 1 xmax = max( xmax, abs( x( k ) )+abs( x( k+n ) ) ) end do end if else ! meet 2 by 2 diagonal block d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) call stdlib${ii}$_slaln2( .false., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, & d, 2_${ik}$, zero, -w, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_sscal( 2_${ik}$*n, scaloc, x, 1_${ik}$ ) scale = scaloc*scale end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ ) x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ ) ! scale x(j1), .... to avoid overflow in ! updating right hand side. xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) )+abs( v( 1_${ik}$, 2_${ik}$ ) ),abs( v( 2_${ik}$, 1_${ik}$ ) )+abs( v( 2_${ik}$, 2_${ik}$ )& ) ) if( xj>one ) then rec = one / xj if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec end if end if ! update the right-hand side. if( j1>1_${ik}$ ) then call stdlib${ii}$_saxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_saxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_saxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) call stdlib${ii}$_saxpy( j1-1, -x( n+j2 ), t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) +b( j2 )*x( n+j2 ) x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -b( j2 )*x( j2 ) xmax = zero do k = 1, j1 - 1 xmax = max( abs( x( k ) )+abs( x( k+n ) ),xmax ) end do end if end if end do loop_70 else ! solve (t + ib)**t*(p+iq) = c+id jnext = 1_${ik}$ loop_80: do j = 1, n if( j<jnext )cycle loop_80 j1 = j j2 = j jnext = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnext = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = abs( x( j1 ) ) + abs( x( j1+n ) ) if( xmax>one ) then rec = one / xmax if( work( j1 )>( bignum-xj )*rec ) then call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) x( n+j1 ) = x( n+j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) if( j1>1_${ik}$ ) then x( j1 ) = x( j1 ) - b( j1 )*x( n+1 ) x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1_${ik}$ ) end if xj = abs( x( j1 ) ) + abs( x( j1+n ) ) z = w if( j1==1_${ik}$ )z = b( 1_${ik}$ ) ! scale if necessary to avoid overflow in ! complex division tjj = abs( t( j1, j1 ) ) + abs( z ) tmp = t( j1, j1 ) if( tjj<sminw ) then tmp = sminw tjj = sminw info = 1_${ik}$ end if if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if call stdlib${ii}$_sladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si ) x( j1 ) = sr x( j1+n ) = si xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax ) else ! 2 by 2 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) ) ) if( xmax>one ) then rec = one / xmax if( max( work( j1 ), work( j2 ) )>( bignum-xj ) / xmax ) then call stdlib${ii}$_sscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ ) d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) - stdlib${ii}$_sdot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) d( 1_${ik}$, 1_${ik}$ ) = d( 1_${ik}$, 1_${ik}$ ) - b( j1 )*x( n+1 ) d( 2_${ik}$, 1_${ik}$ ) = d( 2_${ik}$, 1_${ik}$ ) - b( j2 )*x( n+1 ) d( 1_${ik}$, 2_${ik}$ ) = d( 1_${ik}$, 2_${ik}$ ) + b( j1 )*x( 1_${ik}$ ) d( 2_${ik}$, 2_${ik}$ ) = d( 2_${ik}$, 2_${ik}$ ) + b( j2 )*x( 1_${ik}$ ) call stdlib${ii}$_slaln2( .true., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, d,& 2_${ik}$, zero, w, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_sscal( n2, scaloc, x, 1_${ik}$ ) scale = scaloc*scale end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ ) x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ ) xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) )& , xmax ) end if end do loop_80 end if end if return end subroutine stdlib${ii}$_slaqtr module subroutine stdlib${ii}$_dlaqtr( ltran, lreal, n, t, ldt, b, w, scale, x, work,info ) !! DLAQTR solves the real quasi-triangular system !! op(T)*p = scale*c, if LREAL = .TRUE. !! or the complex quasi-triangular systems !! op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE. !! in real arithmetic, where T is upper quasi-triangular. !! If LREAL = .FALSE., then the first diagonal block of T must be !! 1 by 1, B is the specially structured matrix !! B = [ b(1) b(2) ... b(n) ] !! [ w ] !! [ w ] !! [ . ] !! [ w ] !! op(A) = A or A**T, A**T denotes the transpose of !! matrix A. !! On input, X = [ c ]. On output, X = [ p ]. !! [ d ] [ q ] !! This subroutine is designed for the condition number estimation !! in routine DTRSNA. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: lreal, ltran integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldt, n real(dp), intent(out) :: scale real(dp), intent(in) :: w ! Array Arguments real(dp), intent(in) :: b(*), t(ldt,*) real(dp), intent(out) :: work(*) real(dp), intent(inout) :: x(*) ! ===================================================================== ! Local Scalars logical(lk) :: notran integer(${ik}$) :: i, ierr, j, j1, j2, jnext, k, n1, n2 real(dp) :: bignum, eps, rec, scaloc, si, smin, sminw, smlnum, sr, tjj, tmp, xj, xmax, & xnorm, z ! Local Arrays real(dp) :: d(2_${ik}$,2_${ik}$), v(2_${ik}$,2_${ik}$) ! Intrinsic Functions ! Executable Statements ! do not test the input parameters for errors notran = .not.ltran info = 0_${ik}$ ! quick return if possible if( n==0 )return ! set constants to control overflow eps = stdlib${ii}$_dlamch( 'P' ) smlnum = stdlib${ii}$_dlamch( 'S' ) / eps bignum = one / smlnum xnorm = stdlib${ii}$_dlange( 'M', n, n, t, ldt, d ) if( .not.lreal )xnorm = max( xnorm, abs( w ), stdlib${ii}$_dlange( 'M', n, 1_${ik}$, b, n, d ) ) smin = max( smlnum, eps*xnorm ) ! compute 1-norm of each column of strictly upper triangular ! part of t to control overflow in triangular solver. work( 1_${ik}$ ) = zero do j = 2, n work( j ) = stdlib${ii}$_dasum( j-1, t( 1_${ik}$, j ), 1_${ik}$ ) end do if( .not.lreal ) then do i = 2, n work( i ) = work( i ) + abs( b( i ) ) end do end if n2 = 2_${ik}$*n n1 = n if( .not.lreal )n1 = n2 k = stdlib${ii}$_idamax( n1, x, 1_${ik}$ ) xmax = abs( x( k ) ) scale = one if( xmax>bignum ) then scale = bignum / xmax call stdlib${ii}$_dscal( n1, scale, x, 1_${ik}$ ) xmax = bignum end if if( lreal ) then if( notran ) then ! solve t*p = scale*c jnext = n loop_30: do j = n, 1, -1 if( j>jnext )cycle loop_30 j1 = j j2 = j jnext = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnext = j - 2_${ik}$ end if end if if( j1==j2 ) then ! meet 1 by 1 diagonal block ! scale to avoid overflow when computing ! x(j) = b(j)/t(j,j) xj = abs( x( j1 ) ) tjj = abs( t( j1, j1 ) ) tmp = t( j1, j1 ) if( tjj<smin ) then tmp = smin tjj = smin info = 1_${ik}$ end if if( xj==zero )cycle loop_30 if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) / tmp xj = abs( x( j1 ) ) ! scale x if necessary to avoid overflow when adding a ! multiple of column j1 of t. if( xj>one ) then rec = one / xj if( work( j1 )>( bignum-xmax )*rec ) then call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ ) scale = scale*rec end if end if if( j1>1_${ik}$ ) then call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) k = stdlib${ii}$_idamax( j1-1, x, 1_${ik}$ ) xmax = abs( x( k ) ) end if else ! meet 2 by 2 diagonal block ! call 2 by 2 linear system solve, to take ! care of possible overflow by scaling factor. d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) call stdlib${ii}$_dlaln2( .false., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d,& 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_dscal( n, scaloc, x, 1_${ik}$ ) scale = scale*scaloc end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) ! scale v(1,1) (= x(j1)) and/or v(2,1) (=x(j2)) ! to avoid overflow in updating right-hand side. xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) ), abs( v( 2_${ik}$, 1_${ik}$ ) ) ) if( xj>one ) then rec = one / xj if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ ) scale = scale*rec end if end if ! update right-hand side if( j1>1_${ik}$ ) then call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_daxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ ) k = stdlib${ii}$_idamax( j1-1, x, 1_${ik}$ ) xmax = abs( x( k ) ) end if end if end do loop_30 else ! solve t**t*p = scale*c jnext = 1_${ik}$ loop_40: do j = 1, n if( j<jnext )cycle loop_40 j1 = j j2 = j jnext = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnext = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = abs( x( j1 ) ) if( xmax>one ) then rec = one / xmax if( work( j1 )>( bignum-xj )*rec ) then call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) xj = abs( x( j1 ) ) tjj = abs( t( j1, j1 ) ) tmp = t( j1, j1 ) if( tjj<smin ) then tmp = smin tjj = smin info = 1_${ik}$ end if if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) / tmp xmax = max( xmax, abs( x( j1 ) ) ) else ! 2 by 2 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side elements by inner product. xj = max( abs( x( j1 ) ), abs( x( j2 ) ) ) if( xmax>one ) then rec = one / xmax if( max( work( j2 ), work( j1 ) )>( bignum-xj )*rec ) then call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ ) call stdlib${ii}$_dlaln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d, & 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_dscal( n, scaloc, x, 1_${ik}$ ) scale = scale*scaloc end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax ) end if end do loop_40 end if else sminw = max( eps*abs( w ), smin ) if( notran ) then ! solve (t + ib)*(p+iq) = c+id jnext = n loop_70: do j = n, 1, -1 if( j>jnext )cycle loop_70 j1 = j j2 = j jnext = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnext = j - 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in division z = w if( j1==1_${ik}$ )z = b( 1_${ik}$ ) xj = abs( x( j1 ) ) + abs( x( n+j1 ) ) tjj = abs( t( j1, j1 ) ) + abs( z ) tmp = t( j1, j1 ) if( tjj<sminw ) then tmp = sminw tjj = sminw info = 1_${ik}$ end if if( xj==zero )cycle loop_70 if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if call stdlib${ii}$_dladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si ) x( j1 ) = sr x( n+j1 ) = si xj = abs( x( j1 ) ) + abs( x( n+j1 ) ) ! scale x if necessary to avoid overflow when adding a ! multiple of column j1 of t. if( xj>one ) then rec = one / xj if( work( j1 )>( bignum-xmax )*rec ) then call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec end if end if if( j1>1_${ik}$ ) then call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_daxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) xmax = zero do k = 1, j1 - 1 xmax = max( xmax, abs( x( k ) )+abs( x( k+n ) ) ) end do end if else ! meet 2 by 2 diagonal block d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) call stdlib${ii}$_dlaln2( .false., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, & d, 2_${ik}$, zero, -w, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_dscal( 2_${ik}$*n, scaloc, x, 1_${ik}$ ) scale = scaloc*scale end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ ) x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ ) ! scale x(j1), .... to avoid overflow in ! updating right hand side. xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) )+abs( v( 1_${ik}$, 2_${ik}$ ) ),abs( v( 2_${ik}$, 1_${ik}$ ) )+abs( v( 2_${ik}$, 2_${ik}$ )& ) ) if( xj>one ) then rec = one / xj if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec end if end if ! update the right-hand side. if( j1>1_${ik}$ ) then call stdlib${ii}$_daxpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_daxpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_daxpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) call stdlib${ii}$_daxpy( j1-1, -x( n+j2 ), t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) +b( j2 )*x( n+j2 ) x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -b( j2 )*x( j2 ) xmax = zero do k = 1, j1 - 1 xmax = max( abs( x( k ) )+abs( x( k+n ) ),xmax ) end do end if end if end do loop_70 else ! solve (t + ib)**t*(p+iq) = c+id jnext = 1_${ik}$ loop_80: do j = 1, n if( j<jnext )cycle loop_80 j1 = j j2 = j jnext = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnext = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = abs( x( j1 ) ) + abs( x( j1+n ) ) if( xmax>one ) then rec = one / xmax if( work( j1 )>( bignum-xj )*rec ) then call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) x( n+j1 ) = x( n+j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) if( j1>1_${ik}$ ) then x( j1 ) = x( j1 ) - b( j1 )*x( n+1 ) x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1_${ik}$ ) end if xj = abs( x( j1 ) ) + abs( x( j1+n ) ) z = w if( j1==1_${ik}$ )z = b( 1_${ik}$ ) ! scale if necessary to avoid overflow in ! complex division tjj = abs( t( j1, j1 ) ) + abs( z ) tmp = t( j1, j1 ) if( tjj<sminw ) then tmp = sminw tjj = sminw info = 1_${ik}$ end if if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if call stdlib${ii}$_dladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si ) x( j1 ) = sr x( j1+n ) = si xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax ) else ! 2 by 2 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) ) ) if( xmax>one ) then rec = one / xmax if( max( work( j1 ), work( j2 ) )>( bignum-xj ) / xmax ) then call stdlib${ii}$_dscal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ ) d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) - stdlib${ii}$_ddot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) d( 1_${ik}$, 1_${ik}$ ) = d( 1_${ik}$, 1_${ik}$ ) - b( j1 )*x( n+1 ) d( 2_${ik}$, 1_${ik}$ ) = d( 2_${ik}$, 1_${ik}$ ) - b( j2 )*x( n+1 ) d( 1_${ik}$, 2_${ik}$ ) = d( 1_${ik}$, 2_${ik}$ ) + b( j1 )*x( 1_${ik}$ ) d( 2_${ik}$, 2_${ik}$ ) = d( 2_${ik}$, 2_${ik}$ ) + b( j2 )*x( 1_${ik}$ ) call stdlib${ii}$_dlaln2( .true., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, d,& 2_${ik}$, zero, w, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_dscal( n2, scaloc, x, 1_${ik}$ ) scale = scaloc*scale end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ ) x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ ) xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) )& , xmax ) end if end do loop_80 end if end if return end subroutine stdlib${ii}$_dlaqtr #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] module subroutine stdlib${ii}$_${ri}$laqtr( ltran, lreal, n, t, ldt, b, w, scale, x, work,info ) !! DLAQTR: solves the real quasi-triangular system !! op(T)*p = scale*c, if LREAL = .TRUE. !! or the complex quasi-triangular systems !! op(T + iB)*(p+iq) = scale*(c+id), if LREAL = .FALSE. !! in real arithmetic, where T is upper quasi-triangular. !! If LREAL = .FALSE., then the first diagonal block of T must be !! 1 by 1, B is the specially structured matrix !! B = [ b(1) b(2) ... b(n) ] !! [ w ] !! [ w ] !! [ . ] !! [ w ] !! op(A) = A or A**T, A**T denotes the transpose of !! matrix A. !! On input, X = [ c ]. On output, X = [ p ]. !! [ d ] [ q ] !! This subroutine is designed for the condition number estimation !! in routine DTRSNA. ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: lreal, ltran integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldt, n real(${rk}$), intent(out) :: scale real(${rk}$), intent(in) :: w ! Array Arguments real(${rk}$), intent(in) :: b(*), t(ldt,*) real(${rk}$), intent(out) :: work(*) real(${rk}$), intent(inout) :: x(*) ! ===================================================================== ! Local Scalars logical(lk) :: notran integer(${ik}$) :: i, ierr, j, j1, j2, jnext, k, n1, n2 real(${rk}$) :: bignum, eps, rec, scaloc, si, smin, sminw, smlnum, sr, tjj, tmp, xj, xmax, & xnorm, z ! Local Arrays real(${rk}$) :: d(2_${ik}$,2_${ik}$), v(2_${ik}$,2_${ik}$) ! Intrinsic Functions ! Executable Statements ! do not test the input parameters for errors notran = .not.ltran info = 0_${ik}$ ! quick return if possible if( n==0 )return ! set constants to control overflow eps = stdlib${ii}$_${ri}$lamch( 'P' ) smlnum = stdlib${ii}$_${ri}$lamch( 'S' ) / eps bignum = one / smlnum xnorm = stdlib${ii}$_${ri}$lange( 'M', n, n, t, ldt, d ) if( .not.lreal )xnorm = max( xnorm, abs( w ), stdlib${ii}$_${ri}$lange( 'M', n, 1_${ik}$, b, n, d ) ) smin = max( smlnum, eps*xnorm ) ! compute 1-norm of each column of strictly upper triangular ! part of t to control overflow in triangular solver. work( 1_${ik}$ ) = zero do j = 2, n work( j ) = stdlib${ii}$_${ri}$asum( j-1, t( 1_${ik}$, j ), 1_${ik}$ ) end do if( .not.lreal ) then do i = 2, n work( i ) = work( i ) + abs( b( i ) ) end do end if n2 = 2_${ik}$*n n1 = n if( .not.lreal )n1 = n2 k = stdlib${ii}$_i${ri}$amax( n1, x, 1_${ik}$ ) xmax = abs( x( k ) ) scale = one if( xmax>bignum ) then scale = bignum / xmax call stdlib${ii}$_${ri}$scal( n1, scale, x, 1_${ik}$ ) xmax = bignum end if if( lreal ) then if( notran ) then ! solve t*p = scale*c jnext = n loop_30: do j = n, 1, -1 if( j>jnext )cycle loop_30 j1 = j j2 = j jnext = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnext = j - 2_${ik}$ end if end if if( j1==j2 ) then ! meet 1 by 1 diagonal block ! scale to avoid overflow when computing ! x(j) = b(j)/t(j,j) xj = abs( x( j1 ) ) tjj = abs( t( j1, j1 ) ) tmp = t( j1, j1 ) if( tjj<smin ) then tmp = smin tjj = smin info = 1_${ik}$ end if if( xj==zero )cycle loop_30 if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) / tmp xj = abs( x( j1 ) ) ! scale x if necessary to avoid overflow when adding a ! multiple of column j1 of t. if( xj>one ) then rec = one / xj if( work( j1 )>( bignum-xmax )*rec ) then call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ ) scale = scale*rec end if end if if( j1>1_${ik}$ ) then call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) k = stdlib${ii}$_i${ri}$amax( j1-1, x, 1_${ik}$ ) xmax = abs( x( k ) ) end if else ! meet 2 by 2 diagonal block ! call 2 by 2 linear system solve, to take ! care of possible overflow by scaling factor. d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) call stdlib${ii}$_${ri}$laln2( .false., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d,& 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_${ri}$scal( n, scaloc, x, 1_${ik}$ ) scale = scale*scaloc end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) ! scale v(1,1) (= x(j1)) and/or v(2,1) (=x(j2)) ! to avoid overflow in updating right-hand side. xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) ), abs( v( 2_${ik}$, 1_${ik}$ ) ) ) if( xj>one ) then rec = one / xj if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ ) scale = scale*rec end if end if ! update right-hand side if( j1>1_${ik}$ ) then call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_${ri}$axpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ ) k = stdlib${ii}$_i${ri}$amax( j1-1, x, 1_${ik}$ ) xmax = abs( x( k ) ) end if end if end do loop_30 else ! solve t**t*p = scale*c jnext = 1_${ik}$ loop_40: do j = 1, n if( j<jnext )cycle loop_40 j1 = j j2 = j jnext = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnext = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = abs( x( j1 ) ) if( xmax>one ) then rec = one / xmax if( work( j1 )>( bignum-xj )*rec ) then call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) xj = abs( x( j1 ) ) tjj = abs( t( j1, j1 ) ) tmp = t( j1, j1 ) if( tjj<smin ) then tmp = smin tjj = smin info = 1_${ik}$ end if if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) / tmp xmax = max( xmax, abs( x( j1 ) ) ) else ! 2 by 2 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side elements by inner product. xj = max( abs( x( j1 ) ), abs( x( j2 ) ) ) if( xmax>one ) then rec = one / xmax if( max( work( j2 ), work( j1 ) )>( bignum-xj )*rec ) then call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ ) call stdlib${ii}$_${ri}$laln2( .true., 2_${ik}$, 1_${ik}$, smin, one, t( j1, j1 ),ldt, one, one, d, & 2_${ik}$, zero, zero, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_${ri}$scal( n, scaloc, x, 1_${ik}$ ) scale = scale*scaloc end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax ) end if end do loop_40 end if else sminw = max( eps*abs( w ), smin ) if( notran ) then ! solve (t + ib)*(p+iq) = c+id jnext = n loop_70: do j = n, 1, -1 if( j>jnext )cycle loop_70 j1 = j j2 = j jnext = j - 1_${ik}$ if( j>1_${ik}$ ) then if( t( j, j-1 )/=zero ) then j1 = j - 1_${ik}$ jnext = j - 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in division z = w if( j1==1_${ik}$ )z = b( 1_${ik}$ ) xj = abs( x( j1 ) ) + abs( x( n+j1 ) ) tjj = abs( t( j1, j1 ) ) + abs( z ) tmp = t( j1, j1 ) if( tjj<sminw ) then tmp = sminw tjj = sminw info = 1_${ik}$ end if if( xj==zero )cycle loop_70 if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if call stdlib${ii}$_${ri}$ladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si ) x( j1 ) = sr x( n+j1 ) = si xj = abs( x( j1 ) ) + abs( x( n+j1 ) ) ! scale x if necessary to avoid overflow when adding a ! multiple of column j1 of t. if( xj>one ) then rec = one / xj if( work( j1 )>( bignum-xmax )*rec ) then call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ ) scale = scale*rec end if end if if( j1>1_${ik}$ ) then call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_${ri}$axpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) xmax = zero do k = 1, j1 - 1 xmax = max( xmax, abs( x( k ) )+abs( x( k+n ) ) ) end do end if else ! meet 2 by 2 diagonal block d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) call stdlib${ii}$_${ri}$laln2( .false., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, & d, 2_${ik}$, zero, -w, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_${ri}$scal( 2_${ik}$*n, scaloc, x, 1_${ik}$ ) scale = scaloc*scale end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ ) x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ ) ! scale x(j1), .... to avoid overflow in ! updating right hand side. xj = max( abs( v( 1_${ik}$, 1_${ik}$ ) )+abs( v( 1_${ik}$, 2_${ik}$ ) ),abs( v( 2_${ik}$, 1_${ik}$ ) )+abs( v( 2_${ik}$, 2_${ik}$ )& ) ) if( xj>one ) then rec = one / xj if( max( work( j1 ), work( j2 ) )>( bignum-xmax )*rec ) then call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ ) scale = scale*rec end if end if ! update the right-hand side. if( j1>1_${ik}$ ) then call stdlib${ii}$_${ri}$axpy( j1-1, -x( j1 ), t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_${ri}$axpy( j1-1, -x( j2 ), t( 1_${ik}$, j2 ), 1_${ik}$, x, 1_${ik}$ ) call stdlib${ii}$_${ri}$axpy( j1-1, -x( n+j1 ), t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$axpy( j1-1, -x( n+j2 ), t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) x( 1_${ik}$ ) = x( 1_${ik}$ ) + b( j1 )*x( n+j1 ) +b( j2 )*x( n+j2 ) x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -b( j2 )*x( j2 ) xmax = zero do k = 1, j1 - 1 xmax = max( abs( x( k ) )+abs( x( k+n ) ),xmax ) end do end if end if end do loop_70 else ! solve (t + ib)**t*(p+iq) = c+id jnext = 1_${ik}$ loop_80: do j = 1, n if( j<jnext )cycle loop_80 j1 = j j2 = j jnext = j + 1_${ik}$ if( j<n ) then if( t( j+1, j )/=zero ) then j2 = j + 1_${ik}$ jnext = j + 2_${ik}$ end if end if if( j1==j2 ) then ! 1 by 1 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = abs( x( j1 ) ) + abs( x( j1+n ) ) if( xmax>one ) then rec = one / xmax if( work( j1 )>( bignum-xj )*rec ) then call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if x( j1 ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x, 1_${ik}$ ) x( n+j1 ) = x( n+j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) if( j1>1_${ik}$ ) then x( j1 ) = x( j1 ) - b( j1 )*x( n+1 ) x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1_${ik}$ ) end if xj = abs( x( j1 ) ) + abs( x( j1+n ) ) z = w if( j1==1_${ik}$ )z = b( 1_${ik}$ ) ! scale if necessary to avoid overflow in ! complex division tjj = abs( t( j1, j1 ) ) + abs( z ) tmp = t( j1, j1 ) if( tjj<sminw ) then tmp = sminw tjj = sminw info = 1_${ik}$ end if if( tjj<one ) then if( xj>bignum*tjj ) then rec = one / xj call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if call stdlib${ii}$_${ri}$ladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si ) x( j1 ) = sr x( j1+n ) = si xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax ) else ! 2 by 2 diagonal block ! scale if necessary to avoid overflow in forming the ! right-hand side element by inner product. xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) ) ) if( xmax>one ) then rec = one / xmax if( max( work( j1 ), work( j2 ) )>( bignum-xj ) / xmax ) then call stdlib${ii}$_${ri}$scal( n2, rec, x, 1_${ik}$ ) scale = scale*rec xmax = xmax*rec end if end if d( 1_${ik}$, 1_${ik}$ ) = x( j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$, x,1_${ik}$ ) d( 2_${ik}$, 1_${ik}$ ) = x( j2 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$, x,1_${ik}$ ) d( 1_${ik}$, 2_${ik}$ ) = x( n+j1 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j1 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) d( 2_${ik}$, 2_${ik}$ ) = x( n+j2 ) - stdlib${ii}$_${ri}$dot( j1-1, t( 1_${ik}$, j2 ), 1_${ik}$,x( n+1 ), 1_${ik}$ ) d( 1_${ik}$, 1_${ik}$ ) = d( 1_${ik}$, 1_${ik}$ ) - b( j1 )*x( n+1 ) d( 2_${ik}$, 1_${ik}$ ) = d( 2_${ik}$, 1_${ik}$ ) - b( j2 )*x( n+1 ) d( 1_${ik}$, 2_${ik}$ ) = d( 1_${ik}$, 2_${ik}$ ) + b( j1 )*x( 1_${ik}$ ) d( 2_${ik}$, 2_${ik}$ ) = d( 2_${ik}$, 2_${ik}$ ) + b( j2 )*x( 1_${ik}$ ) call stdlib${ii}$_${ri}$laln2( .true., 2_${ik}$, 2_${ik}$, sminw, one, t( j1, j1 ),ldt, one, one, d,& 2_${ik}$, zero, w, v, 2_${ik}$,scaloc, xnorm, ierr ) if( ierr/=0_${ik}$ )info = 2_${ik}$ if( scaloc/=one ) then call stdlib${ii}$_${ri}$scal( n2, scaloc, x, 1_${ik}$ ) scale = scaloc*scale end if x( j1 ) = v( 1_${ik}$, 1_${ik}$ ) x( j2 ) = v( 2_${ik}$, 1_${ik}$ ) x( n+j1 ) = v( 1_${ik}$, 2_${ik}$ ) x( n+j2 ) = v( 2_${ik}$, 2_${ik}$ ) xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),abs( x( j2 ) )+abs( x( n+j2 ) )& , xmax ) end if end do loop_80 end if end if return end subroutine stdlib${ii}$_${ri}$laqtr #:endif #:endfor pure module subroutine stdlib${ii}$_slahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,iloz, ihiz, z, ldz, & !! SLAHQR is an auxiliary routine called by SHSEQR to update the !! eigenvalues and Schur decomposition already computed by SHSEQR, by !! dealing with the Hessenberg submatrix in rows and columns ILO to !! IHI. info ) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n integer(${ik}$), intent(out) :: info logical(lk), intent(in) :: wantt, wantz ! Array Arguments real(sp), intent(inout) :: h(ldh,*), z(ldz,*) real(sp), intent(out) :: wi(*), wr(*) ! ========================================================= ! Parameters real(sp), parameter :: dat1 = 3.0_sp/4.0_sp real(sp), parameter :: dat2 = -0.4375_sp integer(${ik}$), parameter :: kexsh = 10_${ik}$ ! Local Scalars real(sp) :: aa, ab, ba, bb, cs, det, h11, h12, h21, h21s, h22, rt1i, rt1r, rt2i, rt2r, & rtdisc, s, safmax, safmin, smlnum, sn, sum, t1, t2, t3, tr, tst, ulp, v2, v3 integer(${ik}$) :: i, i1, i2, its, itmax, j, k, l, m, nh, nr, nz, kdefl ! Local Arrays real(sp) :: v(3_${ik}$) ! Intrinsic Functions ! Executable Statements info = 0_${ik}$ ! quick return if possible if( n==0 )return if( ilo==ihi ) then wr( ilo ) = h( ilo, ilo ) wi( ilo ) = zero return end if ! ==== clear out the trash ==== do j = ilo, ihi - 3 h( j+2, j ) = zero h( j+3, j ) = zero end do if( ilo<=ihi-2 )h( ihi, ihi-2 ) = zero nh = ihi - ilo + 1_${ik}$ nz = ihiz - iloz + 1_${ik}$ ! set machine-dependent constants for the stopping criterion. safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) safmax = one / safmin call stdlib${ii}$_slabad( safmin, safmax ) ulp = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = safmin*( real( nh,KIND=sp) / ulp ) ! i1 and i2 are the indices of the first row and last column of h ! to which transformations must be applied. if eigenvalues only are ! being computed, i1 and i2 are set inside the main loop. if( wantt ) then i1 = 1_${ik}$ i2 = n end if ! itmax is the total number of qr iterations allowed. itmax = 30_${ik}$ * max( 10_${ik}$, nh ) ! kdefl counts the number of iterations since a deflation kdefl = 0_${ik}$ ! the main loop begins here. i is the loop index and decreases from ! ihi to ilo in steps of 1 or 2. each iteration of the loop works ! with the active submatrix in rows and columns l to i. ! eigenvalues i+1 to ihi have already converged. either l = ilo or ! h(l,l-1) is negligible so that the matrix splits. i = ihi 20 continue l = ilo if( i<ilo )go to 160 ! perform qr iterations on rows and columns ilo to i until a ! submatrix of order 1 or 2 splits off at the bottom because a ! subdiagonal element has become negligible. loop_140: do its = 0, itmax ! look for a single small subdiagonal element. do k = i, l + 1, -1 if( abs( h( k, k-1 ) )<=smlnum )go to 40 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) ) if( tst==zero ) then if( k-2>=ilo )tst = tst + abs( h( k-1, k-2 ) ) if( k+1<=ihi )tst = tst + abs( h( k+1, k ) ) end if ! ==== the following is a conservative small subdiagonal ! . deflation criterion due to ahues ! . 1997). it has better mathematical foundation and ! . improves accuracy in some cases. ==== if( abs( h( k, k-1 ) )<=ulp*tst ) then ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) ) ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) ) aa = max( abs( h( k, k ) ),abs( h( k-1, k-1 )-h( k, k ) ) ) bb = min( abs( h( k, k ) ),abs( h( k-1, k-1 )-h( k, k ) ) ) s = aa + ab if( ba*( ab / s )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 40 end if end do 40 continue l = k if( l>ilo ) then ! h(l,l-1) is negligible h( l, l-1 ) = zero end if ! exit from loop if a submatrix of order 1 or 2 has split off. if( l>=i-1 )go to 150 kdefl = kdefl + 1_${ik}$ ! now the active submatrix is in rows and columns l to i. if ! eigenvalues only are being computed, only the active submatrix ! need be transformed. if( .not.wantt ) then i1 = l i2 = i end if if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then ! exceptional shift. s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) ) h11 = dat1*s + h( i, i ) h12 = dat2*s h21 = s h22 = h11 else if( mod(kdefl,kexsh)==0_${ik}$ ) then ! exceptional shift. s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) ) h11 = dat1*s + h( l, l ) h12 = dat2*s h21 = s h22 = h11 else ! prepare to use francis' double shift ! (i.e. 2nd degree generalized rayleigh quotient) h11 = h( i-1, i-1 ) h21 = h( i, i-1 ) h12 = h( i-1, i ) h22 = h( i, i ) end if s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 ) if( s==zero ) then rt1r = zero rt1i = zero rt2r = zero rt2i = zero else h11 = h11 / s h21 = h21 / s h12 = h12 / s h22 = h22 / s tr = ( h11+h22 ) / two det = ( h11-tr )*( h22-tr ) - h12*h21 rtdisc = sqrt( abs( det ) ) if( det>=zero ) then ! ==== complex conjugate shifts ==== rt1r = tr*s rt2r = rt1r rt1i = rtdisc*s rt2i = -rt1i else ! ==== realshifts (use only one of them,KIND=sp) ==== rt1r = tr + rtdisc rt2r = tr - rtdisc if( abs( rt1r-h22 )<=abs( rt2r-h22 ) ) then rt1r = rt1r*s rt2r = rt1r else rt2r = rt2r*s rt1r = rt2r end if rt1i = zero rt2i = zero end if end if ! look for two consecutive small subdiagonal elements. do m = i - 2, l, -1 ! determine the effect of starting the double-shift qr ! iteration at row m, and see if this would make h(m,m-1) ! negligible. (the following uses scaling to avoid ! overflows and most underflows.) h21s = h( m+1, m ) s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s ) h21s = h( m+1, m ) / s v( 1_${ik}$ ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*( ( h( m, m )-rt2r ) / s ) - & rt1i*( rt2i / s ) v( 2_${ik}$ ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r ) v( 3_${ik}$ ) = h21s*h( m+2, m+1 ) s = abs( v( 1_${ik}$ ) ) + abs( v( 2_${ik}$ ) ) + abs( v( 3_${ik}$ ) ) v( 1_${ik}$ ) = v( 1_${ik}$ ) / s v( 2_${ik}$ ) = v( 2_${ik}$ ) / s v( 3_${ik}$ ) = v( 3_${ik}$ ) / s if( m==l )go to 60 if( abs( h( m, m-1 ) )*( abs( v( 2_${ik}$ ) )+abs( v( 3_${ik}$ ) ) )<=ulp*abs( v( 1_${ik}$ ) )*( abs( & h( m-1, m-1 ) )+abs( h( m,m ) )+abs( h( m+1, m+1 ) ) ) )go to 60 end do 60 continue ! double-shift qr step loop_130: do k = m, i - 1 ! the first iteration of this loop determines a reflection g ! from the vector v and applies it from left and right to h, ! thus creating a nonzero bulge below the subdiagonal. ! each subsequent iteration determines a reflection g to ! restore the hessenberg form in the (k-1)th column, and thus ! chases the bulge one step toward the bottom of the active ! submatrix. nr is the order of g. nr = min( 3_${ik}$, i-k+1 ) if( k>m )call stdlib${ii}$_scopy( nr, h( k, k-1 ), 1_${ik}$, v, 1_${ik}$ ) call stdlib${ii}$_slarfg( nr, v( 1_${ik}$ ), v( 2_${ik}$ ), 1_${ik}$, t1 ) if( k>m ) then h( k, k-1 ) = v( 1_${ik}$ ) h( k+1, k-1 ) = zero if( k<i-1 )h( k+2, k-1 ) = zero else if( m>l ) then ! ==== use the following instead of ! . h( k, k-1 ) = -h( k, k-1 ) to ! . avoid a bug when v(2) and v(3) ! . underflow. ==== h( k, k-1 ) = h( k, k-1 )*( one-t1 ) end if v2 = v( 2_${ik}$ ) t2 = t1*v2 if( nr==3_${ik}$ ) then v3 = v( 3_${ik}$ ) t3 = t1*v3 ! apply g from the left to transform the rows of the matrix ! in columns k to i2. do j = k, i2 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j ) h( k, j ) = h( k, j ) - sum*t1 h( k+1, j ) = h( k+1, j ) - sum*t2 h( k+2, j ) = h( k+2, j ) - sum*t3 end do ! apply g from the right to transform the columns of the ! matrix in rows i1 to min(k+3,i). do j = i1, min( k+3, i ) sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 ) h( j, k ) = h( j, k ) - sum*t1 h( j, k+1 ) = h( j, k+1 ) - sum*t2 h( j, k+2 ) = h( j, k+2 ) - sum*t3 end do if( wantz ) then ! accumulate transformations in the matrix z do j = iloz, ihiz sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 ) z( j, k ) = z( j, k ) - sum*t1 z( j, k+1 ) = z( j, k+1 ) - sum*t2 z( j, k+2 ) = z( j, k+2 ) - sum*t3 end do end if else if( nr==2_${ik}$ ) then ! apply g from the left to transform the rows of the matrix ! in columns k to i2. do j = k, i2 sum = h( k, j ) + v2*h( k+1, j ) h( k, j ) = h( k, j ) - sum*t1 h( k+1, j ) = h( k+1, j ) - sum*t2 end do ! apply g from the right to transform the columns of the ! matrix in rows i1 to min(k+3,i). do j = i1, i sum = h( j, k ) + v2*h( j, k+1 ) h( j, k ) = h( j, k ) - sum*t1 h( j, k+1 ) = h( j, k+1 ) - sum*t2 end do if( wantz ) then ! accumulate transformations in the matrix z do j = iloz, ihiz sum = z( j, k ) + v2*z( j, k+1 ) z( j, k ) = z( j, k ) - sum*t1 z( j, k+1 ) = z( j, k+1 ) - sum*t2 end do end if end if end do loop_130 end do loop_140 ! failure to converge in remaining number of iterations info = i return 150 continue if( l==i ) then ! h(i,i-1) is negligible: one eigenvalue has converged. wr( i ) = h( i, i ) wi( i ) = zero else if( l==i-1 ) then ! h(i-1,i-2) is negligible: a pair of eigenvalues have converged. ! transform the 2-by-2 submatrix to standard schur form, ! and compute and store the eigenvalues. call stdlib${ii}$_slanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),h( i, i ), wr( i-1 ), & wi( i-1 ), wr( i ), wi( i ),cs, sn ) if( wantt ) then ! apply the transformation to the rest of h. if( i2>i )call stdlib${ii}$_srot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,cs, sn ) call stdlib${ii}$_srot( i-i1-1, h( i1, i-1 ), 1_${ik}$, h( i1, i ), 1_${ik}$, cs, sn ) end if if( wantz ) then ! apply the transformation to z. call stdlib${ii}$_srot( nz, z( iloz, i-1 ), 1_${ik}$, z( iloz, i ), 1_${ik}$, cs, sn ) end if end if ! reset deflation counter kdefl = 0_${ik}$ ! return to start of the main loop with new value of i. i = l - 1_${ik}$ go to 20 160 continue return end subroutine stdlib${ii}$_slahqr pure module subroutine stdlib${ii}$_dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,iloz, ihiz, z, ldz, & !! DLAHQR is an auxiliary routine called by DHSEQR to update the !! eigenvalues and Schur decomposition already computed by DHSEQR, by !! dealing with the Hessenberg submatrix in rows and columns ILO to !! IHI. info ) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n integer(${ik}$), intent(out) :: info logical(lk), intent(in) :: wantt, wantz ! Array Arguments real(dp), intent(inout) :: h(ldh,*), z(ldz,*) real(dp), intent(out) :: wi(*), wr(*) ! ========================================================= ! Parameters real(dp), parameter :: dat1 = 3.0_dp/4.0_dp real(dp), parameter :: dat2 = -0.4375_dp integer(${ik}$), parameter :: kexsh = 10_${ik}$ ! Local Scalars real(dp) :: aa, ab, ba, bb, cs, det, h11, h12, h21, h21s, h22, rt1i, rt1r, rt2i, rt2r, & rtdisc, s, safmax, safmin, smlnum, sn, sum, t1, t2, t3, tr, tst, ulp, v2, v3 integer(${ik}$) :: i, i1, i2, its, itmax, j, k, l, m, nh, nr, nz, kdefl ! Local Arrays real(dp) :: v(3_${ik}$) ! Intrinsic Functions ! Executable Statements info = 0_${ik}$ ! quick return if possible if( n==0 )return if( ilo==ihi ) then wr( ilo ) = h( ilo, ilo ) wi( ilo ) = zero return end if ! ==== clear out the trash ==== do j = ilo, ihi - 3 h( j+2, j ) = zero h( j+3, j ) = zero end do if( ilo<=ihi-2 )h( ihi, ihi-2 ) = zero nh = ihi - ilo + 1_${ik}$ nz = ihiz - iloz + 1_${ik}$ ! set machine-dependent constants for the stopping criterion. safmin = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) safmax = one / safmin call stdlib${ii}$_dlabad( safmin, safmax ) ulp = stdlib${ii}$_dlamch( 'PRECISION' ) smlnum = safmin*( real( nh,KIND=dp) / ulp ) ! i1 and i2 are the indices of the first row and last column of h ! to which transformations must be applied. if eigenvalues only are ! being computed, i1 and i2 are set inside the main loop. if( wantt ) then i1 = 1_${ik}$ i2 = n end if ! itmax is the total number of qr iterations allowed. itmax = 30_${ik}$ * max( 10_${ik}$, nh ) ! kdefl counts the number of iterations since a deflation kdefl = 0_${ik}$ ! the main loop begins here. i is the loop index and decreases from ! ihi to ilo in steps of 1 or 2. each iteration of the loop works ! with the active submatrix in rows and columns l to i. ! eigenvalues i+1 to ihi have already converged. either l = ilo or ! h(l,l-1) is negligible so that the matrix splits. i = ihi 20 continue l = ilo if( i<ilo )go to 160 ! perform qr iterations on rows and columns ilo to i until a ! submatrix of order 1 or 2 splits off at the bottom because a ! subdiagonal element has become negligible. loop_140: do its = 0, itmax ! look for a single small subdiagonal element. do k = i, l + 1, -1 if( abs( h( k, k-1 ) )<=smlnum )go to 40 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) ) if( tst==zero ) then if( k-2>=ilo )tst = tst + abs( h( k-1, k-2 ) ) if( k+1<=ihi )tst = tst + abs( h( k+1, k ) ) end if ! ==== the following is a conservative small subdiagonal ! . deflation criterion due to ahues ! . 1997). it has better mathematical foundation and ! . improves accuracy in some cases. ==== if( abs( h( k, k-1 ) )<=ulp*tst ) then ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) ) ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) ) aa = max( abs( h( k, k ) ),abs( h( k-1, k-1 )-h( k, k ) ) ) bb = min( abs( h( k, k ) ),abs( h( k-1, k-1 )-h( k, k ) ) ) s = aa + ab if( ba*( ab / s )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 40 end if end do 40 continue l = k if( l>ilo ) then ! h(l,l-1) is negligible h( l, l-1 ) = zero end if ! exit from loop if a submatrix of order 1 or 2 has split off. if( l>=i-1 )go to 150 kdefl = kdefl + 1_${ik}$ ! now the active submatrix is in rows and columns l to i. if ! eigenvalues only are being computed, only the active submatrix ! need be transformed. if( .not.wantt ) then i1 = l i2 = i end if if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then ! exceptional shift. s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) ) h11 = dat1*s + h( i, i ) h12 = dat2*s h21 = s h22 = h11 else if( mod(kdefl,kexsh)==0_${ik}$ ) then ! exceptional shift. s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) ) h11 = dat1*s + h( l, l ) h12 = dat2*s h21 = s h22 = h11 else ! prepare to use francis' double shift ! (i.e. 2nd degree generalized rayleigh quotient) h11 = h( i-1, i-1 ) h21 = h( i, i-1 ) h12 = h( i-1, i ) h22 = h( i, i ) end if s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 ) if( s==zero ) then rt1r = zero rt1i = zero rt2r = zero rt2i = zero else h11 = h11 / s h21 = h21 / s h12 = h12 / s h22 = h22 / s tr = ( h11+h22 ) / two det = ( h11-tr )*( h22-tr ) - h12*h21 rtdisc = sqrt( abs( det ) ) if( det>=zero ) then ! ==== complex conjugate shifts ==== rt1r = tr*s rt2r = rt1r rt1i = rtdisc*s rt2i = -rt1i else ! ==== realshifts (use only one of them,KIND=dp) ==== rt1r = tr + rtdisc rt2r = tr - rtdisc if( abs( rt1r-h22 )<=abs( rt2r-h22 ) ) then rt1r = rt1r*s rt2r = rt1r else rt2r = rt2r*s rt1r = rt2r end if rt1i = zero rt2i = zero end if end if ! look for two consecutive small subdiagonal elements. do m = i - 2, l, -1 ! determine the effect of starting the double-shift qr ! iteration at row m, and see if this would make h(m,m-1) ! negligible. (the following uses scaling to avoid ! overflows and most underflows.) h21s = h( m+1, m ) s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s ) h21s = h( m+1, m ) / s v( 1_${ik}$ ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*( ( h( m, m )-rt2r ) / s ) - & rt1i*( rt2i / s ) v( 2_${ik}$ ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r ) v( 3_${ik}$ ) = h21s*h( m+2, m+1 ) s = abs( v( 1_${ik}$ ) ) + abs( v( 2_${ik}$ ) ) + abs( v( 3_${ik}$ ) ) v( 1_${ik}$ ) = v( 1_${ik}$ ) / s v( 2_${ik}$ ) = v( 2_${ik}$ ) / s v( 3_${ik}$ ) = v( 3_${ik}$ ) / s if( m==l )go to 60 if( abs( h( m, m-1 ) )*( abs( v( 2_${ik}$ ) )+abs( v( 3_${ik}$ ) ) )<=ulp*abs( v( 1_${ik}$ ) )*( abs( & h( m-1, m-1 ) )+abs( h( m,m ) )+abs( h( m+1, m+1 ) ) ) )go to 60 end do 60 continue ! double-shift qr step loop_130: do k = m, i - 1 ! the first iteration of this loop determines a reflection g ! from the vector v and applies it from left and right to h, ! thus creating a nonzero bulge below the subdiagonal. ! each subsequent iteration determines a reflection g to ! restore the hessenberg form in the (k-1)th column, and thus ! chases the bulge one step toward the bottom of the active ! submatrix. nr is the order of g. nr = min( 3_${ik}$, i-k+1 ) if( k>m )call stdlib${ii}$_dcopy( nr, h( k, k-1 ), 1_${ik}$, v, 1_${ik}$ ) call stdlib${ii}$_dlarfg( nr, v( 1_${ik}$ ), v( 2_${ik}$ ), 1_${ik}$, t1 ) if( k>m ) then h( k, k-1 ) = v( 1_${ik}$ ) h( k+1, k-1 ) = zero if( k<i-1 )h( k+2, k-1 ) = zero else if( m>l ) then ! ==== use the following instead of ! . h( k, k-1 ) = -h( k, k-1 ) to ! . avoid a bug when v(2) and v(3) ! . underflow. ==== h( k, k-1 ) = h( k, k-1 )*( one-t1 ) end if v2 = v( 2_${ik}$ ) t2 = t1*v2 if( nr==3_${ik}$ ) then v3 = v( 3_${ik}$ ) t3 = t1*v3 ! apply g from the left to transform the rows of the matrix ! in columns k to i2. do j = k, i2 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j ) h( k, j ) = h( k, j ) - sum*t1 h( k+1, j ) = h( k+1, j ) - sum*t2 h( k+2, j ) = h( k+2, j ) - sum*t3 end do ! apply g from the right to transform the columns of the ! matrix in rows i1 to min(k+3,i). do j = i1, min( k+3, i ) sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 ) h( j, k ) = h( j, k ) - sum*t1 h( j, k+1 ) = h( j, k+1 ) - sum*t2 h( j, k+2 ) = h( j, k+2 ) - sum*t3 end do if( wantz ) then ! accumulate transformations in the matrix z do j = iloz, ihiz sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 ) z( j, k ) = z( j, k ) - sum*t1 z( j, k+1 ) = z( j, k+1 ) - sum*t2 z( j, k+2 ) = z( j, k+2 ) - sum*t3 end do end if else if( nr==2_${ik}$ ) then ! apply g from the left to transform the rows of the matrix ! in columns k to i2. do j = k, i2 sum = h( k, j ) + v2*h( k+1, j ) h( k, j ) = h( k, j ) - sum*t1 h( k+1, j ) = h( k+1, j ) - sum*t2 end do ! apply g from the right to transform the columns of the ! matrix in rows i1 to min(k+3,i). do j = i1, i sum = h( j, k ) + v2*h( j, k+1 ) h( j, k ) = h( j, k ) - sum*t1 h( j, k+1 ) = h( j, k+1 ) - sum*t2 end do if( wantz ) then ! accumulate transformations in the matrix z do j = iloz, ihiz sum = z( j, k ) + v2*z( j, k+1 ) z( j, k ) = z( j, k ) - sum*t1 z( j, k+1 ) = z( j, k+1 ) - sum*t2 end do end if end if end do loop_130 end do loop_140 ! failure to converge in remaining number of iterations info = i return 150 continue if( l==i ) then ! h(i,i-1) is negligible: one eigenvalue has converged. wr( i ) = h( i, i ) wi( i ) = zero else if( l==i-1 ) then ! h(i-1,i-2) is negligible: a pair of eigenvalues have converged. ! transform the 2-by-2 submatrix to standard schur form, ! and compute and store the eigenvalues. call stdlib${ii}$_dlanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),h( i, i ), wr( i-1 ), & wi( i-1 ), wr( i ), wi( i ),cs, sn ) if( wantt ) then ! apply the transformation to the rest of h. if( i2>i )call stdlib${ii}$_drot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,cs, sn ) call stdlib${ii}$_drot( i-i1-1, h( i1, i-1 ), 1_${ik}$, h( i1, i ), 1_${ik}$, cs, sn ) end if if( wantz ) then ! apply the transformation to z. call stdlib${ii}$_drot( nz, z( iloz, i-1 ), 1_${ik}$, z( iloz, i ), 1_${ik}$, cs, sn ) end if end if ! reset deflation counter kdefl = 0_${ik}$ ! return to start of the main loop with new value of i. i = l - 1_${ik}$ go to 20 160 continue return end subroutine stdlib${ii}$_dlahqr #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$lahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,iloz, ihiz, z, ldz, & !! DLAHQR: is an auxiliary routine called by DHSEQR to update the !! eigenvalues and Schur decomposition already computed by DHSEQR, by !! dealing with the Hessenberg submatrix in rows and columns ILO to !! IHI. info ) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n integer(${ik}$), intent(out) :: info logical(lk), intent(in) :: wantt, wantz ! Array Arguments real(${rk}$), intent(inout) :: h(ldh,*), z(ldz,*) real(${rk}$), intent(out) :: wi(*), wr(*) ! ========================================================= ! Parameters real(${rk}$), parameter :: dat1 = 3.0_${rk}$/4.0_${rk}$ real(${rk}$), parameter :: dat2 = -0.4375_${rk}$ integer(${ik}$), parameter :: kexsh = 10_${ik}$ ! Local Scalars real(${rk}$) :: aa, ab, ba, bb, cs, det, h11, h12, h21, h21s, h22, rt1i, rt1r, rt2i, rt2r, & rtdisc, s, safmax, safmin, smlnum, sn, sum, t1, t2, t3, tr, tst, ulp, v2, v3 integer(${ik}$) :: i, i1, i2, its, itmax, j, k, l, m, nh, nr, nz, kdefl ! Local Arrays real(${rk}$) :: v(3_${ik}$) ! Intrinsic Functions ! Executable Statements info = 0_${ik}$ ! quick return if possible if( n==0 )return if( ilo==ihi ) then wr( ilo ) = h( ilo, ilo ) wi( ilo ) = zero return end if ! ==== clear out the trash ==== do j = ilo, ihi - 3 h( j+2, j ) = zero h( j+3, j ) = zero end do if( ilo<=ihi-2 )h( ihi, ihi-2 ) = zero nh = ihi - ilo + 1_${ik}$ nz = ihiz - iloz + 1_${ik}$ ! set machine-dependent constants for the stopping criterion. safmin = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) safmax = one / safmin call stdlib${ii}$_${ri}$labad( safmin, safmax ) ulp = stdlib${ii}$_${ri}$lamch( 'PRECISION' ) smlnum = safmin*( real( nh,KIND=${rk}$) / ulp ) ! i1 and i2 are the indices of the first row and last column of h ! to which transformations must be applied. if eigenvalues only are ! being computed, i1 and i2 are set inside the main loop. if( wantt ) then i1 = 1_${ik}$ i2 = n end if ! itmax is the total number of qr iterations allowed. itmax = 30_${ik}$ * max( 10_${ik}$, nh ) ! kdefl counts the number of iterations since a deflation kdefl = 0_${ik}$ ! the main loop begins here. i is the loop index and decreases from ! ihi to ilo in steps of 1 or 2. each iteration of the loop works ! with the active submatrix in rows and columns l to i. ! eigenvalues i+1 to ihi have already converged. either l = ilo or ! h(l,l-1) is negligible so that the matrix splits. i = ihi 20 continue l = ilo if( i<ilo )go to 160 ! perform qr iterations on rows and columns ilo to i until a ! submatrix of order 1 or 2 splits off at the bottom because a ! subdiagonal element has become negligible. loop_140: do its = 0, itmax ! look for a single small subdiagonal element. do k = i, l + 1, -1 if( abs( h( k, k-1 ) )<=smlnum )go to 40 tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) ) if( tst==zero ) then if( k-2>=ilo )tst = tst + abs( h( k-1, k-2 ) ) if( k+1<=ihi )tst = tst + abs( h( k+1, k ) ) end if ! ==== the following is a conservative small subdiagonal ! . deflation criterion due to ahues ! . 1997). it has better mathematical foundation and ! . improves accuracy in some cases. ==== if( abs( h( k, k-1 ) )<=ulp*tst ) then ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) ) ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) ) aa = max( abs( h( k, k ) ),abs( h( k-1, k-1 )-h( k, k ) ) ) bb = min( abs( h( k, k ) ),abs( h( k-1, k-1 )-h( k, k ) ) ) s = aa + ab if( ba*( ab / s )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 40 end if end do 40 continue l = k if( l>ilo ) then ! h(l,l-1) is negligible h( l, l-1 ) = zero end if ! exit from loop if a submatrix of order 1 or 2 has split off. if( l>=i-1 )go to 150 kdefl = kdefl + 1_${ik}$ ! now the active submatrix is in rows and columns l to i. if ! eigenvalues only are being computed, only the active submatrix ! need be transformed. if( .not.wantt ) then i1 = l i2 = i end if if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then ! exceptional shift. s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) ) h11 = dat1*s + h( i, i ) h12 = dat2*s h21 = s h22 = h11 else if( mod(kdefl,kexsh)==0_${ik}$ ) then ! exceptional shift. s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) ) h11 = dat1*s + h( l, l ) h12 = dat2*s h21 = s h22 = h11 else ! prepare to use francis' double shift ! (i.e. 2nd degree generalized rayleigh quotient) h11 = h( i-1, i-1 ) h21 = h( i, i-1 ) h12 = h( i-1, i ) h22 = h( i, i ) end if s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 ) if( s==zero ) then rt1r = zero rt1i = zero rt2r = zero rt2i = zero else h11 = h11 / s h21 = h21 / s h12 = h12 / s h22 = h22 / s tr = ( h11+h22 ) / two det = ( h11-tr )*( h22-tr ) - h12*h21 rtdisc = sqrt( abs( det ) ) if( det>=zero ) then ! ==== complex conjugate shifts ==== rt1r = tr*s rt2r = rt1r rt1i = rtdisc*s rt2i = -rt1i else ! ==== realshifts (use only one of them,KIND=${rk}$) ==== rt1r = tr + rtdisc rt2r = tr - rtdisc if( abs( rt1r-h22 )<=abs( rt2r-h22 ) ) then rt1r = rt1r*s rt2r = rt1r else rt2r = rt2r*s rt1r = rt2r end if rt1i = zero rt2i = zero end if end if ! look for two consecutive small subdiagonal elements. do m = i - 2, l, -1 ! determine the effect of starting the double-shift qr ! iteration at row m, and see if this would make h(m,m-1) ! negligible. (the following uses scaling to avoid ! overflows and most underflows.) h21s = h( m+1, m ) s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s ) h21s = h( m+1, m ) / s v( 1_${ik}$ ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*( ( h( m, m )-rt2r ) / s ) - & rt1i*( rt2i / s ) v( 2_${ik}$ ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r ) v( 3_${ik}$ ) = h21s*h( m+2, m+1 ) s = abs( v( 1_${ik}$ ) ) + abs( v( 2_${ik}$ ) ) + abs( v( 3_${ik}$ ) ) v( 1_${ik}$ ) = v( 1_${ik}$ ) / s v( 2_${ik}$ ) = v( 2_${ik}$ ) / s v( 3_${ik}$ ) = v( 3_${ik}$ ) / s if( m==l )go to 60 if( abs( h( m, m-1 ) )*( abs( v( 2_${ik}$ ) )+abs( v( 3_${ik}$ ) ) )<=ulp*abs( v( 1_${ik}$ ) )*( abs( & h( m-1, m-1 ) )+abs( h( m,m ) )+abs( h( m+1, m+1 ) ) ) )go to 60 end do 60 continue ! double-shift qr step loop_130: do k = m, i - 1 ! the first iteration of this loop determines a reflection g ! from the vector v and applies it from left and right to h, ! thus creating a nonzero bulge below the subdiagonal. ! each subsequent iteration determines a reflection g to ! restore the hessenberg form in the (k-1)th column, and thus ! chases the bulge one step toward the bottom of the active ! submatrix. nr is the order of g. nr = min( 3_${ik}$, i-k+1 ) if( k>m )call stdlib${ii}$_${ri}$copy( nr, h( k, k-1 ), 1_${ik}$, v, 1_${ik}$ ) call stdlib${ii}$_${ri}$larfg( nr, v( 1_${ik}$ ), v( 2_${ik}$ ), 1_${ik}$, t1 ) if( k>m ) then h( k, k-1 ) = v( 1_${ik}$ ) h( k+1, k-1 ) = zero if( k<i-1 )h( k+2, k-1 ) = zero else if( m>l ) then ! ==== use the following instead of ! . h( k, k-1 ) = -h( k, k-1 ) to ! . avoid a bug when v(2) and v(3) ! . underflow. ==== h( k, k-1 ) = h( k, k-1 )*( one-t1 ) end if v2 = v( 2_${ik}$ ) t2 = t1*v2 if( nr==3_${ik}$ ) then v3 = v( 3_${ik}$ ) t3 = t1*v3 ! apply g from the left to transform the rows of the matrix ! in columns k to i2. do j = k, i2 sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j ) h( k, j ) = h( k, j ) - sum*t1 h( k+1, j ) = h( k+1, j ) - sum*t2 h( k+2, j ) = h( k+2, j ) - sum*t3 end do ! apply g from the right to transform the columns of the ! matrix in rows i1 to min(k+3,i). do j = i1, min( k+3, i ) sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 ) h( j, k ) = h( j, k ) - sum*t1 h( j, k+1 ) = h( j, k+1 ) - sum*t2 h( j, k+2 ) = h( j, k+2 ) - sum*t3 end do if( wantz ) then ! accumulate transformations in the matrix z do j = iloz, ihiz sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 ) z( j, k ) = z( j, k ) - sum*t1 z( j, k+1 ) = z( j, k+1 ) - sum*t2 z( j, k+2 ) = z( j, k+2 ) - sum*t3 end do end if else if( nr==2_${ik}$ ) then ! apply g from the left to transform the rows of the matrix ! in columns k to i2. do j = k, i2 sum = h( k, j ) + v2*h( k+1, j ) h( k, j ) = h( k, j ) - sum*t1 h( k+1, j ) = h( k+1, j ) - sum*t2 end do ! apply g from the right to transform the columns of the ! matrix in rows i1 to min(k+3,i). do j = i1, i sum = h( j, k ) + v2*h( j, k+1 ) h( j, k ) = h( j, k ) - sum*t1 h( j, k+1 ) = h( j, k+1 ) - sum*t2 end do if( wantz ) then ! accumulate transformations in the matrix z do j = iloz, ihiz sum = z( j, k ) + v2*z( j, k+1 ) z( j, k ) = z( j, k ) - sum*t1 z( j, k+1 ) = z( j, k+1 ) - sum*t2 end do end if end if end do loop_130 end do loop_140 ! failure to converge in remaining number of iterations info = i return 150 continue if( l==i ) then ! h(i,i-1) is negligible: one eigenvalue has converged. wr( i ) = h( i, i ) wi( i ) = zero else if( l==i-1 ) then ! h(i-1,i-2) is negligible: a pair of eigenvalues have converged. ! transform the 2-by-2 submatrix to standard schur form, ! and compute and store the eigenvalues. call stdlib${ii}$_${ri}$lanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),h( i, i ), wr( i-1 ), & wi( i-1 ), wr( i ), wi( i ),cs, sn ) if( wantt ) then ! apply the transformation to the rest of h. if( i2>i )call stdlib${ii}$_${ri}$rot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,cs, sn ) call stdlib${ii}$_${ri}$rot( i-i1-1, h( i1, i-1 ), 1_${ik}$, h( i1, i ), 1_${ik}$, cs, sn ) end if if( wantz ) then ! apply the transformation to z. call stdlib${ii}$_${ri}$rot( nz, z( iloz, i-1 ), 1_${ik}$, z( iloz, i ), 1_${ik}$, cs, sn ) end if end if ! reset deflation counter kdefl = 0_${ik}$ ! return to start of the main loop with new value of i. i = l - 1_${ik}$ go to 20 160 continue return end subroutine stdlib${ii}$_${ri}$lahqr #:endif #:endfor pure module subroutine stdlib${ii}$_clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,ihiz, z, ldz, info & !! CLAHQR is an auxiliary routine called by CHSEQR to update the !! eigenvalues and Schur decomposition already computed by CHSEQR, by !! dealing with the Hessenberg submatrix in rows and columns ILO to !! IHI. ) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: ihi, ihiz, ilo, iloz, ldh, ldz, n integer(${ik}$), intent(out) :: info logical(lk), intent(in) :: wantt, wantz ! Array Arguments complex(sp), intent(inout) :: h(ldh,*), z(ldz,*) complex(sp), intent(out) :: w(*) ! ========================================================= ! Parameters real(sp), parameter :: dat1 = 3.0_sp/4.0_sp integer(${ik}$), parameter :: kexsh = 10_${ik}$ ! Local Scalars complex(sp) :: cdum, h11, h11s, h22, sc, sum, t, t1, temp, u, v2, x, y real(sp) :: aa, ab, ba, bb, h10, h21, rtemp, s, safmax, safmin, smlnum, sx, t2, tst, & ulp integer(${ik}$) :: i, i1, i2, its, itmax, j, jhi, jlo, k, l, m, nh, nz, kdefl ! Local Arrays complex(sp) :: v(2_${ik}$) ! Statement Functions real(sp) :: cabs1 ! Intrinsic Functions ! Statement Function Definitions cabs1( cdum ) = abs( real( cdum,KIND=sp) ) + abs( aimag( cdum ) ) ! Executable Statements info = 0_${ik}$ ! quick return if possible if( n==0 )return if( ilo==ihi ) then w( ilo ) = h( ilo, ilo ) return end if ! ==== clear out the trash ==== do j = ilo, ihi - 3 h( j+2, j ) = czero h( j+3, j ) = czero end do if( ilo<=ihi-2 )h( ihi, ihi-2 ) = czero ! ==== ensure that subdiagonal entries are real ==== if( wantt ) then jlo = 1_${ik}$ jhi = n else jlo = ilo jhi = ihi end if do i = ilo + 1, ihi if( aimag( h( i, i-1 ) )/=zero ) then ! ==== the following redundant normalization ! . avoids problems with both gradual and ! . sudden underflow in abs(h(i,i-1)) ==== sc = h( i, i-1 ) / cabs1( h( i, i-1 ) ) sc = conjg( sc ) / abs( sc ) h( i, i-1 ) = abs( h( i, i-1 ) ) call stdlib${ii}$_cscal( jhi-i+1, sc, h( i, i ), ldh ) call stdlib${ii}$_cscal( min( jhi, i+1 )-jlo+1, conjg( sc ), h( jlo, i ),1_${ik}$ ) if( wantz )call stdlib${ii}$_cscal( ihiz-iloz+1, conjg( sc ), z( iloz, i ), 1_${ik}$ ) end if end do nh = ihi - ilo + 1_${ik}$ nz = ihiz - iloz + 1_${ik}$ ! set machine-dependent constants for the stopping criterion. safmin = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) safmax = one / safmin call stdlib${ii}$_slabad( safmin, safmax ) ulp = stdlib${ii}$_slamch( 'PRECISION' ) smlnum = safmin*( real( nh,KIND=sp) / ulp ) ! i1 and i2 are the indices of the first row and last column of h ! to which transformations must be applied. if eigenvalues only are ! being computed, i1 and i2 are set inside the main loop. if( wantt ) then i1 = 1_${ik}$ i2 = n end if ! itmax is the total number of qr iterations allowed. itmax = 30_${ik}$ * max( 10_${ik}$, nh ) ! kdefl counts the number of iterations since a deflation kdefl = 0_${ik}$ ! the main loop begins here. i is the loop index and decreases from ! ihi to ilo in steps of 1. each iteration of the loop works ! with the active submatrix in rows and columns l to i. ! eigenvalues i+1 to ihi have already converged. either l = ilo, or ! h(l,l-1) is negligible so that the matrix splits. i = ihi 30 continue if( i<ilo )go to 150 ! perform qr iterations on rows and columns ilo to i until a ! submatrix of order 1 splits off at the bottom because a ! subdiagonal element has become negligible. l = ilo loop_130: do its = 0, itmax ! look for a single small subdiagonal element. do k = i, l + 1, -1 if( cabs1( h( k, k-1 ) )<=smlnum )go to 50 tst = cabs1( h( k-1, k-1 ) ) + cabs1( h( k, k ) ) if( tst==czero ) then if( k-2>=ilo )tst = tst + abs( real( h( k-1, k-2 ),KIND=sp) ) if( k+1<=ihi )tst = tst + abs( real( h( k+1, k ),KIND=sp) ) end if ! ==== the following is a conservative small subdiagonal ! . deflation criterion due to ahues ! . 1997). it has better mathematical foundation and ! . improves accuracy in some examples. ==== if( abs( real( h( k, k-1 ),KIND=sp) )<=ulp*tst ) then ab = max( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) ) ba = min( cabs1( h( k, k-1 ) ), cabs1( h( k-1, k ) ) ) aa = max( cabs1( h( k, k ) ),cabs1( h( k-1, k-1 )-h( k, k ) ) ) bb = min( cabs1( h( k, k ) ),cabs1( h( k-1, k-1 )-h( k, k ) ) ) s = aa + ab if( ba*( ab / s )<=max( smlnum,ulp*( bb*( aa / s ) ) ) )go to 50 end if end do 50 continue l = k if( l>ilo ) then ! h(l,l-1) is negligible h( l, l-1 ) = czero end if ! exit from loop if a submatrix of order 1 has split off. if( l>=i )go to 140 kdefl = kdefl + 1_${ik}$ ! now the active submatrix is in rows and columns l to i. if ! eigenvalues only are being computed, only the active submatrix ! need be transformed. if( .not.wantt ) then i1 = l i2 = i end if if( mod(kdefl,2_${ik}$*kexsh)==0_${ik}$ ) then ! exceptional shift. s = dat1*abs( real( h( i, i-1 ),KIND=sp) ) t = s + h( i, i ) else if( mod(kdefl,kexsh)==0_${ik}$ ) then ! exceptional shift. s = dat1*abs( real( h( l+1, l ),KIND=sp) ) t = s + h( l, l ) else ! wilkinson's shift. t = h( i, i ) u = sqrt( h( i-1, i ) )*sqrt( h( i, i-1 ) ) s = cabs1( u ) if( s/=zero ) then x = half*( h( i-1, i-1 )-t ) sx = cabs1( x ) s = max( s, cabs1( x ) ) y = s*sqrt( ( x / s )**2_${ik}$+( u / s )**2_${ik}$ ) if( sx>zero ) then if( real( x / sx,KIND=sp)*real( y,KIND=sp)+aimag( x / sx )*aimag( y )& <zero )y = -y end if t = t - u*stdlib${ii}$_cladiv( u, ( x+y ) ) end if end if ! look for two consecutive small subdiagonal elements. do m = i - 1, l + 1, -1 ! determine the effect of starting the single-shift qr ! iteration at row m, and see if this would make h(m,m-1) ! negligible. h11 = h( m, m ) h22 = h( m+1, m+1 ) h11s = h11 - t h21 = real( h( m+1, m ),KIND=sp) s = cabs1( h11s ) + abs( h21 ) h11s = h11s / s h21