#:include "common.fypp" submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_sym_comp implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_ssygst( itype, uplo, n, a, lda, b, ldb, info ) !! SSYGST reduces a real symmetric-definite generalized eigenproblem !! to standard form. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. !! B must have been previously factorized as U**T*U or L*L**T by SPOTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, lda, ldb, n ! Array Arguments real(sp), intent(inout) :: a(lda,*) real(sp), intent(in) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: k, kb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSYGST', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSYGST', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code call stdlib${ii}$_ssygs2( itype, uplo, n, a, lda, b, ldb, info ) else ! use blocked code if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) do k = 1, n, nb kb = min( n-k+1, nb ) ! update the upper triangle of a(k:n,k:n) call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) if( k+kb<=n ) then call stdlib${ii}$_strsm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT',kb, n-k-kb+1, & one, b( k, k ), ldb,a( k, k+kb ), lda ) call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( & k, k+kb ), ldb, one,a( k, k+kb ), lda ) call stdlib${ii}$_ssyr2k( uplo, 'TRANSPOSE', n-k-kb+1, kb, -one,a( k, k+kb ), & lda, b( k, k+kb ), ldb,one, a( k+kb, k+kb ), lda ) call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( & k, k+kb ), ldb, one,a( k, k+kb ), lda ) call stdlib${ii}$_strsm( 'RIGHT', uplo, 'NO TRANSPOSE','NON-UNIT', kb, n-k-kb+& 1_${ik}$, one,b( k+kb, k+kb ), ldb, a( k, k+kb ),lda ) end if end do else ! compute inv(l)*a*inv(l**t) do k = 1, n, nb kb = min( n-k+1, nb ) ! update the lower triangle of a(k:n,k:n) call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) if( k+kb<=n ) then call stdlib${ii}$_strsm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',n-k-kb+1, kb, & one, b( k, k ), ldb,a( k+kb, k ), lda ) call stdlib${ii}$_ssymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(& k+kb, k ), ldb, one,a( k+kb, k ), lda ) call stdlib${ii}$_ssyr2k( uplo, 'NO TRANSPOSE', n-k-kb+1, kb,-one, a( k+kb, k & ), lda, b( k+kb, k ),ldb, one, a( k+kb, k+kb ), lda ) call stdlib${ii}$_ssymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(& k+kb, k ), ldb, one,a( k+kb, k ), lda ) call stdlib${ii}$_strsm( 'LEFT', uplo, 'NO TRANSPOSE','NON-UNIT', n-k-kb+1, & kb, one,b( k+kb, k+kb ), ldb, a( k+kb, k ),lda ) end if end do end if else if( upper ) then ! compute u*a*u**t do k = 1, n, nb kb = min( n-k+1, nb ) ! update the upper triangle of a(1:k+kb-1,1:k+kb-1) call stdlib${ii}$_strmm( 'LEFT', uplo, 'NO TRANSPOSE', 'NON-UNIT',k-1, kb, one, & b, ldb, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_ssymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), & ldb, one, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_ssyr2k( uplo, 'NO TRANSPOSE', k-1, kb, one,a( 1_${ik}$, k ), lda, b( & 1_${ik}$, k ), ldb, one, a,lda ) call stdlib${ii}$_ssymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), & ldb, one, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_strmm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',k-1, kb, one, b( & k, k ), ldb, a( 1_${ik}$, k ),lda ) call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) end do else ! compute l**t*a*l do k = 1, n, nb kb = min( n-k+1, nb ) ! update the lower triangle of a(1:k+kb-1,1:k+kb-1) call stdlib${ii}$_strmm( 'RIGHT', uplo, 'NO TRANSPOSE', 'NON-UNIT',kb, k-1, one, & b, ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), & ldb, one, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_ssyr2k( uplo, 'TRANSPOSE', k-1, kb, one,a( k, 1_${ik}$ ), lda, b( k, & 1_${ik}$ ), ldb, one, a,lda ) call stdlib${ii}$_ssymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), & ldb, one, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_strmm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT', kb,k-1, one, b( & k, k ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_ssygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) end do end if end if end if return end subroutine stdlib${ii}$_ssygst pure module subroutine stdlib${ii}$_dsygst( itype, uplo, n, a, lda, b, ldb, info ) !! DSYGST reduces a real symmetric-definite generalized eigenproblem !! to standard form. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. !! B must have been previously factorized as U**T*U or L*L**T by DPOTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, lda, ldb, n ! Array Arguments real(dp), intent(inout) :: a(lda,*) real(dp), intent(in) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: k, kb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSYGST', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYGST', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code call stdlib${ii}$_dsygs2( itype, uplo, n, a, lda, b, ldb, info ) else ! use blocked code if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) do k = 1, n, nb kb = min( n-k+1, nb ) ! update the upper triangle of a(k:n,k:n) call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) if( k+kb<=n ) then call stdlib${ii}$_dtrsm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT',kb, n-k-kb+1, & one, b( k, k ), ldb,a( k, k+kb ), lda ) call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( & k, k+kb ), ldb, one,a( k, k+kb ), lda ) call stdlib${ii}$_dsyr2k( uplo, 'TRANSPOSE', n-k-kb+1, kb, -one,a( k, k+kb ), & lda, b( k, k+kb ), ldb,one, a( k+kb, k+kb ), lda ) call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( & k, k+kb ), ldb, one,a( k, k+kb ), lda ) call stdlib${ii}$_dtrsm( 'RIGHT', uplo, 'NO TRANSPOSE','NON-UNIT', kb, n-k-kb+& 1_${ik}$, one,b( k+kb, k+kb ), ldb, a( k, k+kb ),lda ) end if end do else ! compute inv(l)*a*inv(l**t) do k = 1, n, nb kb = min( n-k+1, nb ) ! update the lower triangle of a(k:n,k:n) call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) if( k+kb<=n ) then call stdlib${ii}$_dtrsm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',n-k-kb+1, kb, & one, b( k, k ), ldb,a( k+kb, k ), lda ) call stdlib${ii}$_dsymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(& k+kb, k ), ldb, one,a( k+kb, k ), lda ) call stdlib${ii}$_dsyr2k( uplo, 'NO TRANSPOSE', n-k-kb+1, kb,-one, a( k+kb, k & ), lda, b( k+kb, k ),ldb, one, a( k+kb, k+kb ), lda ) call stdlib${ii}$_dsymm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(& k+kb, k ), ldb, one,a( k+kb, k ), lda ) call stdlib${ii}$_dtrsm( 'LEFT', uplo, 'NO TRANSPOSE','NON-UNIT', n-k-kb+1, & kb, one,b( k+kb, k+kb ), ldb, a( k+kb, k ),lda ) end if end do end if else if( upper ) then ! compute u*a*u**t do k = 1, n, nb kb = min( n-k+1, nb ) ! update the upper triangle of a(1:k+kb-1,1:k+kb-1) call stdlib${ii}$_dtrmm( 'LEFT', uplo, 'NO TRANSPOSE', 'NON-UNIT',k-1, kb, one, & b, ldb, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_dsymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), & ldb, one, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_dsyr2k( uplo, 'NO TRANSPOSE', k-1, kb, one,a( 1_${ik}$, k ), lda, b( & 1_${ik}$, k ), ldb, one, a,lda ) call stdlib${ii}$_dsymm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), & ldb, one, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_dtrmm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',k-1, kb, one, b( & k, k ), ldb, a( 1_${ik}$, k ),lda ) call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) end do else ! compute l**t*a*l do k = 1, n, nb kb = min( n-k+1, nb ) ! update the lower triangle of a(1:k+kb-1,1:k+kb-1) call stdlib${ii}$_dtrmm( 'RIGHT', uplo, 'NO TRANSPOSE', 'NON-UNIT',kb, k-1, one, & b, ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), & ldb, one, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_dsyr2k( uplo, 'TRANSPOSE', k-1, kb, one,a( k, 1_${ik}$ ), lda, b( k, & 1_${ik}$ ), ldb, one, a,lda ) call stdlib${ii}$_dsymm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), & ldb, one, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_dtrmm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT', kb,k-1, one, b( & k, k ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_dsygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) end do end if end if end if return end subroutine stdlib${ii}$_dsygst #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$sygst( itype, uplo, n, a, lda, b, ldb, info ) !! DSYGST: reduces a real symmetric-definite generalized eigenproblem !! to standard form. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. !! B must have been previously factorized as U**T*U or L*L**T by DPOTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, lda, ldb, n ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(in) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: k, kb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSYGST', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYGST', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code call stdlib${ii}$_${ri}$sygs2( itype, uplo, n, a, lda, b, ldb, info ) else ! use blocked code if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) do k = 1, n, nb kb = min( n-k+1, nb ) ! update the upper triangle of a(k:n,k:n) call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) if( k+kb<=n ) then call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT',kb, n-k-kb+1, & one, b( k, k ), ldb,a( k, k+kb ), lda ) call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( & k, k+kb ), ldb, one,a( k, k+kb ), lda ) call stdlib${ii}$_${ri}$syr2k( uplo, 'TRANSPOSE', n-k-kb+1, kb, -one,a( k, k+kb ), & lda, b( k, k+kb ), ldb,one, a( k+kb, k+kb ), lda ) call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, n-k-kb+1, -half,a( k, k ), lda, b( & k, k+kb ), ldb, one,a( k, k+kb ), lda ) call stdlib${ii}$_${ri}$trsm( 'RIGHT', uplo, 'NO TRANSPOSE','NON-UNIT', kb, n-k-kb+& 1_${ik}$, one,b( k+kb, k+kb ), ldb, a( k, k+kb ),lda ) end if end do else ! compute inv(l)*a*inv(l**t) do k = 1, n, nb kb = min( n-k+1, nb ) ! update the lower triangle of a(k:n,k:n) call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) if( k+kb<=n ) then call stdlib${ii}$_${ri}$trsm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',n-k-kb+1, kb, & one, b( k, k ), ldb,a( k+kb, k ), lda ) call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(& k+kb, k ), ldb, one,a( k+kb, k ), lda ) call stdlib${ii}$_${ri}$syr2k( uplo, 'NO TRANSPOSE', n-k-kb+1, kb,-one, a( k+kb, k & ), lda, b( k+kb, k ),ldb, one, a( k+kb, k+kb ), lda ) call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, n-k-kb+1, kb, -half,a( k, k ), lda, b(& k+kb, k ), ldb, one,a( k+kb, k ), lda ) call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, 'NO TRANSPOSE','NON-UNIT', n-k-kb+1, & kb, one,b( k+kb, k+kb ), ldb, a( k+kb, k ),lda ) end if end do end if else if( upper ) then ! compute u*a*u**t do k = 1, n, nb kb = min( n-k+1, nb ) ! update the upper triangle of a(1:k+kb-1,1:k+kb-1) call stdlib${ii}$_${ri}$trmm( 'LEFT', uplo, 'NO TRANSPOSE', 'NON-UNIT',k-1, kb, one, & b, ldb, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), & ldb, one, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_${ri}$syr2k( uplo, 'NO TRANSPOSE', k-1, kb, one,a( 1_${ik}$, k ), lda, b( & 1_${ik}$, k ), ldb, one, a,lda ) call stdlib${ii}$_${ri}$symm( 'RIGHT', uplo, k-1, kb, half, a( k, k ),lda, b( 1_${ik}$, k ), & ldb, one, a( 1_${ik}$, k ), lda ) call stdlib${ii}$_${ri}$trmm( 'RIGHT', uplo, 'TRANSPOSE', 'NON-UNIT',k-1, kb, one, b( & k, k ), ldb, a( 1_${ik}$, k ),lda ) call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) end do else ! compute l**t*a*l do k = 1, n, nb kb = min( n-k+1, nb ) ! update the lower triangle of a(1:k+kb-1,1:k+kb-1) call stdlib${ii}$_${ri}$trmm( 'RIGHT', uplo, 'NO TRANSPOSE', 'NON-UNIT',kb, k-1, one, & b, ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), & ldb, one, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_${ri}$syr2k( uplo, 'TRANSPOSE', k-1, kb, one,a( k, 1_${ik}$ ), lda, b( k, & 1_${ik}$ ), ldb, one, a,lda ) call stdlib${ii}$_${ri}$symm( 'LEFT', uplo, kb, k-1, half, a( k, k ),lda, b( k, 1_${ik}$ ), & ldb, one, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_${ri}$trmm( 'LEFT', uplo, 'TRANSPOSE', 'NON-UNIT', kb,k-1, one, b( & k, k ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_${ri}$sygs2( itype, uplo, kb, a( k, k ), lda,b( k, k ), ldb, info ) end do end if end if end if return end subroutine stdlib${ii}$_${ri}$sygst #:endif #:endfor pure module subroutine stdlib${ii}$_ssygs2( itype, uplo, n, a, lda, b, ldb, info ) !! SSYGS2 reduces a real symmetric-definite generalized eigenproblem !! to standard form. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T *A*L. !! B must have been previously factorized as U**T *U or L*L**T by SPOTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, lda, ldb, n ! Array Arguments real(sp), intent(inout) :: a(lda,*) real(sp), intent(in) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: k real(sp) :: akk, bkk, ct ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSYGS2', -info ) return end if if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) do k = 1, n ! update the upper triangle of a(k:n,k:n) akk = a( k, k ) bkk = b( k, k ) akk = akk / bkk**2_${ik}$ a( k, k ) = akk if( k<n ) then call stdlib${ii}$_sscal( n-k, one / bkk, a( k, k+1 ), lda ) ct = -half*akk call stdlib${ii}$_saxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda ) call stdlib${ii}$_ssyr2( uplo, n-k, -one, a( k, k+1 ), lda,b( k, k+1 ), ldb, a( & k+1, k+1 ), lda ) call stdlib${ii}$_saxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda ) call stdlib${ii}$_strsv( uplo, 'TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), ldb, & a( k, k+1 ), lda ) end if end do else ! compute inv(l)*a*inv(l**t) do k = 1, n ! update the lower triangle of a(k:n,k:n) akk = a( k, k ) bkk = b( k, k ) akk = akk / bkk**2_${ik}$ a( k, k ) = akk if( k<n ) then call stdlib${ii}$_sscal( n-k, one / bkk, a( k+1, k ), 1_${ik}$ ) ct = -half*akk call stdlib${ii}$_saxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ ) call stdlib${ii}$_ssyr2( uplo, n-k, -one, a( k+1, k ), 1_${ik}$,b( k+1, k ), 1_${ik}$, a( k+1, & k+1 ), lda ) call stdlib${ii}$_saxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ ) call stdlib${ii}$_strsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), & ldb, a( k+1, k ), 1_${ik}$ ) end if end do end if else if( upper ) then ! compute u*a*u**t do k = 1, n ! update the upper triangle of a(1:k,1:k) akk = a( k, k ) bkk = b( k, k ) call stdlib${ii}$_strmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, b,ldb, a( 1_${ik}$, k ), 1_${ik}$ & ) ct = half*akk call stdlib${ii}$_saxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ ) call stdlib${ii}$_ssyr2( uplo, k-1, one, a( 1_${ik}$, k ), 1_${ik}$, b( 1_${ik}$, k ), 1_${ik}$,a, lda ) call stdlib${ii}$_saxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ ) call stdlib${ii}$_sscal( k-1, bkk, a( 1_${ik}$, k ), 1_${ik}$ ) a( k, k ) = akk*bkk**2_${ik}$ end do else ! compute l**t *a*l do k = 1, n ! update the lower triangle of a(1:k,1:k) akk = a( k, k ) bkk = b( k, k ) call stdlib${ii}$_strmv( uplo, 'TRANSPOSE', 'NON-UNIT', k-1, b, ldb,a( k, 1_${ik}$ ), lda ) ct = half*akk call stdlib${ii}$_saxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_ssyr2( uplo, k-1, one, a( k, 1_${ik}$ ), lda, b( k, 1_${ik}$ ),ldb, a, lda ) call stdlib${ii}$_saxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_sscal( k-1, bkk, a( k, 1_${ik}$ ), lda ) a( k, k ) = akk*bkk**2_${ik}$ end do end if end if return end subroutine stdlib${ii}$_ssygs2 pure module subroutine stdlib${ii}$_dsygs2( itype, uplo, n, a, lda, b, ldb, info ) !! DSYGS2 reduces a real symmetric-definite generalized eigenproblem !! to standard form. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T *A*L. !! B must have been previously factorized as U**T *U or L*L**T by DPOTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, lda, ldb, n ! Array Arguments real(dp), intent(inout) :: a(lda,*) real(dp), intent(in) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: k real(dp) :: akk, bkk, ct ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSYGS2', -info ) return end if if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) do k = 1, n ! update the upper triangle of a(k:n,k:n) akk = a( k, k ) bkk = b( k, k ) akk = akk / bkk**2_${ik}$ a( k, k ) = akk if( k<n ) then call stdlib${ii}$_dscal( n-k, one / bkk, a( k, k+1 ), lda ) ct = -half*akk call stdlib${ii}$_daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda ) call stdlib${ii}$_dsyr2( uplo, n-k, -one, a( k, k+1 ), lda,b( k, k+1 ), ldb, a( & k+1, k+1 ), lda ) call stdlib${ii}$_daxpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda ) call stdlib${ii}$_dtrsv( uplo, 'TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), ldb, & a( k, k+1 ), lda ) end if end do else ! compute inv(l)*a*inv(l**t) do k = 1, n ! update the lower triangle of a(k:n,k:n) akk = a( k, k ) bkk = b( k, k ) akk = akk / bkk**2_${ik}$ a( k, k ) = akk if( k<n ) then call stdlib${ii}$_dscal( n-k, one / bkk, a( k+1, k ), 1_${ik}$ ) ct = -half*akk call stdlib${ii}$_daxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ ) call stdlib${ii}$_dsyr2( uplo, n-k, -one, a( k+1, k ), 1_${ik}$,b( k+1, k ), 1_${ik}$, a( k+1, & k+1 ), lda ) call stdlib${ii}$_daxpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ ) call stdlib${ii}$_dtrsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), & ldb, a( k+1, k ), 1_${ik}$ ) end if end do end if else if( upper ) then ! compute u*a*u**t do k = 1, n ! update the upper triangle of a(1:k,1:k) akk = a( k, k ) bkk = b( k, k ) call stdlib${ii}$_dtrmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, b,ldb, a( 1_${ik}$, k ), 1_${ik}$ & ) ct = half*akk call stdlib${ii}$_daxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ ) call stdlib${ii}$_dsyr2( uplo, k-1, one, a( 1_${ik}$, k ), 1_${ik}$, b( 1_${ik}$, k ), 1_${ik}$,a, lda ) call stdlib${ii}$_daxpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ ) call stdlib${ii}$_dscal( k-1, bkk, a( 1_${ik}$, k ), 1_${ik}$ ) a( k, k ) = akk*bkk**2_${ik}$ end do else ! compute l**t *a*l do k = 1, n ! update the lower triangle of a(1:k,1:k) akk = a( k, k ) bkk = b( k, k ) call stdlib${ii}$_dtrmv( uplo, 'TRANSPOSE', 'NON-UNIT', k-1, b, ldb,a( k, 1_${ik}$ ), lda ) ct = half*akk call stdlib${ii}$_daxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_dsyr2( uplo, k-1, one, a( k, 1_${ik}$ ), lda, b( k, 1_${ik}$ ),ldb, a, lda ) call stdlib${ii}$_daxpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_dscal( k-1, bkk, a( k, 1_${ik}$ ), lda ) a( k, k ) = akk*bkk**2_${ik}$ end do end if end if return end subroutine stdlib${ii}$_dsygs2 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$sygs2( itype, uplo, n, a, lda, b, ldb, info ) !! DSYGS2: reduces a real symmetric-definite generalized eigenproblem !! to standard form. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T *A*L. !! B must have been previously factorized as U**T *U or L*L**T by DPOTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, lda, ldb, n ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(in) :: b(ldb,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: k real(${rk}$) :: akk, bkk, ct ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -5_${ik}$ else if( ldb<max( 1_${ik}$, n ) ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSYGS2', -info ) return end if if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) do k = 1, n ! update the upper triangle of a(k:n,k:n) akk = a( k, k ) bkk = b( k, k ) akk = akk / bkk**2_${ik}$ a( k, k ) = akk if( k<n ) then call stdlib${ii}$_${ri}$scal( n-k, one / bkk, a( k, k+1 ), lda ) ct = -half*akk call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda ) call stdlib${ii}$_${ri}$syr2( uplo, n-k, -one, a( k, k+1 ), lda,b( k, k+1 ), ldb, a( & k+1, k+1 ), lda ) call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k, k+1 ), ldb, a( k, k+1 ),lda ) call stdlib${ii}$_${ri}$trsv( uplo, 'TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), ldb, & a( k, k+1 ), lda ) end if end do else ! compute inv(l)*a*inv(l**t) do k = 1, n ! update the lower triangle of a(k:n,k:n) akk = a( k, k ) bkk = b( k, k ) akk = akk / bkk**2_${ik}$ a( k, k ) = akk if( k<n ) then call stdlib${ii}$_${ri}$scal( n-k, one / bkk, a( k+1, k ), 1_${ik}$ ) ct = -half*akk call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ ) call stdlib${ii}$_${ri}$syr2( uplo, n-k, -one, a( k+1, k ), 1_${ik}$,b( k+1, k ), 1_${ik}$, a( k+1, & k+1 ), lda ) call stdlib${ii}$_${ri}$axpy( n-k, ct, b( k+1, k ), 1_${ik}$, a( k+1, k ), 1_${ik}$ ) call stdlib${ii}$_${ri}$trsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,b( k+1, k+1 ), & ldb, a( k+1, k ), 1_${ik}$ ) end if end do end if else if( upper ) then ! compute u*a*u**t do k = 1, n ! update the upper triangle of a(1:k,1:k) akk = a( k, k ) bkk = b( k, k ) call stdlib${ii}$_${ri}$trmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, b,ldb, a( 1_${ik}$, k ), 1_${ik}$ & ) ct = half*akk call stdlib${ii}$_${ri}$axpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ ) call stdlib${ii}$_${ri}$syr2( uplo, k-1, one, a( 1_${ik}$, k ), 1_${ik}$, b( 1_${ik}$, k ), 1_${ik}$,a, lda ) call stdlib${ii}$_${ri}$axpy( k-1, ct, b( 1_${ik}$, k ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( k-1, bkk, a( 1_${ik}$, k ), 1_${ik}$ ) a( k, k ) = akk*bkk**2_${ik}$ end do else ! compute l**t *a*l do k = 1, n ! update the lower triangle of a(1:k,1:k) akk = a( k, k ) bkk = b( k, k ) call stdlib${ii}$_${ri}$trmv( uplo, 'TRANSPOSE', 'NON-UNIT', k-1, b, ldb,a( k, 1_${ik}$ ), lda ) ct = half*akk call stdlib${ii}$_${ri}$axpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_${ri}$syr2( uplo, k-1, one, a( k, 1_${ik}$ ), lda, b( k, 1_${ik}$ ),ldb, a, lda ) call stdlib${ii}$_${ri}$axpy( k-1, ct, b( k, 1_${ik}$ ), ldb, a( k, 1_${ik}$ ), lda ) call stdlib${ii}$_${ri}$scal( k-1, bkk, a( k, 1_${ik}$ ), lda ) a( k, k ) = akk*bkk**2_${ik}$ end do end if end if return end subroutine stdlib${ii}$_${ri}$sygs2 #:endif #:endfor pure module subroutine stdlib${ii}$_sspgst( itype, uplo, n, ap, bp, info ) !! SSPGST reduces a real symmetric-definite generalized eigenproblem !! to standard form, using packed storage. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. !! B must have been previously factorized as U**T*U or L*L**T by SPPTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, n ! Array Arguments real(sp), intent(inout) :: ap(*) real(sp), intent(in) :: bp(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, j1, j1j1, jj, k, k1, k1k1, kk real(sp) :: ajj, akk, bjj, bkk, ct ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSPGST', -info ) return end if if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) ! j1 and jj are the indices of a(1,j) and a(j,j) jj = 0_${ik}$ do j = 1, n j1 = jj + 1_${ik}$ jj = jj + j ! compute the j-th column of the upper triangle of a bjj = bp( jj ) call stdlib${ii}$_stpsv( uplo, 'TRANSPOSE', 'NONUNIT', j, bp,ap( j1 ), 1_${ik}$ ) call stdlib${ii}$_sspmv( uplo, j-1, -one, ap, bp( j1 ), 1_${ik}$, one,ap( j1 ), 1_${ik}$ ) call stdlib${ii}$_sscal( j-1, one / bjj, ap( j1 ), 1_${ik}$ ) ap( jj ) = ( ap( jj )-stdlib${ii}$_sdot( j-1, ap( j1 ), 1_${ik}$, bp( j1 ),1_${ik}$ ) ) / & bjj end do else ! compute inv(l)*a*inv(l**t) ! kk and k1k1 are the indices of a(k,k) and a(k+1,k+1) kk = 1_${ik}$ do k = 1, n k1k1 = kk + n - k + 1_${ik}$ ! update the lower triangle of a(k:n,k:n) akk = ap( kk ) bkk = bp( kk ) akk = akk / bkk**2_${ik}$ ap( kk ) = akk if( k<n ) then call stdlib${ii}$_sscal( n-k, one / bkk, ap( kk+1 ), 1_${ik}$ ) ct = -half*akk call stdlib${ii}$_saxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ ) call stdlib${ii}$_sspr2( uplo, n-k, -one, ap( kk+1 ), 1_${ik}$,bp( kk+1 ), 1_${ik}$, ap( k1k1 )& ) call stdlib${ii}$_saxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ ) call stdlib${ii}$_stpsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,bp( k1k1 ), ap( & kk+1 ), 1_${ik}$ ) end if kk = k1k1 end do end if else if( upper ) then ! compute u*a*u**t ! k1 and kk are the indices of a(1,k) and a(k,k) kk = 0_${ik}$ do k = 1, n k1 = kk + 1_${ik}$ kk = kk + k ! update the upper triangle of a(1:k,1:k) akk = ap( kk ) bkk = bp( kk ) call stdlib${ii}$_stpmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, bp,ap( k1 ), 1_${ik}$ ) ct = half*akk call stdlib${ii}$_saxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ ) call stdlib${ii}$_sspr2( uplo, k-1, one, ap( k1 ), 1_${ik}$, bp( k1 ), 1_${ik}$,ap ) call stdlib${ii}$_saxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ ) call stdlib${ii}$_sscal( k-1, bkk, ap( k1 ), 1_${ik}$ ) ap( kk ) = akk*bkk**2_${ik}$ end do else ! compute l**t *a*l ! jj and j1j1 are the indices of a(j,j) and a(j+1,j+1) jj = 1_${ik}$ do j = 1, n j1j1 = jj + n - j + 1_${ik}$ ! compute the j-th column of the lower triangle of a ajj = ap( jj ) bjj = bp( jj ) ap( jj ) = ajj*bjj + stdlib${ii}$_sdot( n-j, ap( jj+1 ), 1_${ik}$,bp( jj+1 ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-j, bjj, ap( jj+1 ), 1_${ik}$ ) call stdlib${ii}$_sspmv( uplo, n-j, one, ap( j1j1 ), bp( jj+1 ), 1_${ik}$,one, ap( jj+1 ), & 1_${ik}$ ) call stdlib${ii}$_stpmv( uplo, 'TRANSPOSE', 'NON-UNIT', n-j+1,bp( jj ), ap( jj ), 1_${ik}$ & ) jj = j1j1 end do end if end if return end subroutine stdlib${ii}$_sspgst pure module subroutine stdlib${ii}$_dspgst( itype, uplo, n, ap, bp, info ) !! DSPGST reduces a real symmetric-definite generalized eigenproblem !! to standard form, using packed storage. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. !! B must have been previously factorized as U**T*U or L*L**T by DPPTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, n ! Array Arguments real(dp), intent(inout) :: ap(*) real(dp), intent(in) :: bp(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, j1, j1j1, jj, k, k1, k1k1, kk real(dp) :: ajj, akk, bjj, bkk, ct ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSPGST', -info ) return end if if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) ! j1 and jj are the indices of a(1,j) and a(j,j) jj = 0_${ik}$ do j = 1, n j1 = jj + 1_${ik}$ jj = jj + j ! compute the j-th column of the upper triangle of a bjj = bp( jj ) call stdlib${ii}$_dtpsv( uplo, 'TRANSPOSE', 'NONUNIT', j, bp,ap( j1 ), 1_${ik}$ ) call stdlib${ii}$_dspmv( uplo, j-1, -one, ap, bp( j1 ), 1_${ik}$, one,ap( j1 ), 1_${ik}$ ) call stdlib${ii}$_dscal( j-1, one / bjj, ap( j1 ), 1_${ik}$ ) ap( jj ) = ( ap( jj )-stdlib${ii}$_ddot( j-1, ap( j1 ), 1_${ik}$, bp( j1 ),1_${ik}$ ) ) / & bjj end do else ! compute inv(l)*a*inv(l**t) ! kk and k1k1 are the indices of a(k,k) and a(k+1,k+1) kk = 1_${ik}$ do k = 1, n k1k1 = kk + n - k + 1_${ik}$ ! update the lower triangle of a(k:n,k:n) akk = ap( kk ) bkk = bp( kk ) akk = akk / bkk**2_${ik}$ ap( kk ) = akk if( k<n ) then call stdlib${ii}$_dscal( n-k, one / bkk, ap( kk+1 ), 1_${ik}$ ) ct = -half*akk call stdlib${ii}$_daxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ ) call stdlib${ii}$_dspr2( uplo, n-k, -one, ap( kk+1 ), 1_${ik}$,bp( kk+1 ), 1_${ik}$, ap( k1k1 )& ) call stdlib${ii}$_daxpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ ) call stdlib${ii}$_dtpsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,bp( k1k1 ), ap( & kk+1 ), 1_${ik}$ ) end if kk = k1k1 end do end if else if( upper ) then ! compute u*a*u**t ! k1 and kk are the indices of a(1,k) and a(k,k) kk = 0_${ik}$ do k = 1, n k1 = kk + 1_${ik}$ kk = kk + k ! update the upper triangle of a(1:k,1:k) akk = ap( kk ) bkk = bp( kk ) call stdlib${ii}$_dtpmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, bp,ap( k1 ), 1_${ik}$ ) ct = half*akk call stdlib${ii}$_daxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ ) call stdlib${ii}$_dspr2( uplo, k-1, one, ap( k1 ), 1_${ik}$, bp( k1 ), 1_${ik}$,ap ) call stdlib${ii}$_daxpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ ) call stdlib${ii}$_dscal( k-1, bkk, ap( k1 ), 1_${ik}$ ) ap( kk ) = akk*bkk**2_${ik}$ end do else ! compute l**t *a*l ! jj and j1j1 are the indices of a(j,j) and a(j+1,j+1) jj = 1_${ik}$ do j = 1, n j1j1 = jj + n - j + 1_${ik}$ ! compute the j-th column of the lower triangle of a ajj = ap( jj ) bjj = bp( jj ) ap( jj ) = ajj*bjj + stdlib${ii}$_ddot( n-j, ap( jj+1 ), 1_${ik}$,bp( jj+1 ), 1_${ik}$ ) call stdlib${ii}$_dscal( n-j, bjj, ap( jj+1 ), 1_${ik}$ ) call stdlib${ii}$_dspmv( uplo, n-j, one, ap( j1j1 ), bp( jj+1 ), 1_${ik}$,one, ap( jj+1 ), & 1_${ik}$ ) call stdlib${ii}$_dtpmv( uplo, 'TRANSPOSE', 'NON-UNIT', n-j+1,bp( jj ), ap( jj ), 1_${ik}$ & ) jj = j1j1 end do end if end if return end subroutine stdlib${ii}$_dspgst #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$spgst( itype, uplo, n, ap, bp, info ) !! DSPGST: reduces a real symmetric-definite generalized eigenproblem !! to standard form, using packed storage. !! If ITYPE = 1, the problem is A*x = lambda*B*x, !! and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) !! If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or !! B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. !! B must have been previously factorized as U**T*U or L*L**T by DPPTRF. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: itype, n ! Array Arguments real(${rk}$), intent(inout) :: ap(*) real(${rk}$), intent(in) :: bp(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, j1, j1j1, jj, k, k1, k1k1, kk real(${rk}$) :: ajj, akk, bjj, bkk, ct ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( itype<1_${ik}$ .or. itype>3_${ik}$ ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSPGST', -info ) return end if if( itype==1_${ik}$ ) then if( upper ) then ! compute inv(u**t)*a*inv(u) ! j1 and jj are the indices of a(1,j) and a(j,j) jj = 0_${ik}$ do j = 1, n j1 = jj + 1_${ik}$ jj = jj + j ! compute the j-th column of the upper triangle of a bjj = bp( jj ) call stdlib${ii}$_${ri}$tpsv( uplo, 'TRANSPOSE', 'NONUNIT', j, bp,ap( j1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$spmv( uplo, j-1, -one, ap, bp( j1 ), 1_${ik}$, one,ap( j1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( j-1, one / bjj, ap( j1 ), 1_${ik}$ ) ap( jj ) = ( ap( jj )-stdlib${ii}$_${ri}$dot( j-1, ap( j1 ), 1_${ik}$, bp( j1 ),1_${ik}$ ) ) / & bjj end do else ! compute inv(l)*a*inv(l**t) ! kk and k1k1 are the indices of a(k,k) and a(k+1,k+1) kk = 1_${ik}$ do k = 1, n k1k1 = kk + n - k + 1_${ik}$ ! update the lower triangle of a(k:n,k:n) akk = ap( kk ) bkk = bp( kk ) akk = akk / bkk**2_${ik}$ ap( kk ) = akk if( k<n ) then call stdlib${ii}$_${ri}$scal( n-k, one / bkk, ap( kk+1 ), 1_${ik}$ ) ct = -half*akk call stdlib${ii}$_${ri}$axpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$spr2( uplo, n-k, -one, ap( kk+1 ), 1_${ik}$,bp( kk+1 ), 1_${ik}$, ap( k1k1 )& ) call stdlib${ii}$_${ri}$axpy( n-k, ct, bp( kk+1 ), 1_${ik}$, ap( kk+1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$tpsv( uplo, 'NO TRANSPOSE', 'NON-UNIT', n-k,bp( k1k1 ), ap( & kk+1 ), 1_${ik}$ ) end if kk = k1k1 end do end if else if( upper ) then ! compute u*a*u**t ! k1 and kk are the indices of a(1,k) and a(k,k) kk = 0_${ik}$ do k = 1, n k1 = kk + 1_${ik}$ kk = kk + k ! update the upper triangle of a(1:k,1:k) akk = ap( kk ) bkk = bp( kk ) call stdlib${ii}$_${ri}$tpmv( uplo, 'NO TRANSPOSE', 'NON-UNIT', k-1, bp,ap( k1 ), 1_${ik}$ ) ct = half*akk call stdlib${ii}$_${ri}$axpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$spr2( uplo, k-1, one, ap( k1 ), 1_${ik}$, bp( k1 ), 1_${ik}$,ap ) call stdlib${ii}$_${ri}$axpy( k-1, ct, bp( k1 ), 1_${ik}$, ap( k1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( k-1, bkk, ap( k1 ), 1_${ik}$ ) ap( kk ) = akk*bkk**2_${ik}$ end do else ! compute l**t *a*l ! jj and j1j1 are the indices of a(j,j) and a(j+1,j+1) jj = 1_${ik}$ do j = 1, n j1j1 = jj + n - j + 1_${ik}$ ! compute the j-th column of the lower triangle of a ajj = ap( jj ) bjj = bp( jj ) ap( jj ) = ajj*bjj + stdlib${ii}$_${ri}$dot( n-j, ap( jj+1 ), 1_${ik}$,bp( jj+1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( n-j, bjj, ap( jj+1 ), 1_${ik}$ ) call stdlib${ii}$_${ri}$spmv( uplo, n-j, one, ap( j1j1 ), bp( jj+1 ), 1_${ik}$,one, ap( jj+1 ), & 1_${ik}$ ) call stdlib${ii}$_${ri}$tpmv( uplo, 'TRANSPOSE', 'NON-UNIT', n-j+1,bp( jj ), ap( jj ), 1_${ik}$ & ) jj = j1j1 end do end if end if return end subroutine stdlib${ii}$_${ri}$spgst #:endif #:endfor pure module subroutine stdlib${ii}$_ssbgst( vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x,ldx, work, info ) !! SSBGST reduces a real symmetric-definite banded generalized !! eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y, !! such that C has the same bandwidth as A. !! B must have been previously factorized as S**T*S by SPBSTF, using a !! split Cholesky factorization. A is overwritten by C = X**T*A*X, where !! X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the !! bandwidth of A. ! -- 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) :: uplo, vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldx, n ! Array Arguments real(sp), intent(inout) :: ab(ldab,*) real(sp), intent(in) :: bb(ldbb,*) real(sp), intent(out) :: work(*), x(ldx,*) ! ===================================================================== ! Local Scalars logical(lk) :: update, upper, wantx integer(${ik}$) :: i, i0, i1, i2, inca, j, j1, j1t, j2, j2t, k, ka1, kb1, kbt, l, m, nr, & nrt, nx real(sp) :: bii, ra, ra1, t ! Intrinsic Functions ! Executable Statements ! test the input parameters wantx = stdlib_lsame( vect, 'V' ) upper = stdlib_lsame( uplo, 'U' ) ka1 = ka + 1_${ik}$ kb1 = kb + 1_${ik}$ info = 0_${ik}$ if( .not.wantx .and. .not.stdlib_lsame( vect, 'N' ) ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ka<0_${ik}$ ) then info = -4_${ik}$ else if( kb<0_${ik}$ .or. kb>ka ) then info = -5_${ik}$ else if( ldab<ka+1 ) then info = -7_${ik}$ else if( ldbb<kb+1 ) then info = -9_${ik}$ else if( ldx<1_${ik}$ .or. wantx .and. ldx<max( 1_${ik}$, n ) ) then info = -11_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SSBGST', -info ) return end if ! quick return if possible if( n==0 )return inca = ldab*ka1 ! initialize x to the unit matrix, if needed if( wantx )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, x, ldx ) ! set m to the splitting point m. it must be the same value as is ! used in stdlib${ii}$_spbstf. the chosen value allows the arrays work and rwork ! to be of dimension (n). m = ( n+kb ) / 2_${ik}$ ! the routine works in two phases, corresponding to the two halves ! of the split cholesky factorization of b as s**t*s where ! s = ( u ) ! ( m l ) ! with u upper triangular of order m, and l lower triangular of ! order n-m. s has the same bandwidth as b. ! s is treated as a product of elementary matrices: ! s = s(m)*s(m-1)*...*s(2)*s(1)*s(m+1)*s(m+2)*...*s(n-1)*s(n) ! where s(i) is determined by the i-th row of s. ! in phase 1, the index i takes the values n, n-1, ... , m+1; ! in phase 2, it takes the values 1, 2, ... , m. ! for each value of i, the current matrix a is updated by forming ! inv(s(i))**t*a*inv(s(i)). this creates a triangular bulge outside ! the band of a. the bulge is then pushed down toward the bottom of ! a in phase 1, and up toward the top of a in phase 2, by applying ! plane rotations. ! there are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1 ! of them are linearly independent, so annihilating a bulge requires ! only 2*kb-1 plane rotations. the rotations are divided into a 1st ! set of kb-1 rotations, and a 2nd set of kb rotations. ! wherever possible, rotations are generated and applied in vector ! operations of length nr between the indices j1 and j2 (sometimes ! replaced by modified values nrt, j1t or j2t). ! the cosines and sines of the rotations are stored in the array ! work. the cosines of the 1st set of rotations are stored in ! elements n+2:n+m-kb-1 and the sines of the 1st set in elements ! 2:m-kb-1; the cosines of the 2nd set are stored in elements ! n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n. ! the bulges are not formed explicitly; nonzero elements outside the ! band are created only when they are required for generating new ! rotations; they are stored in the array work, in positions where ! they are later overwritten by the sines of the rotations which ! annihilate them. ! **************************** phase 1 ***************************** ! the logical structure of this phase is: ! update = .true. ! do i = n, m + 1, -1 ! use s(i) to update a and create a new bulge ! apply rotations to push all bulges ka positions downward ! end do ! update = .false. ! do i = m + ka + 1, n - 1 ! apply rotations to push all bulges ka positions downward ! end do ! to avoid duplicating code, the two loops are merged. update = .true. i = n + 1_${ik}$ 10 continue if( update ) then i = i - 1_${ik}$ kbt = min( kb, i-1 ) i0 = i - 1_${ik}$ i1 = min( n, i+ka ) i2 = i - kbt + ka1 if( i<m+1 ) then update = .false. i = i + 1_${ik}$ i0 = m if( ka==0 )go to 480 go to 10 end if else i = i + ka if( i>n-1 )go to 480 end if if( upper ) then ! transform a, working with the upper triangle if( update ) then ! form inv(s(i))**t * a * inv(s(i)) bii = bb( kb1, i ) do j = i, i1 ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii end do do j = max( 1, i-ka ), i ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii end do do k = i - kbt, i - 1 do j = i - kbt, k ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( j-i+kb1, i )*ab( k-i+ka1, i ) -bb(& k-i+kb1, i )*ab( j-i+ka1, i ) +ab( ka1, i )*bb( j-i+kb1, i )*bb( k-i+kb1, & i ) end do do j = max( 1, i-ka ), i - kbt - 1 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( k-i+kb1, i )*ab( j-i+ka1, i ) end do end do do j = i, i1 do k = max( j-ka, i-kbt ), i - 1 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( k-i+kb1, i )*ab( i-j+ka1, j ) end do end do if( wantx ) then ! post-multiply x by inv(s(i)) call stdlib${ii}$_sscal( n-m, one / bii, x( m+1, i ), 1_${ik}$ ) if( kbt>0_${ik}$ )call stdlib${ii}$_sger( n-m, kbt, -one, x( m+1, i ), 1_${ik}$,bb( kb1-kbt, i ), & 1_${ik}$, x( m+1, i-kbt ), ldx ) end if ! store a(i,i1) in ra1 for use in next loop over k ra1 = ab( i-i1+ka1, i1 ) end if ! generate and apply vectors of rotations to chase all the ! existing bulges ka positions down toward the bottom of the ! band loop_130: do k = 1, kb - 1 if( update ) then ! determine the rotations which would annihilate the bulge ! which has in theory just been created if( i-k+ka<n .and. i-k>1_${ik}$ ) then ! generate rotation to annihilate a(i,i-k+ka+1) call stdlib${ii}$_slartg( ab( k+1, i-k+ka ), ra1,work( n+i-k+ka-m ), work( i-k+& ka-m ),ra ) ! create nonzero element a(i-k,i-k+ka+1) outside the ! band and store it in work(i-k) t = -bb( kb1-k, i )*ra1 work( i-k ) = work( n+i-k+ka-m )*t -work( i-k+ka-m )*ab( 1_${ik}$, i-k+ka ) ab( 1_${ik}$, i-k+ka ) = work( i-k+ka-m )*t +work( n+i-k+ka-m )*ab( 1_${ik}$, i-k+ka ) ra1 = ra end if end if j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1 nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 if( update ) then j2t = max( j2, i+2*ka-k+1 ) else j2t = j2 end if nrt = ( n-j2t+ka ) / ka1 do j = j2t, j1, ka1 ! create nonzero element a(j-ka,j+1) outside the band ! and store it in work(j-m) work( j-m ) = work( j-m )*ab( 1_${ik}$, j+1 ) ab( 1_${ik}$, j+1 ) = work( n+j-m )*ab( 1_${ik}$, j+1 ) end do ! generate rotations in 1st set to annihilate elements which ! have been created outside the band if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( 1_${ik}$, j2t ), inca, work( j2t-m ), ka1,work( & n+j2t-m ), ka1 ) if( nr>0_${ik}$ ) then ! apply rotations in 1st set from the right do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( ka1-l, j2 ), inca,ab( ka-l, j2+1 ), inca, work(& n+j2-m ),work( j2-m ), ka1 ) end do ! apply rotations in 1st set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),ab( ka, j2+1 ), inca, & work( n+j2-m ),work( j2-m ), ka1 ) end if ! start applying rotations in 1st set from the left do l = ka - 1, kb - k + 1, -1 nrt = ( n-j2+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l & ), inca,work( n+j2-m ), work( j2-m ), ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 1st set do j = j2, j1, ka1 call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j-m ), & work( j-m ) ) end do end if end do loop_130 if( update ) then if( i2<=n .and. kbt>0_${ik}$ ) then ! create nonzero element a(i-kbt,i-kbt+ka+1) outside the ! band and store it in work(i-kbt) work( i-kbt ) = -bb( kb1-kbt, i )*ra1 end if end if loop_170: do k = kb, 1, -1 if( update ) then j2 = i - k - 1_${ik}$ + max( 2_${ik}$, k-i0+1 )*ka1 else j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1 end if ! finish applying rotations in 2nd set from the left do l = kb - k, 1, -1 nrt = ( n-j2+ka+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2-l+1 ), inca,ab( l+1, j2-l+1 ), & inca, work( n+j2-ka ),work( j2-ka ), ka1 ) end do nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 do j = j1, j2, -ka1 work( j ) = work( j-ka ) work( n+j ) = work( n+j-ka ) end do do j = j2, j1, ka1 ! create nonzero element a(j-ka,j+1) outside the band ! and store it in work(j) work( j ) = work( j )*ab( 1_${ik}$, j+1 ) ab( 1_${ik}$, j+1 ) = work( n+j )*ab( 1_${ik}$, j+1 ) end do if( update ) then if( i-k<n-ka .and. k<=kbt )work( i-k+ka ) = work( i-k ) end if end do loop_170 loop_210: do k = kb, 1, -1 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1 nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 if( nr>0_${ik}$ ) then ! generate rotations in 2nd set to annihilate elements ! which have been created outside the band call stdlib${ii}$_slargv( nr, ab( 1_${ik}$, j2 ), inca, work( j2 ), ka1,work( n+j2 ), ka1 ) ! apply rotations in 2nd set from the right do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( ka1-l, j2 ), inca,ab( ka-l, j2+1 ), inca, work(& n+j2 ),work( j2 ), ka1 ) end do ! apply rotations in 2nd set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),ab( ka, j2+1 ), inca, & work( n+j2 ),work( j2 ), ka1 ) end if ! start applying rotations in 2nd set from the left do l = ka - 1, kb - k + 1, -1 nrt = ( n-j2+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l & ), inca, work( n+j2 ),work( j2 ), ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 2nd set do j = j2, j1, ka1 call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j ), work( & j ) ) end do end if end do loop_210 do k = 1, kb - 1 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1 ! finish applying rotations in 1st set from the left do l = kb - k, 1, -1 nrt = ( n-j2+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l & ), inca,work( n+j2-m ), work( j2-m ), ka1 ) end do end do if( kb>1_${ik}$ ) then do j = n - 1, i - kb + 2*ka + 1, -1 work( n+j-m ) = work( n+j-ka-m ) work( j-m ) = work( j-ka-m ) end do end if else ! transform a, working with the lower triangle if( update ) then ! form inv(s(i))**t * a * inv(s(i)) bii = bb( 1_${ik}$, i ) do j = i, i1 ab( j-i+1, i ) = ab( j-i+1, i ) / bii end do do j = max( 1, i-ka ), i ab( i-j+1, j ) = ab( i-j+1, j ) / bii end do do k = i - kbt, i - 1 do j = i - kbt, k ab( k-j+1, j ) = ab( k-j+1, j ) -bb( i-j+1, j )*ab( i-k+1, k ) -bb( i-k+1, & k )*ab( i-j+1, j ) +ab( 1_${ik}$, i )*bb( i-j+1, j )*bb( i-k+1, k ) end do do j = max( 1, i-ka ), i - kbt - 1 ab( k-j+1, j ) = ab( k-j+1, j ) -bb( i-k+1, k )*ab( i-j+1, j ) end do end do do j = i, i1 do k = max( j-ka, i-kbt ), i - 1 ab( j-k+1, k ) = ab( j-k+1, k ) -bb( i-k+1, k )*ab( j-i+1, i ) end do end do if( wantx ) then ! post-multiply x by inv(s(i)) call stdlib${ii}$_sscal( n-m, one / bii, x( m+1, i ), 1_${ik}$ ) if( kbt>0_${ik}$ )call stdlib${ii}$_sger( n-m, kbt, -one, x( m+1, i ), 1_${ik}$,bb( kbt+1, i-kbt )& , ldbb-1,x( m+1, i-kbt ), ldx ) end if ! store a(i1,i) in ra1 for use in next loop over k ra1 = ab( i1-i+1, i ) end if ! generate and apply vectors of rotations to chase all the ! existing bulges ka positions down toward the bottom of the ! band loop_360: do k = 1, kb - 1 if( update ) then ! determine the rotations which would annihilate the bulge ! which has in theory just been created if( i-k+ka<n .and. i-k>1_${ik}$ ) then ! generate rotation to annihilate a(i-k+ka+1,i) call stdlib${ii}$_slartg( ab( ka1-k, i ), ra1, work( n+i-k+ka-m ),work( i-k+ka-m & ), ra ) ! create nonzero element a(i-k+ka+1,i-k) outside the ! band and store it in work(i-k) t = -bb( k+1, i-k )*ra1 work( i-k ) = work( n+i-k+ka-m )*t -work( i-k+ka-m )*ab( ka1, i-k ) ab( ka1, i-k ) = work( i-k+ka-m )*t +work( n+i-k+ka-m )*ab( ka1, i-k ) ra1 = ra end if end if j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1 nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 if( update ) then j2t = max( j2, i+2*ka-k+1 ) else j2t = j2 end if nrt = ( n-j2t+ka ) / ka1 do j = j2t, j1, ka1 ! create nonzero element a(j+1,j-ka) outside the band ! and store it in work(j-m) work( j-m ) = work( j-m )*ab( ka1, j-ka+1 ) ab( ka1, j-ka+1 ) = work( n+j-m )*ab( ka1, j-ka+1 ) end do ! generate rotations in 1st set to annihilate elements which ! have been created outside the band if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( ka1, j2t-ka ), inca, work( j2t-m ),ka1, & work( n+j2t-m ), ka1 ) if( nr>0_${ik}$ ) then ! apply rotations in 1st set from the left do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( l+1, j2-l ), inca,ab( l+2, j2-l ), inca, work( & n+j2-m ),work( j2-m ), ka1 ) end do ! apply rotations in 1st set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j2 ), ab( 1_${ik}$, j2+1 ), ab( 2_${ik}$, j2 ),inca, work( n+& j2-m ), work( j2-m ), ka1 ) end if ! start applying rotations in 1st set from the right do l = ka - 1, kb - k + 1, -1 nrt = ( n-j2+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2 ), inca,ab( ka1-l, j2+1 ),& inca, work( n+j2-m ),work( j2-m ), ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 1st set do j = j2, j1, ka1 call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j-m ), & work( j-m ) ) end do end if end do loop_360 if( update ) then if( i2<=n .and. kbt>0_${ik}$ ) then ! create nonzero element a(i-kbt+ka+1,i-kbt) outside the ! band and store it in work(i-kbt) work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1 end if end if loop_400: do k = kb, 1, -1 if( update ) then j2 = i - k - 1_${ik}$ + max( 2_${ik}$, k-i0+1 )*ka1 else j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1 end if ! finish applying rotations in 2nd set from the right do l = kb - k, 1, -1 nrt = ( n-j2+ka+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2-ka ), inca,ab( ka1-l, j2-& ka+1 ), inca,work( n+j2-ka ), work( j2-ka ), ka1 ) end do nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 do j = j1, j2, -ka1 work( j ) = work( j-ka ) work( n+j ) = work( n+j-ka ) end do do j = j2, j1, ka1 ! create nonzero element a(j+1,j-ka) outside the band ! and store it in work(j) work( j ) = work( j )*ab( ka1, j-ka+1 ) ab( ka1, j-ka+1 ) = work( n+j )*ab( ka1, j-ka+1 ) end do if( update ) then if( i-k<n-ka .and. k<=kbt )work( i-k+ka ) = work( i-k ) end if end do loop_400 loop_440: do k = kb, 1, -1 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1 nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 if( nr>0_${ik}$ ) then ! generate rotations in 2nd set to annihilate elements ! which have been created outside the band call stdlib${ii}$_slargv( nr, ab( ka1, j2-ka ), inca, work( j2 ), ka1,work( n+j2 ), & ka1 ) ! apply rotations in 2nd set from the left do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( l+1, j2-l ), inca,ab( l+2, j2-l ), inca, work( & n+j2 ),work( j2 ), ka1 ) end do ! apply rotations in 2nd set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j2 ), ab( 1_${ik}$, j2+1 ), ab( 2_${ik}$, j2 ),inca, work( n+& j2 ), work( j2 ), ka1 ) end if ! start applying rotations in 2nd set from the right do l = ka - 1, kb - k + 1, -1 nrt = ( n-j2+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2 ), inca,ab( ka1-l, j2+1 ),& inca, work( n+j2 ),work( j2 ), ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 2nd set do j = j2, j1, ka1 call stdlib${ii}$_srot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j ), work( & j ) ) end do end if end do loop_440 do k = 1, kb - 1 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1 ! finish applying rotations in 1st set from the right do l = kb - k, 1, -1 nrt = ( n-j2+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j2 ), inca,ab( ka1-l, j2+1 ),& inca, work( n+j2-m ),work( j2-m ), ka1 ) end do end do if( kb>1_${ik}$ ) then do j = n - 1, i - kb + 2*ka + 1, -1 work( n+j-m ) = work( n+j-ka-m ) work( j-m ) = work( j-ka-m ) end do end if end if go to 10 480 continue ! **************************** phase 2 ***************************** ! the logical structure of this phase is: ! update = .true. ! do i = 1, m ! use s(i) to update a and create a new bulge ! apply rotations to push all bulges ka positions upward ! end do ! update = .false. ! do i = m - ka - 1, 2, -1 ! apply rotations to push all bulges ka positions upward ! end do ! to avoid duplicating code, the two loops are merged. update = .true. i = 0_${ik}$ 490 continue if( update ) then i = i + 1_${ik}$ kbt = min( kb, m-i ) i0 = i + 1_${ik}$ i1 = max( 1_${ik}$, i-ka ) i2 = i + kbt - ka1 if( i>m ) then update = .false. i = i - 1_${ik}$ i0 = m + 1_${ik}$ if( ka==0 )return go to 490 end if else i = i - ka if( i<2 )return end if if( i<m-kbt ) then nx = m else nx = n end if if( upper ) then ! transform a, working with the upper triangle if( update ) then ! form inv(s(i))**t * a * inv(s(i)) bii = bb( kb1, i ) do j = i1, i ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii end do do j = i, min( n, i+ka ) ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii end do do k = i + 1, i + kbt do j = k, i + kbt ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( i-j+kb1, j )*ab( i-k+ka1, k ) -bb(& i-k+kb1, k )*ab( i-j+ka1, j ) +ab( ka1, i )*bb( i-j+kb1, j )*bb( i-k+kb1, & k ) end do do j = i + kbt + 1, min( n, i+ka ) ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( i-k+kb1, k )*ab( i-j+ka1, j ) end do end do do j = i1, i do k = i + 1, min( j+ka, i+kbt ) ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( i-k+kb1, k )*ab( j-i+ka1, i ) end do end do if( wantx ) then ! post-multiply x by inv(s(i)) call stdlib${ii}$_sscal( nx, one / bii, x( 1_${ik}$, i ), 1_${ik}$ ) if( kbt>0_${ik}$ )call stdlib${ii}$_sger( nx, kbt, -one, x( 1_${ik}$, i ), 1_${ik}$, bb( kb, i+1 ),ldbb-& 1_${ik}$, x( 1_${ik}$, i+1 ), ldx ) end if ! store a(i1,i) in ra1 for use in next loop over k ra1 = ab( i1-i+ka1, i ) end if ! generate and apply vectors of rotations to chase all the ! existing bulges ka positions up toward the top of the band loop_610: do k = 1, kb - 1 if( update ) then ! determine the rotations which would annihilate the bulge ! which has in theory just been created if( i+k-ka1>0_${ik}$ .and. i+k<m ) then ! generate rotation to annihilate a(i+k-ka-1,i) call stdlib${ii}$_slartg( ab( k+1, i ), ra1, work( n+i+k-ka ),work( i+k-ka ), ra & ) ! create nonzero element a(i+k-ka-1,i+k) outside the ! band and store it in work(m-kb+i+k) t = -bb( kb1-k, i+k )*ra1 work( m-kb+i+k ) = work( n+i+k-ka )*t -work( i+k-ka )*ab( 1_${ik}$, i+k ) ab( 1_${ik}$, i+k ) = work( i+k-ka )*t +work( n+i+k-ka )*ab( 1_${ik}$, i+k ) ra1 = ra end if end if j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1 nr = ( j2+ka-1 ) / ka1 j1 = j2 - ( nr-1 )*ka1 if( update ) then j2t = min( j2, i-2*ka+k-1 ) else j2t = j2 end if nrt = ( j2t+ka-1 ) / ka1 do j = j1, j2t, ka1 ! create nonzero element a(j-1,j+ka) outside the band ! and store it in work(j) work( j ) = work( j )*ab( 1_${ik}$, j+ka-1 ) ab( 1_${ik}$, j+ka-1 ) = work( n+j )*ab( 1_${ik}$, j+ka-1 ) end do ! generate rotations in 1st set to annihilate elements which ! have been created outside the band if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( 1_${ik}$, j1+ka ), inca, work( j1 ), ka1,work( & n+j1 ), ka1 ) if( nr>0_${ik}$ ) then ! apply rotations in 1st set from the left do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( ka1-l, j1+l ), inca,ab( ka-l, j1+l ), inca, & work( n+j1 ),work( j1 ), ka1 ) end do ! apply rotations in 1st set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),ab( ka, j1 ), inca, & work( n+j1 ),work( j1 ), ka1 ) end if ! start applying rotations in 1st set from the right do l = ka - 1, kb - k + 1, -1 nrt = ( j2+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t ), inca,ab( l+1, j1t-1 ), inca,& work( n+j1t ),work( j1t ), ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 1st set do j = j1, j2, ka1 call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+j ), work( j ) ) end do end if end do loop_610 if( update ) then if( i2>0_${ik}$ .and. kbt>0_${ik}$ ) then ! create nonzero element a(i+kbt-ka-1,i+kbt) outside the ! band and store it in work(m-kb+i+kbt) work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1 end if end if loop_650: do k = kb, 1, -1 if( update ) then j2 = i + k + 1_${ik}$ - max( 2_${ik}$, k+i0-m )*ka1 else j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1 end if ! finish applying rotations in 2nd set from the right do l = kb - k, 1, -1 nrt = ( j2+ka+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t+ka ), inca,ab( l+1, j1t+ka-1 ),& inca,work( n+m-kb+j1t+ka ),work( m-kb+j1t+ka ), ka1 ) end do nr = ( j2+ka-1 ) / ka1 j1 = j2 - ( nr-1 )*ka1 do j = j1, j2, ka1 work( m-kb+j ) = work( m-kb+j+ka ) work( n+m-kb+j ) = work( n+m-kb+j+ka ) end do do j = j1, j2, ka1 ! create nonzero element a(j-1,j+ka) outside the band ! and store it in work(m-kb+j) work( m-kb+j ) = work( m-kb+j )*ab( 1_${ik}$, j+ka-1 ) ab( 1_${ik}$, j+ka-1 ) = work( n+m-kb+j )*ab( 1_${ik}$, j+ka-1 ) end do if( update ) then if( i+k>ka1 .and. k<=kbt )work( m-kb+i+k-ka ) = work( m-kb+i+k ) end if end do loop_650 loop_690: do k = kb, 1, -1 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1 nr = ( j2+ka-1 ) / ka1 j1 = j2 - ( nr-1 )*ka1 if( nr>0_${ik}$ ) then ! generate rotations in 2nd set to annihilate elements ! which have been created outside the band call stdlib${ii}$_slargv( nr, ab( 1_${ik}$, j1+ka ), inca, work( m-kb+j1 ),ka1, work( n+m-& kb+j1 ), ka1 ) ! apply rotations in 2nd set from the left do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( ka1-l, j1+l ), inca,ab( ka-l, j1+l ), inca,& work( n+m-kb+j1 ), work( m-kb+j1 ), ka1 ) end do ! apply rotations in 2nd set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),ab( ka, j1 ), inca, & work( n+m-kb+j1 ),work( m-kb+j1 ), ka1 ) end if ! start applying rotations in 2nd set from the right do l = ka - 1, kb - k + 1, -1 nrt = ( j2+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t ), inca,ab( l+1, j1t-1 ), inca,& work( n+m-kb+j1t ), work( m-kb+j1t ),ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 2nd set do j = j1, j2, ka1 call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+m-kb+j ), work( & m-kb+j ) ) end do end if end do loop_690 do k = 1, kb - 1 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1 ! finish applying rotations in 1st set from the right do l = kb - k, 1, -1 nrt = ( j2+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l, j1t ), inca,ab( l+1, j1t-1 ), inca,& work( n+j1t ),work( j1t ), ka1 ) end do end do if( kb>1_${ik}$ ) then do j = 2, min( i+kb, m ) - 2*ka - 1 work( n+j ) = work( n+j+ka ) work( j ) = work( j+ka ) end do end if else ! transform a, working with the lower triangle if( update ) then ! form inv(s(i))**t * a * inv(s(i)) bii = bb( 1_${ik}$, i ) do j = i1, i ab( i-j+1, j ) = ab( i-j+1, j ) / bii end do do j = i, min( n, i+ka ) ab( j-i+1, i ) = ab( j-i+1, i ) / bii end do do k = i + 1, i + kbt do j = k, i + kbt ab( j-k+1, k ) = ab( j-k+1, k ) -bb( j-i+1, i )*ab( k-i+1, i ) -bb( k-i+1, & i )*ab( j-i+1, i ) +ab( 1_${ik}$, i )*bb( j-i+1, i )*bb( k-i+1, i ) end do do j = i + kbt + 1, min( n, i+ka ) ab( j-k+1, k ) = ab( j-k+1, k ) -bb( k-i+1, i )*ab( j-i+1, i ) end do end do do j = i1, i do k = i + 1, min( j+ka, i+kbt ) ab( k-j+1, j ) = ab( k-j+1, j ) -bb( k-i+1, i )*ab( i-j+1, j ) end do end do if( wantx ) then ! post-multiply x by inv(s(i)) call stdlib${ii}$_sscal( nx, one / bii, x( 1_${ik}$, i ), 1_${ik}$ ) if( kbt>0_${ik}$ )call stdlib${ii}$_sger( nx, kbt, -one, x( 1_${ik}$, i ), 1_${ik}$, bb( 2_${ik}$, i ), 1_${ik}$,x( 1_${ik}$, & i+1 ), ldx ) end if ! store a(i,i1) in ra1 for use in next loop over k ra1 = ab( i-i1+1, i1 ) end if ! generate and apply vectors of rotations to chase all the ! existing bulges ka positions up toward the top of the band loop_840: do k = 1, kb - 1 if( update ) then ! determine the rotations which would annihilate the bulge ! which has in theory just been created if( i+k-ka1>0_${ik}$ .and. i+k<m ) then ! generate rotation to annihilate a(i,i+k-ka-1) call stdlib${ii}$_slartg( ab( ka1-k, i+k-ka ), ra1,work( n+i+k-ka ), work( i+k-& ka ), ra ) ! create nonzero element a(i+k,i+k-ka-1) outside the ! band and store it in work(m-kb+i+k) t = -bb( k+1, i )*ra1 work( m-kb+i+k ) = work( n+i+k-ka )*t -work( i+k-ka )*ab( ka1, i+k-ka ) ab( ka1, i+k-ka ) = work( i+k-ka )*t +work( n+i+k-ka )*ab( ka1, i+k-ka ) ra1 = ra end if end if j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1 nr = ( j2+ka-1 ) / ka1 j1 = j2 - ( nr-1 )*ka1 if( update ) then j2t = min( j2, i-2*ka+k-1 ) else j2t = j2 end if nrt = ( j2t+ka-1 ) / ka1 do j = j1, j2t, ka1 ! create nonzero element a(j+ka,j-1) outside the band ! and store it in work(j) work( j ) = work( j )*ab( ka1, j-1 ) ab( ka1, j-1 ) = work( n+j )*ab( ka1, j-1 ) end do ! generate rotations in 1st set to annihilate elements which ! have been created outside the band if( nrt>0_${ik}$ )call stdlib${ii}$_slargv( nrt, ab( ka1, j1 ), inca, work( j1 ), ka1,work( n+& j1 ), ka1 ) if( nr>0_${ik}$ ) then ! apply rotations in 1st set from the right do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),inca, work( n+& j1 ), work( j1 ), ka1 ) end do ! apply rotations in 1st set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j1 ), ab( 1_${ik}$, j1-1 ),ab( 2_${ik}$, j1-1 ), inca, work( & n+j1 ),work( j1 ), ka1 ) end if ! start applying rotations in 1st set from the left do l = ka - 1, kb - k + 1, -1 nrt = ( j2+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,ab( ka1-l, & j1t-ka1+l ), inca,work( n+j1t ), work( j1t ), ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 1st set do j = j1, j2, ka1 call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+j ), work( j ) ) end do end if end do loop_840 if( update ) then if( i2>0_${ik}$ .and. kbt>0_${ik}$ ) then ! create nonzero element a(i+kbt,i+kbt-ka-1) outside the ! band and store it in work(m-kb+i+kbt) work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1 end if end if loop_880: do k = kb, 1, -1 if( update ) then j2 = i + k + 1_${ik}$ - max( 2_${ik}$, k+i0-m )*ka1 else j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1 end if ! finish applying rotations in 2nd set from the left do l = kb - k, 1, -1 nrt = ( j2+ka+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,ab( ka1-l, & j1t+l-1 ), inca,work( n+m-kb+j1t+ka ),work( m-kb+j1t+ka ), ka1 ) end do nr = ( j2+ka-1 ) / ka1 j1 = j2 - ( nr-1 )*ka1 do j = j1, j2, ka1 work( m-kb+j ) = work( m-kb+j+ka ) work( n+m-kb+j ) = work( n+m-kb+j+ka ) end do do j = j1, j2, ka1 ! create nonzero element a(j+ka,j-1) outside the band ! and store it in work(m-kb+j) work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 ) ab( ka1, j-1 ) = work( n+m-kb+j )*ab( ka1, j-1 ) end do if( update ) then if( i+k>ka1 .and. k<=kbt )work( m-kb+i+k-ka ) = work( m-kb+i+k ) end if end do loop_880 loop_920: do k = kb, 1, -1 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m )*ka1 nr = ( j2+ka-1 ) / ka1 j1 = j2 - ( nr-1 )*ka1 if( nr>0_${ik}$ ) then ! generate rotations in 2nd set to annihilate elements ! which have been created outside the band call stdlib${ii}$_slargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),ka1, work( n+m-& kb+j1 ), ka1 ) ! apply rotations in 2nd set from the right do l = 1, ka - 1 call stdlib${ii}$_slartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),inca, work( n+& m-kb+j1 ), work( m-kb+j1 ),ka1 ) end do ! apply rotations in 2nd set from both sides to diagonal ! blocks call stdlib${ii}$_slar2v( nr, ab( 1_${ik}$, j1 ), ab( 1_${ik}$, j1-1 ),ab( 2_${ik}$, j1-1 ), inca, work( & n+m-kb+j1 ),work( m-kb+j1 ), ka1 ) end if ! start applying rotations in 2nd set from the left do l = ka - 1, kb - k + 1, -1 nrt = ( j2+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,ab( ka1-l, & j1t-ka1+l ), inca,work( n+m-kb+j1t ), work( m-kb+j1t ),ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 2nd set do j = j1, j2, ka1 call stdlib${ii}$_srot( nx, x( 1_${ik}$, j ), 1_${ik}$, x( 1_${ik}$, j-1 ), 1_${ik}$,work( n+m-kb+j ), work( & m-kb+j ) ) end do end if end do loop_920 do k = 1, kb - 1 j2 = i + k + 1_${ik}$ - max( 1_${ik}$, k+i0-m+1 )*ka1 ! finish applying rotations in 1st set from the left do l = kb - k, 1, -1 nrt = ( j2+l-1 ) / ka1 j1t = j2 - ( nrt-1 )*ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,ab( ka1-l, & j1t-ka1+l ), inca,work( n+j1t ), work( j1t ), ka1 ) end do end do if( kb>1_${ik}$ ) then do j = 2, min( i+kb, m ) - 2*ka - 1 work( n+j ) = work( n+j+ka ) work( j ) = work( j+ka ) end do end if end if go to 490 end subroutine stdlib${ii}$_ssbgst pure module subroutine stdlib${ii}$_dsbgst( vect, uplo, n, ka, kb, ab, ldab, bb, ldbb, x,ldx, work, info ) !! DSBGST reduces a real symmetric-definite banded generalized !! eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y, !! such that C has the same bandwidth as A. !! B must have been previously factorized as S**T*S by DPBSTF, using a !! split Cholesky factorization. A is overwritten by C = X**T*A*X, where !! X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the !! bandwidth of A. ! -- 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) :: uplo, vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ka, kb, ldab, ldbb, ldx, n ! Array Arguments real(dp), intent(inout) :: ab(ldab,*) real(dp), intent(in) :: bb(ldbb,*) real(dp), intent(out) :: work(*), x(ldx,*) ! ===================================================================== ! Local Scalars logical(lk) :: update, upper, wantx integer(${ik}$) :: i, i0, i1, i2, inca, j, j1, j1t, j2, j2t, k, ka1, kb1, kbt, l, m, nr, & nrt, nx real(dp) :: bii, ra, ra1, t ! Intrinsic Functions ! Executable Statements ! test the input parameters wantx = stdlib_lsame( vect, 'V' ) upper = stdlib_lsame( uplo, 'U' ) ka1 = ka + 1_${ik}$ kb1 = kb + 1_${ik}$ info = 0_${ik}$ if( .not.wantx .and. .not.stdlib_lsame( vect, 'N' ) ) then info = -1_${ik}$ else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ka<0_${ik}$ ) then info = -4_${ik}$ else if( kb<0_${ik}$ .or. kb>ka ) then info = -5_${ik}$ else if( ldab<ka+1 ) then info = -7_${ik}$ else if( ldbb<kb+1 ) then info = -9_${ik}$ else if( ldx<1_${ik}$ .or. wantx .and. ldx<max( 1_${ik}$, n ) ) then info = -11_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DSBGST', -info ) return end if ! quick return if possible if( n==0 )return inca = ldab*ka1 ! initialize x to the unit matrix, if needed if( wantx )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, x, ldx ) ! set m to the splitting point m. it must be the same value as is ! used in stdlib${ii}$_dpbstf. the chosen value allows the arrays work and rwork ! to be of dimension (n). m = ( n+kb ) / 2_${ik}$ ! the routine works in two phases, corresponding to the two halves ! of the split cholesky factorization of b as s**t*s where ! s = ( u ) ! ( m l ) ! with u upper triangular of order m, and l lower triangular of ! order n-m. s has the same bandwidth as b. ! s is treated as a product of elementary matrices: ! s = s(m)*s(m-1)*...*s(2)*s(1)*s(m+1)*s(m+2)*...*s(n-1)*s(n) ! where s(i) is determined by the i-th row of s. ! in phase 1, the index i takes the values n, n-1, ... , m+1; ! in phase 2, it takes the values 1, 2, ... , m. ! for each value of i, the current matrix a is updated by forming ! inv(s(i))**t*a*inv(s(i)). this creates a triangular bulge outside ! the band of a. the bulge is then pushed down toward the bottom of ! a in phase 1, and up toward the top of a in phase 2, by applying ! plane rotations. ! there are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1 ! of them are linearly independent, so annihilating a bulge requires ! only 2*kb-1 plane rotations. the rotations are divided into a 1st ! set of kb-1 rotations, and a 2nd set of kb rotations. ! wherever possible, rotations are generated and applied in vector ! operations of length nr between the indices j1 and j2 (sometimes ! replaced by modified values nrt, j1t or j2t). ! the cosines and sines of the rotations are stored in the array ! work. the cosines of the 1st set of rotations are stored in ! elements n+2:n+m-kb-1 and the sines of the 1st set in elements ! 2:m-kb-1; the cosines of the 2nd set are stored in elements ! n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n. ! the bulges are not formed explicitly; nonzero elements outside the ! band are created only when they are required for generating new ! rotations; they are stored in the array work, in positions where ! they are later overwritten by the sines of the rotations which ! annihilate them. ! **************************** phase 1 ***************************** ! the logical structure of this phase is: ! update = .true. ! do i = n, m + 1, -1 ! use s(i) to update a and create a new bulge ! apply rotations to push all bulges ka positions downward ! end do ! update = .false. ! do i = m + ka + 1, n - 1 ! apply rotations to push all bulges ka positions downward ! end do ! to avoid duplicating code, the two loops are merged. update = .true. i = n + 1_${ik}$ 10 continue if( update ) then i = i - 1_${ik}$ kbt = min( kb, i-1 ) i0 = i - 1_${ik}$ i1 = min( n, i+ka ) i2 = i - kbt + ka1 if( i<m+1 ) then update = .false. i = i + 1_${ik}$ i0 = m if( ka==0 )go to 480 go to 10 end if else i = i + ka if( i>n-1 )go to 480 end if if( upper ) then ! transform a, working with the upper triangle if( update ) then ! form inv(s(i))**t * a * inv(s(i)) bii = bb( kb1, i ) do j = i, i1 ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii end do do j = max( 1, i-ka ), i ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii end do do k = i - kbt, i - 1 do j = i - kbt, k ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( j-i+kb1, i )*ab( k-i+ka1, i ) -bb(& k-i+kb1, i )*ab( j-i+ka1, i ) +ab( ka1, i )*bb( j-i+kb1, i )*bb( k-i+kb1, & i ) end do do j = max( 1, i-ka ), i - kbt - 1 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -bb( k-i+kb1, i )*ab( j-i+ka1, i ) end do end do do j = i, i1 do k = max( j-ka, i-kbt ), i - 1 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -bb( k-i+kb1, i )*ab( i-j+ka1, j ) end do end do if( wantx ) then ! post-multiply x by inv(s(i)) call stdlib${ii}$_dscal( n-m, one / bii, x( m+1, i ), 1_${ik}$ ) if( kbt>0_${ik}$ )call stdlib${ii}$_dger( n-m, kbt, -one, x( m+1, i ), 1_${ik}$,bb( kb1-kbt, i ), & 1_${ik}$, x( m+1, i-kbt ), ldx ) end if ! store a(i,i1) in ra1 for use in next loop over k ra1 = ab( i-i1+ka1, i1 ) end if ! generate and apply vectors of rotations to chase all the ! existing bulges ka positions down toward the bottom of the ! band loop_130: do k = 1, kb - 1 if( update ) then ! determine the rotations which would annihilate the bulge ! which has in theory just been created if( i-k+ka<n .and. i-k>1_${ik}$ ) then ! generate rotation to annihilate a(i,i-k+ka+1) call stdlib${ii}$_dlartg( ab( k+1, i-k+ka ), ra1,work( n+i-k+ka-m ), work( i-k+& ka-m ),ra ) ! create nonzero element a(i-k,i-k+ka+1) outside the ! band and store it in work(i-k) t = -bb( kb1-k, i )*ra1 work( i-k ) = work( n+i-k+ka-m )*t -work( i-k+ka-m )*ab( 1_${ik}$, i-k+ka ) ab( 1_${ik}$, i-k+ka ) = work( i-k+ka-m )*t +work( n+i-k+ka-m )*ab( 1_${ik}$, i-k+ka ) ra1 = ra end if end if j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+2 )*ka1 nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 if( update ) then j2t = max( j2, i+2*ka-k+1 ) else j2t = j2 end if nrt = ( n-j2t+ka ) / ka1 do j = j2t, j1, ka1 ! create nonzero element a(j-ka,j+1) outside the band ! and store it in work(j-m) work( j-m ) = work( j-m )*ab( 1_${ik}$, j+1 ) ab( 1_${ik}$, j+1 ) = work( n+j-m )*ab( 1_${ik}$, j+1 ) end do ! generate rotations in 1st set to annihilate elements which ! have been created outside the band if( nrt>0_${ik}$ )call stdlib${ii}$_dlargv( nrt, ab( 1_${ik}$, j2t ), inca, work( j2t-m ), ka1,work( & n+j2t-m ), ka1 ) if( nr>0_${ik}$ ) then ! apply rotations in 1st set from the right do l = 1, ka - 1 call stdlib${ii}$_dlartv( nr, ab( ka1-l, j2 ), inca,ab( ka-l, j2+1 ), inca, work(& n+j2-m ),work( j2-m ), ka1 ) end do ! apply rotations in 1st set from both sides to diagonal ! blocks call stdlib${ii}$_dlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),ab( ka, j2+1 ), inca, & work( n+j2-m ),work( j2-m ), ka1 ) end if ! start applying rotations in 1st set from the left do l = ka - 1, kb - k + 1, -1 nrt = ( n-j2+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( l, j2+ka1-l ), inca,ab( l+1, j2+ka1-l & ), inca,work( n+j2-m ), work( j2-m ), ka1 ) end do if( wantx ) then ! post-multiply x by product of rotations in 1st set do j = j2, j1, ka1 call stdlib${ii}$_drot( n-m, x( m+1, j ), 1_${ik}$, x( m+1, j+1 ), 1_${ik}$,work( n+j-m ), & work( j-m ) ) end do end if end do loop_130 if( update ) then if( i2<=n .and. kbt>0_${ik}$ ) then ! create nonzero element a(i-kbt,i-kbt+ka+1) outside the ! band and store it in work(i-kbt) work( i-kbt ) = -bb( kb1-kbt, i )*ra1 end if end if loop_170: do k = kb, 1, -1 if( update ) then j2 = i - k - 1_${ik}$ + max( 2_${ik}$, k-i0+1 )*ka1 else j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1 end if ! finish applying rotations in 2nd set from the left do l = kb - k, 1, -1 nrt = ( n-j2+ka+l ) / ka1 if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( l, j2-l+1 ), inca,ab( l+1, j2-l+1 ), & inca, work( n+j2-ka ),work( j2-ka ), ka1 ) end do nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 do j = j1, j2, -ka1 work( j ) = work( j-ka ) work( n+j ) = work( n+j-ka ) end do do j = j2, j1, ka1 ! create nonzero element a(j-ka,j+1) outside the band ! and store it in work(j) work( j ) = work( j )*ab( 1_${ik}$, j+1 ) ab( 1_${ik}$, j+1 ) = work( n+j )*ab( 1_${ik}$, j+1 ) end do if( update ) then if( i-k<n-ka .and. k<=kbt )work( i-k+ka ) = work( i-k ) end if end do loop_170 loop_210: do k = kb, 1, -1 j2 = i - k - 1_${ik}$ + max( 1_${ik}$, k-i0+1 )*ka1 nr = ( n-j2+ka ) / ka1 j1 = j2 + ( nr-1 )*ka1 if( nr>0_${ik}$ ) then ! generate rotations in 2nd set to annihilate elements ! which have been created outside the band call stdlib${ii}$_dlargv( nr, ab( 1_${ik}$, j2