#:include "common.fypp" submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_comp2 implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_stgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alphar, alphai, & !! STGSEN reorders the generalized real Schur decomposition of a real !! matrix pair (A, B) (in terms of an orthonormal equivalence trans- !! formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues !! appears in the leading diagonal blocks of the upper quasi-triangular !! matrix A and the upper triangular B. The leading columns of Q and !! Z form orthonormal bases of the corresponding left and right eigen- !! spaces (deflating subspaces). (A, B) must be in generalized real !! Schur canonical form (as returned by SGGES), i.e. A is block upper !! triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper !! triangular. !! STGSEN also computes the generalized eigenvalues !! w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) !! of the reordered matrix pair (A, B). !! Optionally, STGSEN computes the estimates of reciprocal condition !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to !! the selected cluster and the eigenvalues outside the cluster, resp., !! and norms of "projections" onto left and right eigenspaces w.r.t. !! the selected cluster in the (1,1)-block. beta, q, ldq, z, ldz, m, pl,pr, dif, work, lwork, iwork, liwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: wantq, wantz integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(sp), intent(out) :: pl, pr ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*) real(sp), intent(out) :: alphai(*), alphar(*), beta(*), dif(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: idifjb = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, pair, swap, wantd, wantd1, wantd2, wantp integer(${ik}$) :: i, ierr, ijb, k, kase, kk, ks, liwmin, lwmin, mn2, n1, n2 real(sp) :: dscale, dsum, eps, rdscal, smlnum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -9_${ik}$ else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then info = -14_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'STGSEN', -info ) return end if ! get machine constants eps = stdlib${ii}$_slamch( 'P' ) smlnum = stdlib${ii}$_slamch( 'S' ) / eps ierr = 0_${ik}$ wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$ wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$ wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$ wantd = wantd1 .or. wantd2 ! set m to the dimension of the specified pair of deflating ! subspaces. m = 0_${ik}$ pair = .false. if( .not.lquery .or. ijob/=0_${ik}$ ) then do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$ end if else if( select( n ) )m = m + 1_${ik}$ end if end if end do end if if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 2_${ik}$*m*(n-m) ) liwmin = max( 1_${ik}$, n+6 ) else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 4_${ik}$*m*(n-m) ) liwmin = max( 1_${ik}$, 2_${ik}$*m*(n-m), n+6 ) else lwmin = max( 1_${ik}$, 4_${ik}$*n+16 ) liwmin = 1_${ik}$ end if work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -22_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -24_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'STGSEN', -info ) return else if( lquery ) then return end if ! quick return if possible. if( m==n .or. m==0_${ik}$ ) then if( wantp ) then pl = one pr = one end if if( wantd ) then dscale = zero dsum = one do i = 1, n call stdlib${ii}$_slassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) call stdlib${ii}$_slassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) end do dif( 1_${ik}$ ) = dscale*sqrt( dsum ) dif( 2_${ik}$ ) = dif( 1_${ik}$ ) end if go to 60 end if ! collect the selected blocks at the top-left corner of (a, b). ks = 0_${ik}$ pair = .false. loop_30: do k = 1, n if( pair ) then pair = .false. else swap = select( k ) if( k<n ) then if( a( k+1, k )/=zero ) then pair = .true. swap = swap .or. select( k+1 ) end if end if if( swap ) then ks = ks + 1_${ik}$ ! swap the k-th block to position ks. ! perform the reordering of diagonal blocks in (a, b) ! by orthogonal transformation matrices and update ! q and z accordingly (if requested): kk = k if( k/=ks )call stdlib${ii}$_stgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,z, ldz,& kk, ks, work, lwork, ierr ) if( ierr>0_${ik}$ ) then ! swap is rejected: exit. info = 1_${ik}$ if( wantp ) then pl = zero pr = zero end if if( wantd ) then dif( 1_${ik}$ ) = zero dif( 2_${ik}$ ) = zero end if go to 60 end if if( pair )ks = ks + 1_${ik}$ end if end if end do loop_30 if( wantp ) then ! solve generalized sylvester equation for r and l ! and compute pl and pr. n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ call stdlib${ii}$_slacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 ) call stdlib${ii}$_slacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 ) call stdlib${ii}$_stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,& i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, & iwork, ierr ) ! estimate the reciprocal of norms of "projections" onto left ! and right eigenspaces. rdscal = zero dsum = one call stdlib${ii}$_slassq( n1*n2, work, 1_${ik}$, rdscal, dsum ) pl = rdscal*sqrt( dsum ) if( pl==zero ) then pl = one else pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) ) end if rdscal = zero dsum = one call stdlib${ii}$_slassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum ) pr = rdscal*sqrt( dsum ) if( pr==zero ) then pr = one else pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) ) end if end if if( wantd ) then ! compute estimates of difu and difl. if( wantd1 ) then n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = idifjb ! frobenius norm-based difu-estimate. call stdlib${ii}$_stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(& i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) ! frobenius norm-based difl-estimate. call stdlib${ii}$_stgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),& ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) else ! compute 1-norm-based estimates of difu and difl using ! reversed communication with stdlib${ii}$_slacn2. in each step a ! generalized sylvester equation or a transposed variant ! is solved. kase = 0_${ik}$ n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ mn2 = 2_${ik}$*n1*n2 ! 1-norm-based estimate of difu. 40 continue call stdlib${ii}$_slacn2( mn2, work( mn2+1 ), work, iwork, dif( 1_${ik}$ ),kase, isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation. call stdlib${ii}$_stgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_stgsyl( 'T', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 40 end if dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ ) ! 1-norm-based estimate of difl. 50 continue call stdlib${ii}$_slacn2( mn2, work( mn2+1 ), work, iwork, dif( 2_${ik}$ ),kase, isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation. call stdlib${ii}$_stgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_stgsyl( 'T', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 50 end if dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ ) end if end if 60 continue ! compute generalized eigenvalues of reordered pair (a, b) and ! normalize the generalized schur form. pair = .false. loop_70: do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )/=zero ) then pair = .true. end if end if if( pair ) then ! compute the eigenvalue(s) at position k. work( 1_${ik}$ ) = a( k, k ) work( 2_${ik}$ ) = a( k+1, k ) work( 3_${ik}$ ) = a( k, k+1 ) work( 4_${ik}$ ) = a( k+1, k+1 ) work( 5_${ik}$ ) = b( k, k ) work( 6_${ik}$ ) = b( k+1, k ) work( 7_${ik}$ ) = b( k, k+1 ) work( 8_${ik}$ ) = b( k+1, k+1 ) call stdlib${ii}$_slag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta( k ),beta( k+1 ), & alphar( k ), alphar( k+1 ),alphai( k ) ) alphai( k+1 ) = -alphai( k ) else if( sign( one, b( k, k ) )<zero ) then ! if b(k,k) is negative, make it positive do i = 1, n a( k, i ) = -a( k, i ) b( k, i ) = -b( k, i ) if( wantq ) q( i, k ) = -q( i, k ) end do end if alphar( k ) = a( k, k ) alphai( k ) = zero beta( k ) = b( k, k ) end if end if end do loop_70 work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_stgsen pure module subroutine stdlib${ii}$_dtgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alphar, alphai, & !! DTGSEN reorders the generalized real Schur decomposition of a real !! matrix pair (A, B) (in terms of an orthonormal equivalence trans- !! formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues !! appears in the leading diagonal blocks of the upper quasi-triangular !! matrix A and the upper triangular B. The leading columns of Q and !! Z form orthonormal bases of the corresponding left and right eigen- !! spaces (deflating subspaces). (A, B) must be in generalized real !! Schur canonical form (as returned by DGGES), i.e. A is block upper !! triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper !! triangular. !! DTGSEN also computes the generalized eigenvalues !! w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) !! of the reordered matrix pair (A, B). !! Optionally, DTGSEN computes the estimates of reciprocal condition !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to !! the selected cluster and the eigenvalues outside the cluster, resp., !! and norms of "projections" onto left and right eigenspaces w.r.t. !! the selected cluster in the (1,1)-block. beta, q, ldq, z, ldz, m, pl,pr, dif, work, lwork, iwork, liwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: wantq, wantz integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(dp), intent(out) :: pl, pr ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*) real(dp), intent(out) :: alphai(*), alphar(*), beta(*), dif(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: idifjb = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, pair, swap, wantd, wantd1, wantd2, wantp integer(${ik}$) :: i, ierr, ijb, k, kase, kk, ks, liwmin, lwmin, mn2, n1, n2 real(dp) :: dscale, dsum, eps, rdscal, smlnum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -9_${ik}$ else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then info = -14_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTGSEN', -info ) return end if ! get machine constants eps = stdlib${ii}$_dlamch( 'P' ) smlnum = stdlib${ii}$_dlamch( 'S' ) / eps ierr = 0_${ik}$ wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$ wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$ wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$ wantd = wantd1 .or. wantd2 ! set m to the dimension of the specified pair of deflating ! subspaces. m = 0_${ik}$ pair = .false. if( .not.lquery .or. ijob/=0_${ik}$ ) then do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$ end if else if( select( n ) )m = m + 1_${ik}$ end if end if end do end if if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 2_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, n+6 ) else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 4_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+6 ) else lwmin = max( 1_${ik}$, 4_${ik}$*n+16 ) liwmin = 1_${ik}$ end if work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -22_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -24_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTGSEN', -info ) return else if( lquery ) then return end if ! quick return if possible. if( m==n .or. m==0_${ik}$ ) then if( wantp ) then pl = one pr = one end if if( wantd ) then dscale = zero dsum = one do i = 1, n call stdlib${ii}$_dlassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) call stdlib${ii}$_dlassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) end do dif( 1_${ik}$ ) = dscale*sqrt( dsum ) dif( 2_${ik}$ ) = dif( 1_${ik}$ ) end if go to 60 end if ! collect the selected blocks at the top-left corner of (a, b). ks = 0_${ik}$ pair = .false. loop_30: do k = 1, n if( pair ) then pair = .false. else swap = select( k ) if( k<n ) then if( a( k+1, k )/=zero ) then pair = .true. swap = swap .or. select( k+1 ) end if end if if( swap ) then ks = ks + 1_${ik}$ ! swap the k-th block to position ks. ! perform the reordering of diagonal blocks in (a, b) ! by orthogonal transformation matrices and update ! q and z accordingly (if requested): kk = k if( k/=ks )call stdlib${ii}$_dtgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,z, ldz,& kk, ks, work, lwork, ierr ) if( ierr>0_${ik}$ ) then ! swap is rejected: exit. info = 1_${ik}$ if( wantp ) then pl = zero pr = zero end if if( wantd ) then dif( 1_${ik}$ ) = zero dif( 2_${ik}$ ) = zero end if go to 60 end if if( pair )ks = ks + 1_${ik}$ end if end if end do loop_30 if( wantp ) then ! solve generalized sylvester equation for r and l ! and compute pl and pr. n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ call stdlib${ii}$_dlacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 ) call stdlib${ii}$_dlacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 ) call stdlib${ii}$_dtgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,& i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, & iwork, ierr ) ! estimate the reciprocal of norms of "projections" onto left ! and right eigenspaces. rdscal = zero dsum = one call stdlib${ii}$_dlassq( n1*n2, work, 1_${ik}$, rdscal, dsum ) pl = rdscal*sqrt( dsum ) if( pl==zero ) then pl = one else pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) ) end if rdscal = zero dsum = one call stdlib${ii}$_dlassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum ) pr = rdscal*sqrt( dsum ) if( pr==zero ) then pr = one else pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) ) end if end if if( wantd ) then ! compute estimates of difu and difl. if( wantd1 ) then n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = idifjb ! frobenius norm-based difu-estimate. call stdlib${ii}$_dtgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(& i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) ! frobenius norm-based difl-estimate. call stdlib${ii}$_dtgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),& ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) else ! compute 1-norm-based estimates of difu and difl using ! reversed communication with stdlib${ii}$_dlacn2. in each step a ! generalized sylvester equation or a transposed variant ! is solved. kase = 0_${ik}$ n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ mn2 = 2_${ik}$*n1*n2 ! 1-norm-based estimate of difu. 40 continue call stdlib${ii}$_dlacn2( mn2, work( mn2+1 ), work, iwork, dif( 1_${ik}$ ),kase, isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation. call stdlib${ii}$_dtgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_dtgsyl( 'T', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 40 end if dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ ) ! 1-norm-based estimate of difl. 50 continue call stdlib${ii}$_dlacn2( mn2, work( mn2+1 ), work, iwork, dif( 2_${ik}$ ),kase, isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation. call stdlib${ii}$_dtgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_dtgsyl( 'T', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 50 end if dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ ) end if end if 60 continue ! compute generalized eigenvalues of reordered pair (a, b) and ! normalize the generalized schur form. pair = .false. loop_80: do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )/=zero ) then pair = .true. end if end if if( pair ) then ! compute the eigenvalue(s) at position k. work( 1_${ik}$ ) = a( k, k ) work( 2_${ik}$ ) = a( k+1, k ) work( 3_${ik}$ ) = a( k, k+1 ) work( 4_${ik}$ ) = a( k+1, k+1 ) work( 5_${ik}$ ) = b( k, k ) work( 6_${ik}$ ) = b( k+1, k ) work( 7_${ik}$ ) = b( k, k+1 ) work( 8_${ik}$ ) = b( k+1, k+1 ) call stdlib${ii}$_dlag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta( k ),beta( k+1 ), & alphar( k ), alphar( k+1 ),alphai( k ) ) alphai( k+1 ) = -alphai( k ) else if( sign( one, b( k, k ) )<zero ) then ! if b(k,k) is negative, make it positive do i = 1, n a( k, i ) = -a( k, i ) b( k, i ) = -b( k, i ) if( wantq ) q( i, k ) = -q( i, k ) end do end if alphar( k ) = a( k, k ) alphai( k ) = zero beta( k ) = b( k, k ) end if end if end do loop_80 work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_dtgsen #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$tgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alphar, alphai, & !! DTGSEN: reorders the generalized real Schur decomposition of a real !! matrix pair (A, B) (in terms of an orthonormal equivalence trans- !! formation Q**T * (A, B) * Z), so that a selected cluster of eigenvalues !! appears in the leading diagonal blocks of the upper quasi-triangular !! matrix A and the upper triangular B. The leading columns of Q and !! Z form orthonormal bases of the corresponding left and right eigen- !! spaces (deflating subspaces). (A, B) must be in generalized real !! Schur canonical form (as returned by DGGES), i.e. A is block upper !! triangular with 1-by-1 and 2-by-2 diagonal blocks. B is upper !! triangular. !! DTGSEN also computes the generalized eigenvalues !! w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j) !! of the reordered matrix pair (A, B). !! Optionally, DTGSEN computes the estimates of reciprocal condition !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to !! the selected cluster and the eigenvalues outside the cluster, resp., !! and norms of "projections" onto left and right eigenspaces w.r.t. !! the selected cluster in the (1,1)-block. beta, q, ldq, z, ldz, m, pl,pr, dif, work, lwork, iwork, liwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: wantq, wantz integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(${rk}$), intent(out) :: pl, pr ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*) real(${rk}$), intent(out) :: alphai(*), alphar(*), beta(*), dif(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: idifjb = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, pair, swap, wantd, wantd1, wantd2, wantp integer(${ik}$) :: i, ierr, ijb, k, kase, kk, ks, liwmin, lwmin, mn2, n1, n2 real(${rk}$) :: dscale, dsum, eps, rdscal, smlnum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -9_${ik}$ else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then info = -14_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTGSEN', -info ) return end if ! get machine constants eps = stdlib${ii}$_${ri}$lamch( 'P' ) smlnum = stdlib${ii}$_${ri}$lamch( 'S' ) / eps ierr = 0_${ik}$ wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$ wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$ wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$ wantd = wantd1 .or. wantd2 ! set m to the dimension of the specified pair of deflating ! subspaces. m = 0_${ik}$ pair = .false. if( .not.lquery .or. ijob/=0_${ik}$ ) then do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$ end if else if( select( n ) )m = m + 1_${ik}$ end if end if end do end if if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 2_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, n+6 ) else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*n+16, 4_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+6 ) else lwmin = max( 1_${ik}$, 4_${ik}$*n+16 ) liwmin = 1_${ik}$ end if work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -22_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -24_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTGSEN', -info ) return else if( lquery ) then return end if ! quick return if possible. if( m==n .or. m==0_${ik}$ ) then if( wantp ) then pl = one pr = one end if if( wantd ) then dscale = zero dsum = one do i = 1, n call stdlib${ii}$_${ri}$lassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) call stdlib${ii}$_${ri}$lassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) end do dif( 1_${ik}$ ) = dscale*sqrt( dsum ) dif( 2_${ik}$ ) = dif( 1_${ik}$ ) end if go to 60 end if ! collect the selected blocks at the top-left corner of (a, b). ks = 0_${ik}$ pair = .false. loop_30: do k = 1, n if( pair ) then pair = .false. else swap = select( k ) if( k<n ) then if( a( k+1, k )/=zero ) then pair = .true. swap = swap .or. select( k+1 ) end if end if if( swap ) then ks = ks + 1_${ik}$ ! swap the k-th block to position ks. ! perform the reordering of diagonal blocks in (a, b) ! by orthogonal transformation matrices and update ! q and z accordingly (if requested): kk = k if( k/=ks )call stdlib${ii}$_${ri}$tgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq,z, ldz,& kk, ks, work, lwork, ierr ) if( ierr>0_${ik}$ ) then ! swap is rejected: exit. info = 1_${ik}$ if( wantp ) then pl = zero pr = zero end if if( wantd ) then dif( 1_${ik}$ ) = zero dif( 2_${ik}$ ) = zero end if go to 60 end if if( pair )ks = ks + 1_${ik}$ end if end if end do loop_30 if( wantp ) then ! solve generalized sylvester equation for r and l ! and compute pl and pr. n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ call stdlib${ii}$_${ri}$lacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 ) call stdlib${ii}$_${ri}$lacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 ) call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,& i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, & iwork, ierr ) ! estimate the reciprocal of norms of "projections" onto left ! and right eigenspaces. rdscal = zero dsum = one call stdlib${ii}$_${ri}$lassq( n1*n2, work, 1_${ik}$, rdscal, dsum ) pl = rdscal*sqrt( dsum ) if( pl==zero ) then pl = one else pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) ) end if rdscal = zero dsum = one call stdlib${ii}$_${ri}$lassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum ) pr = rdscal*sqrt( dsum ) if( pr==zero ) then pr = one else pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) ) end if end if if( wantd ) then ! compute estimates of difu and difl. if( wantd1 ) then n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = idifjb ! frobenius norm-based difu-estimate. call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(& i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) ! frobenius norm-based difl-estimate. call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),& ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( 2_${ik}$*n1*n2+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) else ! compute 1-norm-based estimates of difu and difl using ! reversed communication with stdlib${ii}$_${ri}$lacn2. in each step a ! generalized sylvester equation or a transposed variant ! is solved. kase = 0_${ik}$ n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ mn2 = 2_${ik}$*n1*n2 ! 1-norm-based estimate of difu. 40 continue call stdlib${ii}$_${ri}$lacn2( mn2, work( mn2+1 ), work, iwork, dif( 1_${ik}$ ),kase, isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation. call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_${ri}$tgsyl( 'T', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 40 end if dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ ) ! 1-norm-based estimate of difl. 50 continue call stdlib${ii}$_${ri}$lacn2( mn2, work( mn2+1 ), work, iwork, dif( 2_${ik}$ ),kase, isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation. call stdlib${ii}$_${ri}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_${ri}$tgsyl( 'T', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( 2_${ik}$*n1*n2+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 50 end if dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ ) end if end if 60 continue ! compute generalized eigenvalues of reordered pair (a, b) and ! normalize the generalized schur form. pair = .false. loop_80: do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )/=zero ) then pair = .true. end if end if if( pair ) then ! compute the eigenvalue(s) at position k. work( 1_${ik}$ ) = a( k, k ) work( 2_${ik}$ ) = a( k+1, k ) work( 3_${ik}$ ) = a( k, k+1 ) work( 4_${ik}$ ) = a( k+1, k+1 ) work( 5_${ik}$ ) = b( k, k ) work( 6_${ik}$ ) = b( k+1, k ) work( 7_${ik}$ ) = b( k, k+1 ) work( 8_${ik}$ ) = b( k+1, k+1 ) call stdlib${ii}$_${ri}$lag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta( k ),beta( k+1 ), & alphar( k ), alphar( k+1 ),alphai( k ) ) alphai( k+1 ) = -alphai( k ) else if( sign( one, b( k, k ) )<zero ) then ! if b(k,k) is negative, make it positive do i = 1, n a( k, i ) = -a( k, i ) b( k, i ) = -b( k, i ) if( wantq ) q( i, k ) = -q( i, k ) end do end if alphar( k ) = a( k, k ) alphai( k ) = zero beta( k ) = b( k, k ) end if end if end do loop_80 work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_${ri}$tgsen #:endif #:endfor pure module subroutine stdlib${ii}$_ctgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alpha, beta, q, & !! CTGSEN reorders the generalized Schur decomposition of a complex !! matrix pair (A, B) (in terms of an unitary equivalence trans- !! formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues !! appears in the leading diagonal blocks of the pair (A,B). The leading !! columns of Q and Z form unitary bases of the corresponding left and !! right eigenspaces (deflating subspaces). (A, B) must be in !! generalized Schur canonical form, that is, A and B are both upper !! triangular. !! CTGSEN also computes the generalized eigenvalues !! w(j)= ALPHA(j) / BETA(j) !! of the reordered matrix pair (A, B). !! Optionally, the routine computes estimates of reciprocal condition !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to !! the selected cluster and the eigenvalues outside the cluster, resp., !! and norms of "projections" onto left and right eigenspaces w.r.t. !! the selected cluster in the (1,1)-block. ldq, z, ldz, m, pl, pr, dif,work, lwork, iwork, liwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: wantq, wantz integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(sp), intent(out) :: pl, pr ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(out) :: dif(*) complex(sp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*) complex(sp), intent(out) :: alpha(*), beta(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: idifjb = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, swap, wantd, wantd1, wantd2, wantp integer(${ik}$) :: i, ierr, ijb, k, kase, ks, liwmin, lwmin, mn2, n1, n2 real(sp) :: dscale, dsum, rdscal, safmin complex(sp) :: temp1, temp2 ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -9_${ik}$ else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then info = -13_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -15_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CTGSEN', -info ) return end if ierr = 0_${ik}$ wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$ wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$ wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$ wantd = wantd1 .or. wantd2 ! set m to the dimension of the specified pair of deflating ! subspaces. m = 0_${ik}$ if( .not.lquery .or. ijob/=0_${ik}$ ) then do k = 1, n alpha( k ) = a( k, k ) beta( k ) = b( k, k ) if( k<n ) then if( select( k ) )m = m + 1_${ik}$ else if( select( n ) )m = m + 1_${ik}$ end if end do end if if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then lwmin = max( 1_${ik}$, 2_${ik}$*m*(n-m) ) liwmin = max( 1_${ik}$, n+2 ) else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*m*(n-m) ) liwmin = max( 1_${ik}$, 2_${ik}$*m*(n-m), n+2 ) else lwmin = 1_${ik}$ liwmin = 1_${ik}$ end if work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -21_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -23_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CTGSEN', -info ) return else if( lquery ) then return end if ! quick return if possible. if( m==n .or. m==0_${ik}$ ) then if( wantp ) then pl = one pr = one end if if( wantd ) then dscale = zero dsum = one do i = 1, n call stdlib${ii}$_classq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) call stdlib${ii}$_classq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) end do dif( 1_${ik}$ ) = dscale*sqrt( dsum ) dif( 2_${ik}$ ) = dif( 1_${ik}$ ) end if go to 70 end if ! get machine constant safmin = stdlib${ii}$_slamch( 'S' ) ! collect the selected blocks at the top-left corner of (a, b). ks = 0_${ik}$ do k = 1, n swap = select( k ) if( swap ) then ks = ks + 1_${ik}$ ! swap the k-th block to position ks. compute unitary q ! and z that will swap adjacent diagonal blocks in (a, b). if( k/=ks )call stdlib${ii}$_ctgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,ldz, k,& ks, ierr ) if( ierr>0_${ik}$ ) then ! swap is rejected: exit. info = 1_${ik}$ if( wantp ) then pl = zero pr = zero end if if( wantd ) then dif( 1_${ik}$ ) = zero dif( 2_${ik}$ ) = zero end if go to 70 end if end if end do if( wantp ) then ! solve generalized sylvester equation for r and l: ! a11 * r - l * a22 = a12 ! b11 * r - l * b22 = b12 n1 = m n2 = n - m i = n1 + 1_${ik}$ call stdlib${ii}$_clacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 ) call stdlib${ii}$_clacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 ) ijb = 0_${ik}$ call stdlib${ii}$_ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,& i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, & iwork, ierr ) ! estimate the reciprocal of norms of "projections" onto ! left and right eigenspaces rdscal = zero dsum = one call stdlib${ii}$_classq( n1*n2, work, 1_${ik}$, rdscal, dsum ) pl = rdscal*sqrt( dsum ) if( pl==zero ) then pl = one else pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) ) end if rdscal = zero dsum = one call stdlib${ii}$_classq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum ) pr = rdscal*sqrt( dsum ) if( pr==zero ) then pr = one else pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) ) end if end if if( wantd ) then ! compute estimates difu and difl. if( wantd1 ) then n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = idifjb ! frobenius norm-based difu estimate. call stdlib${ii}$_ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(& i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) ! frobenius norm-based difl estimate. call stdlib${ii}$_ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),& ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) else ! compute 1-norm-based estimates of difu and difl using ! reversed communication with stdlib${ii}$_clacn2. in each step a ! generalized sylvester equation or a transposed variant ! is solved. kase = 0_${ik}$ n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ mn2 = 2_${ik}$*n1*n2 ! 1-norm-based estimate of difu. 40 continue call stdlib${ii}$_clacn2( mn2, work( mn2+1 ), work, dif( 1_${ik}$ ), kase,isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation call stdlib${ii}$_ctgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_ctgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 40 end if dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ ) ! 1-norm-based estimate of difl. 50 continue call stdlib${ii}$_clacn2( mn2, work( mn2+1 ), work, dif( 2_${ik}$ ), kase,isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation call stdlib${ii}$_ctgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_ctgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 50 end if dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ ) end if end if ! if b(k,k) is complex, make it real and positive (normalization ! of the generalized schur form) and store the generalized ! eigenvalues of reordered pair (a, b) do k = 1, n dscale = abs( b( k, k ) ) if( dscale>safmin ) then temp1 = conjg( b( k, k ) / dscale ) temp2 = b( k, k ) / dscale b( k, k ) = dscale call stdlib${ii}$_cscal( n-k, temp1, b( k, k+1 ), ldb ) call stdlib${ii}$_cscal( n-k+1, temp1, a( k, k ), lda ) if( wantq )call stdlib${ii}$_cscal( n, temp2, q( 1_${ik}$, k ), 1_${ik}$ ) else b( k, k ) = cmplx( zero, zero,KIND=sp) end if alpha( k ) = a( k, k ) beta( k ) = b( k, k ) end do 70 continue work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_ctgsen pure module subroutine stdlib${ii}$_ztgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alpha, beta, q, & !! ZTGSEN reorders the generalized Schur decomposition of a complex !! matrix pair (A, B) (in terms of an unitary equivalence trans- !! formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues !! appears in the leading diagonal blocks of the pair (A,B). The leading !! columns of Q and Z form unitary bases of the corresponding left and !! right eigenspaces (deflating subspaces). (A, B) must be in !! generalized Schur canonical form, that is, A and B are both upper !! triangular. !! ZTGSEN also computes the generalized eigenvalues !! w(j)= ALPHA(j) / BETA(j) !! of the reordered matrix pair (A, B). !! Optionally, the routine computes estimates of reciprocal condition !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to !! the selected cluster and the eigenvalues outside the cluster, resp., !! and norms of "projections" onto left and right eigenspaces w.r.t. !! the selected cluster in the (1,1)-block. ldq, z, ldz, m, pl, pr, dif,work, lwork, iwork, liwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: wantq, wantz integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(dp), intent(out) :: pl, pr ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(out) :: dif(*) complex(dp), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*) complex(dp), intent(out) :: alpha(*), beta(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: idifjb = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, swap, wantd, wantd1, wantd2, wantp integer(${ik}$) :: i, ierr, ijb, k, kase, ks, liwmin, lwmin, mn2, n1, n2 real(dp) :: dscale, dsum, rdscal, safmin complex(dp) :: temp1, temp2 ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -9_${ik}$ else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then info = -13_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -15_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZTGSEN', -info ) return end if ierr = 0_${ik}$ wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$ wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$ wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$ wantd = wantd1 .or. wantd2 ! set m to the dimension of the specified pair of deflating ! subspaces. m = 0_${ik}$ if( .not.lquery .or. ijob/=0_${ik}$ ) then do k = 1, n alpha( k ) = a( k, k ) beta( k ) = b( k, k ) if( k<n ) then if( select( k ) )m = m + 1_${ik}$ else if( select( n ) )m = m + 1_${ik}$ end if end do end if if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then lwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, n+2 ) else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+2 ) else lwmin = 1_${ik}$ liwmin = 1_${ik}$ end if work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -21_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -23_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZTGSEN', -info ) return else if( lquery ) then return end if ! quick return if possible. if( m==n .or. m==0_${ik}$ ) then if( wantp ) then pl = one pr = one end if if( wantd ) then dscale = zero dsum = one do i = 1, n call stdlib${ii}$_zlassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) call stdlib${ii}$_zlassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) end do dif( 1_${ik}$ ) = dscale*sqrt( dsum ) dif( 2_${ik}$ ) = dif( 1_${ik}$ ) end if go to 70 end if ! get machine constant safmin = stdlib${ii}$_dlamch( 'S' ) ! collect the selected blocks at the top-left corner of (a, b). ks = 0_${ik}$ do k = 1, n swap = select( k ) if( swap ) then ks = ks + 1_${ik}$ ! swap the k-th block to position ks. compute unitary q ! and z that will swap adjacent diagonal blocks in (a, b). if( k/=ks )call stdlib${ii}$_ztgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,ldz, k,& ks, ierr ) if( ierr>0_${ik}$ ) then ! swap is rejected: exit. info = 1_${ik}$ if( wantp ) then pl = zero pr = zero end if if( wantd ) then dif( 1_${ik}$ ) = zero dif( 2_${ik}$ ) = zero end if go to 70 end if end if end do if( wantp ) then ! solve generalized sylvester equation for r and l: ! a11 * r - l * a22 = a12 ! b11 * r - l * b22 = b12 n1 = m n2 = n - m i = n1 + 1_${ik}$ call stdlib${ii}$_zlacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 ) call stdlib${ii}$_zlacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 ) ijb = 0_${ik}$ call stdlib${ii}$_ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,& i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, & iwork, ierr ) ! estimate the reciprocal of norms of "projections" onto ! left and right eigenspaces rdscal = zero dsum = one call stdlib${ii}$_zlassq( n1*n2, work, 1_${ik}$, rdscal, dsum ) pl = rdscal*sqrt( dsum ) if( pl==zero ) then pl = one else pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) ) end if rdscal = zero dsum = one call stdlib${ii}$_zlassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum ) pr = rdscal*sqrt( dsum ) if( pr==zero ) then pr = one else pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) ) end if end if if( wantd ) then ! compute estimates difu and difl. if( wantd1 ) then n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = idifjb ! frobenius norm-based difu estimate. call stdlib${ii}$_ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(& i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) ! frobenius norm-based difl estimate. call stdlib${ii}$_ztgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),& ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) else ! compute 1-norm-based estimates of difu and difl using ! reversed communication with stdlib${ii}$_zlacn2. in each step a ! generalized sylvester equation or a transposed variant ! is solved. kase = 0_${ik}$ n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ mn2 = 2_${ik}$*n1*n2 ! 1-norm-based estimate of difu. 40 continue call stdlib${ii}$_zlacn2( mn2, work( mn2+1 ), work, dif( 1_${ik}$ ), kase,isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation call stdlib${ii}$_ztgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_ztgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 40 end if dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ ) ! 1-norm-based estimate of difl. 50 continue call stdlib${ii}$_zlacn2( mn2, work( mn2+1 ), work, dif( 2_${ik}$ ), kase,isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation call stdlib${ii}$_ztgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_ztgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 50 end if dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ ) end if end if ! if b(k,k) is complex, make it real and positive (normalization ! of the generalized schur form) and store the generalized ! eigenvalues of reordered pair (a, b) do k = 1, n dscale = abs( b( k, k ) ) if( dscale>safmin ) then temp1 = conjg( b( k, k ) / dscale ) temp2 = b( k, k ) / dscale b( k, k ) = dscale call stdlib${ii}$_zscal( n-k, temp1, b( k, k+1 ), ldb ) call stdlib${ii}$_zscal( n-k+1, temp1, a( k, k ), lda ) if( wantq )call stdlib${ii}$_zscal( n, temp2, q( 1_${ik}$, k ), 1_${ik}$ ) else b( k, k ) = cmplx( zero, zero,KIND=dp) end if alpha( k ) = a( k, k ) beta( k ) = b( k, k ) end do 70 continue work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_ztgsen #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$tgsen( ijob, wantq, wantz, select, n, a, lda, b, ldb,alpha, beta, q, & !! ZTGSEN: reorders the generalized Schur decomposition of a complex !! matrix pair (A, B) (in terms of an unitary equivalence trans- !! formation Q**H * (A, B) * Z), so that a selected cluster of eigenvalues !! appears in the leading diagonal blocks of the pair (A,B). The leading !! columns of Q and Z form unitary bases of the corresponding left and !! right eigenspaces (deflating subspaces). (A, B) must be in !! generalized Schur canonical form, that is, A and B are both upper !! triangular. !! ZTGSEN also computes the generalized eigenvalues !! w(j)= ALPHA(j) / BETA(j) !! of the reordered matrix pair (A, B). !! Optionally, the routine computes estimates of reciprocal condition !! numbers for eigenvalues and eigenspaces. These are Difu[(A11,B11), !! (A22,B22)] and Difl[(A11,B11), (A22,B22)], i.e. the separation(s) !! between the matrix pairs (A11, B11) and (A22,B22) that correspond to !! the selected cluster and the eigenvalues outside the cluster, resp., !! and norms of "projections" onto left and right eigenspaces w.r.t. !! the selected cluster in the (1,1)-block. ldq, z, ldz, m, pl, pr, dif,work, lwork, iwork, liwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments logical(lk), intent(in) :: wantq, wantz integer(${ik}$), intent(in) :: ijob, lda, ldb, ldq, ldz, liwork, lwork, n integer(${ik}$), intent(out) :: info, m real(${ck}$), intent(out) :: pl, pr ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(${ck}$), intent(out) :: dif(*) complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*), q(ldq,*), z(ldz,*) complex(${ck}$), intent(out) :: alpha(*), beta(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: idifjb = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, swap, wantd, wantd1, wantd2, wantp integer(${ik}$) :: i, ierr, ijb, k, kase, ks, liwmin, lwmin, mn2, n1, n2 real(${ck}$) :: dscale, dsum, rdscal, safmin complex(${ck}$) :: temp1, temp2 ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ .or. liwork==-1_${ik}$ ) if( ijob<0_${ik}$ .or. ijob>5_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -5_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -9_${ik}$ else if( ldq<1_${ik}$ .or. ( wantq .and. ldq<n ) ) then info = -13_${ik}$ else if( ldz<1_${ik}$ .or. ( wantz .and. ldz<n ) ) then info = -15_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZTGSEN', -info ) return end if ierr = 0_${ik}$ wantp = ijob==1_${ik}$ .or. ijob>=4_${ik}$ wantd1 = ijob==2_${ik}$ .or. ijob==4_${ik}$ wantd2 = ijob==3_${ik}$ .or. ijob==5_${ik}$ wantd = wantd1 .or. wantd2 ! set m to the dimension of the specified pair of deflating ! subspaces. m = 0_${ik}$ if( .not.lquery .or. ijob/=0_${ik}$ ) then do k = 1, n alpha( k ) = a( k, k ) beta( k ) = b( k, k ) if( k<n ) then if( select( k ) )m = m + 1_${ik}$ else if( select( n ) )m = m + 1_${ik}$ end if end do end if if( ijob==1_${ik}$ .or. ijob==2_${ik}$ .or. ijob==4_${ik}$ ) then lwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, n+2 ) else if( ijob==3_${ik}$ .or. ijob==5_${ik}$ ) then lwmin = max( 1_${ik}$, 4_${ik}$*m*( n-m ) ) liwmin = max( 1_${ik}$, 2_${ik}$*m*( n-m ), n+2 ) else lwmin = 1_${ik}$ liwmin = 1_${ik}$ end if work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin if( lwork<lwmin .and. .not.lquery ) then info = -21_${ik}$ else if( liwork<liwmin .and. .not.lquery ) then info = -23_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZTGSEN', -info ) return else if( lquery ) then return end if ! quick return if possible. if( m==n .or. m==0_${ik}$ ) then if( wantp ) then pl = one pr = one end if if( wantd ) then dscale = zero dsum = one do i = 1, n call stdlib${ii}$_${ci}$lassq( n, a( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) call stdlib${ii}$_${ci}$lassq( n, b( 1_${ik}$, i ), 1_${ik}$, dscale, dsum ) end do dif( 1_${ik}$ ) = dscale*sqrt( dsum ) dif( 2_${ik}$ ) = dif( 1_${ik}$ ) end if go to 70 end if ! get machine constant safmin = stdlib${ii}$_${c2ri(ci)}$lamch( 'S' ) ! collect the selected blocks at the top-left corner of (a, b). ks = 0_${ik}$ do k = 1, n swap = select( k ) if( swap ) then ks = ks + 1_${ik}$ ! swap the k-th block to position ks. compute unitary q ! and z that will swap adjacent diagonal blocks in (a, b). if( k/=ks )call stdlib${ii}$_${ci}$tgexc( wantq, wantz, n, a, lda, b, ldb, q, ldq, z,ldz, k,& ks, ierr ) if( ierr>0_${ik}$ ) then ! swap is rejected: exit. info = 1_${ik}$ if( wantp ) then pl = zero pr = zero end if if( wantd ) then dif( 1_${ik}$ ) = zero dif( 2_${ik}$ ) = zero end if go to 70 end if end if end do if( wantp ) then ! solve generalized sylvester equation for r and l: ! a11 * r - l * a22 = a12 ! b11 * r - l * b22 = b12 n1 = m n2 = n - m i = n1 + 1_${ik}$ call stdlib${ii}$_${ci}$lacpy( 'FULL', n1, n2, a( 1_${ik}$, i ), lda, work, n1 ) call stdlib${ii}$_${ci}$lacpy( 'FULL', n1, n2, b( 1_${ik}$, i ), ldb, work( n1*n2+1 ),n1 ) ijb = 0_${ik}$ call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b( i,& i ), ldb, work( n1*n2+1 ), n1,dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-2*n1*n2, & iwork, ierr ) ! estimate the reciprocal of norms of "projections" onto ! left and right eigenspaces rdscal = zero dsum = one call stdlib${ii}$_${ci}$lassq( n1*n2, work, 1_${ik}$, rdscal, dsum ) pl = rdscal*sqrt( dsum ) if( pl==zero ) then pl = one else pl = dscale / ( sqrt( dscale*dscale / pl+pl )*sqrt( pl ) ) end if rdscal = zero dsum = one call stdlib${ii}$_${ci}$lassq( n1*n2, work( n1*n2+1 ), 1_${ik}$, rdscal, dsum ) pr = rdscal*sqrt( dsum ) if( pr==zero ) then pr = one else pr = dscale / ( sqrt( dscale*dscale / pr+pr )*sqrt( pr ) ) end if end if if( wantd ) then ! compute estimates difu and difl. if( wantd1 ) then n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = idifjb ! frobenius norm-based difu estimate. call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda, work,n1, b, ldb, b(& i, i ), ldb, work( n1*n2+1 ),n1, dscale, dif( 1_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) ! frobenius norm-based difl estimate. call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda, work,n2, b( i, i ),& ldb, b, ldb, work( n1*n2+1 ),n2, dscale, dif( 2_${ik}$ ), work( n1*n2*2_${ik}$+1 ),lwork-& 2_${ik}$*n1*n2, iwork, ierr ) else ! compute 1-norm-based estimates of difu and difl using ! reversed communication with stdlib${ii}$_${ci}$lacn2. in each step a ! generalized sylvester equation or a transposed variant ! is solved. kase = 0_${ik}$ n1 = m n2 = n - m i = n1 + 1_${ik}$ ijb = 0_${ik}$ mn2 = 2_${ik}$*n1*n2 ! 1-norm-based estimate of difu. 40 continue call stdlib${ii}$_${ci}$lacn2( mn2, work( mn2+1 ), work, dif( 1_${ik}$ ), kase,isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_${ci}$tgsyl( 'C', ijb, n1, n2, a, lda, a( i, i ), lda,work, n1, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n1, dscale, dif( 1_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 40 end if dif( 1_${ik}$ ) = dscale / dif( 1_${ik}$ ) ! 1-norm-based estimate of difl. 50 continue call stdlib${ii}$_${ci}$lacn2( mn2, work( mn2+1 ), work, dif( 2_${ik}$ ), kase,isave ) if( kase/=0_${ik}$ ) then if( kase==1_${ik}$ ) then ! solve generalized sylvester equation call stdlib${ii}$_${ci}$tgsyl( 'N', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b( & i, i ), ldb, b, ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) else ! solve the transposed variant. call stdlib${ii}$_${ci}$tgsyl( 'C', ijb, n2, n1, a( i, i ), lda, a, lda,work, n2, b, & ldb, b( i, i ), ldb,work( n1*n2+1 ), n2, dscale, dif( 2_${ik}$ ),work( n1*n2*2_${ik}$+1 )& , lwork-2*n1*n2, iwork,ierr ) end if go to 50 end if dif( 2_${ik}$ ) = dscale / dif( 2_${ik}$ ) end if end if ! if b(k,k) is complex, make it real and positive (normalization ! of the generalized schur form) and store the generalized ! eigenvalues of reordered pair (a, b) do k = 1, n dscale = abs( b( k, k ) ) if( dscale>safmin ) then temp1 = conjg( b( k, k ) / dscale ) temp2 = b( k, k ) / dscale b( k, k ) = dscale call stdlib${ii}$_${ci}$scal( n-k, temp1, b( k, k+1 ), ldb ) call stdlib${ii}$_${ci}$scal( n-k+1, temp1, a( k, k ), lda ) if( wantq )call stdlib${ii}$_${ci}$scal( n, temp2, q( 1_${ik}$, k ), 1_${ik}$ ) else b( k, k ) = cmplx( zero, zero,KIND=${ck}$) end if alpha( k ) = a( k, k ) beta( k ) = b( k, k ) end do 70 continue work( 1_${ik}$ ) = lwmin iwork( 1_${ik}$ ) = liwmin return end subroutine stdlib${ii}$_${ci}$tgsen #:endif #:endfor pure module subroutine stdlib${ii}$_stgsna( job, howmny, select, n, a, lda, b, ldb, vl,ldvl, vr, ldvr, s, & !! STGSNA estimates reciprocal condition numbers for specified !! eigenvalues and/or eigenvectors of a matrix pair (A, B) in !! generalized real Schur canonical form (or of any matrix pair !! (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where !! Z**T denotes the transpose of Z. !! (A, B) must be in generalized real Schur form (as returned by SGGES), !! i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal !! blocks. B is upper triangular. dif, mm, m, work, lwork,iwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: howmny, job integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: lda, ldb, ldvl, ldvr, lwork, mm, n ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(in) :: a(lda,*), b(ldb,*), vl(ldvl,*), vr(ldvr,*) real(sp), intent(out) :: dif(*), s(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: difdri = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, pair, somcon, wantbh, wantdf, wants integer(${ik}$) :: i, ierr, ifst, ilst, iz, k, ks, lwmin, n1, n2 real(sp) :: alphai, alphar, alprqt, beta, c1, c2, cond, eps, lnrm, rnrm, root1, root2, & scale, smlnum, tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv, uhbvi ! Local Arrays real(sp) :: dummy(1_${ik}$), dummy1(1_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters wantbh = stdlib_lsame( job, 'B' ) wants = stdlib_lsame( job, 'E' ) .or. wantbh wantdf = stdlib_lsame( job, 'V' ) .or. wantbh somcon = stdlib_lsame( howmny, 'S' ) info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ ) if( .not.wants .and. .not.wantdf ) then info = -1_${ik}$ else if( .not.stdlib_lsame( howmny, 'A' ) .and. .not.somcon ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -4_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( wants .and. ldvl<n ) then info = -10_${ik}$ else if( wants .and. ldvr<n ) then info = -12_${ik}$ else ! set m to the number of eigenpairs for which condition numbers ! are required, and test mm. if( somcon ) then m = 0_${ik}$ pair = .false. do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$ end if else if( select( n ) )m = m + 1_${ik}$ end if end if end do else m = n end if if( n==0_${ik}$ ) then lwmin = 1_${ik}$ else if( stdlib_lsame( job, 'V' ) .or. stdlib_lsame( job, 'B' ) ) then lwmin = 2_${ik}$*n*( n + 2_${ik}$ ) + 16_${ik}$ else lwmin = n end if work( 1_${ik}$ ) = lwmin if( mm<m ) then info = -15_${ik}$ else if( lwork<lwmin .and. .not.lquery ) then info = -18_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'STGSNA', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return ! get machine constants eps = stdlib${ii}$_slamch( 'P' ) smlnum = stdlib${ii}$_slamch( 'S' ) / eps ks = 0_${ik}$ pair = .false. loop_20: do k = 1, n ! determine whether a(k,k) begins a 1-by-1 or 2-by-2 block. if( pair ) then pair = .false. cycle loop_20 else if( k<n )pair = a( k+1, k )/=zero end if ! determine whether condition numbers are required for the k-th ! eigenpair. if( somcon ) then if( pair ) then if( .not.select( k ) .and. .not.select( k+1 ) )cycle loop_20 else if( .not.select( k ) )cycle loop_20 end if end if ks = ks + 1_${ik}$ if( wants ) then ! compute the reciprocal condition number of the k-th ! eigenvalue. if( pair ) then ! complex eigenvalue pair. rnrm = stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, vr( & 1_${ik}$, ks+1 ), 1_${ik}$ ) ) lnrm = stdlib${ii}$_slapy2( stdlib${ii}$_snrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_snrm2( n, vl( & 1_${ik}$, ks+1 ), 1_${ik}$ ) ) call stdlib${ii}$_sgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) tmprr = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) tmpri = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) call stdlib${ii}$_sgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ ) tmpii = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) tmpir = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) uhav = tmprr + tmpii uhavi = tmpir - tmpri call stdlib${ii}$_sgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) tmprr = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) tmpri = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) call stdlib${ii}$_sgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ ) tmpii = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) tmpir = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) uhbv = tmprr + tmpii uhbvi = tmpir - tmpri uhav = stdlib${ii}$_slapy2( uhav, uhavi ) uhbv = stdlib${ii}$_slapy2( uhbv, uhbvi ) cond = stdlib${ii}$_slapy2( uhav, uhbv ) s( ks ) = cond / ( rnrm*lnrm ) s( ks+1 ) = s( ks ) else ! real eigenvalue. rnrm = stdlib${ii}$_snrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ) lnrm = stdlib${ii}$_snrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ) call stdlib${ii}$_sgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) uhav = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) call stdlib${ii}$_sgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) uhbv = stdlib${ii}$_sdot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) cond = stdlib${ii}$_slapy2( uhav, uhbv ) if( cond==zero ) then s( ks ) = -one else s( ks ) = cond / ( rnrm*lnrm ) end if end if end if if( wantdf ) then if( n==1_${ik}$ ) then dif( ks ) = stdlib${ii}$_slapy2( a( 1_${ik}$, 1_${ik}$ ), b( 1_${ik}$, 1_${ik}$ ) ) cycle loop_20 end if ! estimate the reciprocal condition number of the k-th ! eigenvectors. if( pair ) then ! copy the 2-by 2 pencil beginning at (a(k,k), b(k, k)). ! compute the eigenvalue(s) at position k. work( 1_${ik}$ ) = a( k, k ) work( 2_${ik}$ ) = a( k+1, k ) work( 3_${ik}$ ) = a( k, k+1 ) work( 4_${ik}$ ) = a( k+1, k+1 ) work( 5_${ik}$ ) = b( k, k ) work( 6_${ik}$ ) = b( k+1, k ) work( 7_${ik}$ ) = b( k, k+1 ) work( 8_${ik}$ ) = b( k+1, k+1 ) call stdlib${ii}$_slag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta,dummy1( 1_${ik}$ ), & alphar, dummy( 1_${ik}$ ), alphai ) alprqt = one c1 = two*( alphar*alphar+alphai*alphai+beta*beta ) c2 = four*beta*beta*alphai*alphai root1 = c1 + sqrt( c1*c1-4.0_sp*c2 ) root2 = c2 / root1 root1 = root1 / two cond = min( sqrt( root1 ), sqrt( root2 ) ) end if ! copy the matrix (a, b) to the array work and swap the ! diagonal block beginning at a(k,k) to the (1,1) position. call stdlib${ii}$_slacpy( 'FULL', n, n, a, lda, work, n ) call stdlib${ii}$_slacpy( 'FULL', n, n, b, ldb, work( n*n+1 ), n ) ifst = k ilst = 1_${ik}$ call stdlib${ii}$_stgexc( .false., .false., n, work, n, work( n*n+1 ), n,dummy, 1_${ik}$, & dummy1, 1_${ik}$, ifst, ilst,work( n*n*2_${ik}$+1 ), lwork-2*n*n, ierr ) if( ierr>0_${ik}$ ) then ! ill-conditioned problem - swap rejected. dif( ks ) = zero else ! reordering successful, solve generalized sylvester ! equation for r and l, ! a22 * r - l * a11 = a12 ! b22 * r - l * b11 = b12, ! and compute estimate of difl((a11,b11), (a22, b22)). n1 = 1_${ik}$ if( work( 2_${ik}$ )/=zero )n1 = 2_${ik}$ n2 = n - n1 if( n2==0_${ik}$ ) then dif( ks ) = cond else i = n*n + 1_${ik}$ iz = 2_${ik}$*n*n + 1_${ik}$ call stdlib${ii}$_stgsyl( 'N', difdri, n2, n1, work( n*n1+n1+1 ),n, work, n, & work( n1+1 ), n,work( n*n1+n1+i ), n, work( i ), n,work( n1+i ), n, scale, & dif( ks ),work( iz+1 ), lwork-2*n*n, iwork, ierr ) if( pair )dif( ks ) = min( max( one, alprqt )*dif( ks ),cond ) end if end if if( pair )dif( ks+1 ) = dif( ks ) end if if( pair )ks = ks + 1_${ik}$ end do loop_20 work( 1_${ik}$ ) = lwmin return end subroutine stdlib${ii}$_stgsna pure module subroutine stdlib${ii}$_dtgsna( job, howmny, select, n, a, lda, b, ldb, vl,ldvl, vr, ldvr, s, & !! DTGSNA estimates reciprocal condition numbers for specified !! eigenvalues and/or eigenvectors of a matrix pair (A, B) in !! generalized real Schur canonical form (or of any matrix pair !! (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where !! Z**T denotes the transpose of Z. !! (A, B) must be in generalized real Schur form (as returned by DGGES), !! i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal !! blocks. B is upper triangular. dif, mm, m, work, lwork,iwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: howmny, job integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: lda, ldb, ldvl, ldvr, lwork, mm, n ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(in) :: a(lda,*), b(ldb,*), vl(ldvl,*), vr(ldvr,*) real(dp), intent(out) :: dif(*), s(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: difdri = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, pair, somcon, wantbh, wantdf, wants integer(${ik}$) :: i, ierr, ifst, ilst, iz, k, ks, lwmin, n1, n2 real(dp) :: alphai, alphar, alprqt, beta, c1, c2, cond, eps, lnrm, rnrm, root1, root2, & scale, smlnum, tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv, uhbvi ! Local Arrays real(dp) :: dummy(1_${ik}$), dummy1(1_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters wantbh = stdlib_lsame( job, 'B' ) wants = stdlib_lsame( job, 'E' ) .or. wantbh wantdf = stdlib_lsame( job, 'V' ) .or. wantbh somcon = stdlib_lsame( howmny, 'S' ) info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ ) if( .not.wants .and. .not.wantdf ) then info = -1_${ik}$ else if( .not.stdlib_lsame( howmny, 'A' ) .and. .not.somcon ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -4_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( wants .and. ldvl<n ) then info = -10_${ik}$ else if( wants .and. ldvr<n ) then info = -12_${ik}$ else ! set m to the number of eigenpairs for which condition numbers ! are required, and test mm. if( somcon ) then m = 0_${ik}$ pair = .false. do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$ end if else if( select( n ) )m = m + 1_${ik}$ end if end if end do else m = n end if if( n==0_${ik}$ ) then lwmin = 1_${ik}$ else if( stdlib_lsame( job, 'V' ) .or. stdlib_lsame( job, 'B' ) ) then lwmin = 2_${ik}$*n*( n + 2_${ik}$ ) + 16_${ik}$ else lwmin = n end if work( 1_${ik}$ ) = lwmin if( mm<m ) then info = -15_${ik}$ else if( lwork<lwmin .and. .not.lquery ) then info = -18_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTGSNA', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return ! get machine constants eps = stdlib${ii}$_dlamch( 'P' ) smlnum = stdlib${ii}$_dlamch( 'S' ) / eps ks = 0_${ik}$ pair = .false. loop_20: do k = 1, n ! determine whether a(k,k) begins a 1-by-1 or 2-by-2 block. if( pair ) then pair = .false. cycle loop_20 else if( k<n )pair = a( k+1, k )/=zero end if ! determine whether condition numbers are required for the k-th ! eigenpair. if( somcon ) then if( pair ) then if( .not.select( k ) .and. .not.select( k+1 ) )cycle loop_20 else if( .not.select( k ) )cycle loop_20 end if end if ks = ks + 1_${ik}$ if( wants ) then ! compute the reciprocal condition number of the k-th ! eigenvalue. if( pair ) then ! complex eigenvalue pair. rnrm = stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, vr( & 1_${ik}$, ks+1 ), 1_${ik}$ ) ) lnrm = stdlib${ii}$_dlapy2( stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_dnrm2( n, vl( & 1_${ik}$, ks+1 ), 1_${ik}$ ) ) call stdlib${ii}$_dgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) tmprr = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) tmpri = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) call stdlib${ii}$_dgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ ) tmpii = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) tmpir = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) uhav = tmprr + tmpii uhavi = tmpir - tmpri call stdlib${ii}$_dgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) tmprr = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) tmpri = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) call stdlib${ii}$_dgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ ) tmpii = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) tmpir = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) uhbv = tmprr + tmpii uhbvi = tmpir - tmpri uhav = stdlib${ii}$_dlapy2( uhav, uhavi ) uhbv = stdlib${ii}$_dlapy2( uhbv, uhbvi ) cond = stdlib${ii}$_dlapy2( uhav, uhbv ) s( ks ) = cond / ( rnrm*lnrm ) s( ks+1 ) = s( ks ) else ! real eigenvalue. rnrm = stdlib${ii}$_dnrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ) lnrm = stdlib${ii}$_dnrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ) call stdlib${ii}$_dgemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) uhav = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) call stdlib${ii}$_dgemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) uhbv = stdlib${ii}$_ddot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) cond = stdlib${ii}$_dlapy2( uhav, uhbv ) if( cond==zero ) then s( ks ) = -one else s( ks ) = cond / ( rnrm*lnrm ) end if end if end if if( wantdf ) then if( n==1_${ik}$ ) then dif( ks ) = stdlib${ii}$_dlapy2( a( 1_${ik}$, 1_${ik}$ ), b( 1_${ik}$, 1_${ik}$ ) ) cycle loop_20 end if ! estimate the reciprocal condition number of the k-th ! eigenvectors. if( pair ) then ! copy the 2-by 2 pencil beginning at (a(k,k), b(k, k)). ! compute the eigenvalue(s) at position k. work( 1_${ik}$ ) = a( k, k ) work( 2_${ik}$ ) = a( k+1, k ) work( 3_${ik}$ ) = a( k, k+1 ) work( 4_${ik}$ ) = a( k+1, k+1 ) work( 5_${ik}$ ) = b( k, k ) work( 6_${ik}$ ) = b( k+1, k ) work( 7_${ik}$ ) = b( k, k+1 ) work( 8_${ik}$ ) = b( k+1, k+1 ) call stdlib${ii}$_dlag2( work, 2_${ik}$, work( 5_${ik}$ ), 2_${ik}$, smlnum*eps, beta,dummy1( 1_${ik}$ ), & alphar, dummy( 1_${ik}$ ), alphai ) alprqt = one c1 = two*( alphar*alphar+alphai*alphai+beta*beta ) c2 = four*beta*beta*alphai*alphai root1 = c1 + sqrt( c1*c1-4.0_dp*c2 ) root2 = c2 / root1 root1 = root1 / two cond = min( sqrt( root1 ), sqrt( root2 ) ) end if ! copy the matrix (a, b) to the array work and swap the ! diagonal block beginning at a(k,k) to the (1,1) position. call stdlib${ii}$_dlacpy( 'FULL', n, n, a, lda, work, n ) call stdlib${ii}$_dlacpy( 'FULL', n, n, b, ldb, work( n*n+1 ), n ) ifst = k ilst = 1_${ik}$ call stdlib${ii}$_dtgexc( .false., .false., n, work, n, work( n*n+1 ), n,dummy, 1_${ik}$, & dummy1, 1_${ik}$, ifst, ilst,work( n*n*2_${ik}$+1 ), lwork-2*n*n, ierr ) if( ierr>0_${ik}$ ) then ! ill-conditioned problem - swap rejected. dif( ks ) = zero else ! reordering successful, solve generalized sylvester ! equation for r and l, ! a22 * r - l * a11 = a12 ! b22 * r - l * b11 = b12, ! and compute estimate of difl((a11,b11), (a22, b22)). n1 = 1_${ik}$ if( work( 2_${ik}$ )/=zero )n1 = 2_${ik}$ n2 = n - n1 if( n2==0_${ik}$ ) then dif( ks ) = cond else i = n*n + 1_${ik}$ iz = 2_${ik}$*n*n + 1_${ik}$ call stdlib${ii}$_dtgsyl( 'N', difdri, n2, n1, work( n*n1+n1+1 ),n, work, n, & work( n1+1 ), n,work( n*n1+n1+i ), n, work( i ), n,work( n1+i ), n, scale, & dif( ks ),work( iz+1 ), lwork-2*n*n, iwork, ierr ) if( pair )dif( ks ) = min( max( one, alprqt )*dif( ks ),cond ) end if end if if( pair )dif( ks+1 ) = dif( ks ) end if if( pair )ks = ks + 1_${ik}$ end do loop_20 work( 1_${ik}$ ) = lwmin return end subroutine stdlib${ii}$_dtgsna #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$tgsna( job, howmny, select, n, a, lda, b, ldb, vl,ldvl, vr, ldvr, s, & !! DTGSNA: estimates reciprocal condition numbers for specified !! eigenvalues and/or eigenvectors of a matrix pair (A, B) in !! generalized real Schur canonical form (or of any matrix pair !! (Q*A*Z**T, Q*B*Z**T) with orthogonal matrices Q and Z, where !! Z**T denotes the transpose of Z. !! (A, B) must be in generalized real Schur form (as returned by DGGES), !! i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal !! blocks. B is upper triangular. dif, mm, m, work, lwork,iwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: howmny, job integer(${ik}$), intent(out) :: info, m integer(${ik}$), intent(in) :: lda, ldb, ldvl, ldvr, lwork, mm, n ! Array Arguments logical(lk), intent(in) :: select(*) integer(${ik}$), intent(out) :: iwork(*) real(${rk}$), intent(in) :: a(lda,*), b(ldb,*), vl(ldvl,*), vr(ldvr,*) real(${rk}$), intent(out) :: dif(*), s(*), work(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: difdri = 3_${ik}$ ! Local Scalars logical(lk) :: lquery, pair, somcon, wantbh, wantdf, wants integer(${ik}$) :: i, ierr, ifst, ilst, iz, k, ks, lwmin, n1, n2 real(${rk}$) :: alphai, alphar, alprqt, beta, c1, c2, cond, eps, lnrm, rnrm, root1, root2, & scale, smlnum, tmpii, tmpir, tmpri, tmprr, uhav, uhavi, uhbv, uhbvi ! Local Arrays real(${rk}$) :: dummy(1_${ik}$), dummy1(1_${ik}$) ! Intrinsic Functions ! Executable Statements ! decode and test the input parameters wantbh = stdlib_lsame( job, 'B' ) wants = stdlib_lsame( job, 'E' ) .or. wantbh wantdf = stdlib_lsame( job, 'V' ) .or. wantbh somcon = stdlib_lsame( howmny, 'S' ) info = 0_${ik}$ lquery = ( lwork==-1_${ik}$ ) if( .not.wants .and. .not.wantdf ) then info = -1_${ik}$ else if( .not.stdlib_lsame( howmny, 'A' ) .and. .not.somcon ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -4_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -8_${ik}$ else if( wants .and. ldvl<n ) then info = -10_${ik}$ else if( wants .and. ldvr<n ) then info = -12_${ik}$ else ! set m to the number of eigenpairs for which condition numbers ! are required, and test mm. if( somcon ) then m = 0_${ik}$ pair = .false. do k = 1, n if( pair ) then pair = .false. else if( k<n ) then if( a( k+1, k )==zero ) then if( select( k ) )m = m + 1_${ik}$ else pair = .true. if( select( k ) .or. select( k+1 ) )m = m + 2_${ik}$ end if else if( select( n ) )m = m + 1_${ik}$ end if end if end do else m = n end if if( n==0_${ik}$ ) then lwmin = 1_${ik}$ else if( stdlib_lsame( job, 'V' ) .or. stdlib_lsame( job, 'B' ) ) then lwmin = 2_${ik}$*n*( n + 2_${ik}$ ) + 16_${ik}$ else lwmin = n end if work( 1_${ik}$ ) = lwmin if( mm<m ) then info = -15_${ik}$ else if( lwork<lwmin .and. .not.lquery ) then info = -18_${ik}$ end if end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DTGSNA', -info ) return else if( lquery ) then return end if ! quick return if possible if( n==0 )return ! get machine constants eps = stdlib${ii}$_${ri}$lamch( 'P' ) smlnum = stdlib${ii}$_${ri}$lamch( 'S' ) / eps ks = 0_${ik}$ pair = .false. loop_20: do k = 1, n ! determine whether a(k,k) begins a 1-by-1 or 2-by-2 block. if( pair ) then pair = .false. cycle loop_20 else if( k<n )pair = a( k+1, k )/=zero end if ! determine whether condition numbers are required for the k-th ! eigenpair. if( somcon ) then if( pair ) then if( .not.select( k ) .and. .not.select( k+1 ) )cycle loop_20 else if( .not.select( k ) )cycle loop_20 end if end if ks = ks + 1_${ik}$ if( wants ) then ! compute the reciprocal condition number of the k-th ! eigenvalue. if( pair ) then ! complex eigenvalue pair. rnrm = stdlib${ii}$_${ri}$lapy2( stdlib${ii}$_${ri}$nrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_${ri}$nrm2( n, vr( & 1_${ik}$, ks+1 ), 1_${ik}$ ) ) lnrm = stdlib${ii}$_${ri}$lapy2( stdlib${ii}$_${ri}$nrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ),stdlib${ii}$_${ri}$nrm2( n, vl( & 1_${ik}$, ks+1 ), 1_${ik}$ ) ) call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) tmprr = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) tmpri = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ ) tmpii = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) tmpir = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) uhav = tmprr + tmpii uhavi = tmpir - tmpri call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) tmprr = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) tmpri = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks+1 ), 1_${ik}$,zero, work, 1_${ik}$ ) tmpii = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks+1 ), 1_${ik}$ ) tmpir = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) uhbv = tmprr + tmpii uhbvi = tmpir - tmpri uhav = stdlib${ii}$_${ri}$lapy2( uhav, uhavi ) uhbv = stdlib${ii}$_${ri}$lapy2( uhbv, uhbvi ) cond = stdlib${ii}$_${ri}$lapy2( uhav, uhbv ) s( ks ) = cond / ( rnrm*lnrm ) s( ks+1 ) = s( ks ) else ! real eigenvalue. rnrm = stdlib${ii}$_${ri}$nrm2( n, vr( 1_${ik}$, ks ), 1_${ik}$ ) lnrm = stdlib${ii}$_${ri}$nrm2( n, vl( 1_${ik}$, ks ), 1_${ik}$ ) call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, a, lda, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) uhav = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) call stdlib${ii}$_${ri}$gemv( 'N', n, n, one, b, ldb, vr( 1_${ik}$, ks ), 1_${ik}$, zero,work, 1_${ik}$ ) uhbv = stdlib${ii}$_${ri}$dot( n, work, 1_${ik}$, vl( 1_${ik}$, ks ), 1_${ik}$ ) cond = stdlib${ii}$_${ri}$lapy2( uhav, uhbv ) if( cond==zero ) then s( ks ) = -one else s( ks ) = cond / ( rnrm*lnrm ) end if end if end if if( wantdf ) then if( n==1_${ik}$ ) then dif( ks ) = stdlib${ii}$_${ri}$lapy2( a( 1_${ik}$, 1_${ik}$ ), b( 1_${ik}$, 1_${ik}$ ) ) cycle loop_20 end if ! estimate the reciprocal condition number of the k-th ! eigenvectors. if( pair ) then ! copy the 2-by 2 pencil beginning at (a(k,k), b(k, k)). ! compute the eigenvalue(s) at position k. work( 1_${ik}$ ) = a( k, k ) work( 2_${ik}$ ) = a( k+1, k ) work( 3_${ik}$ ) = a( k, k