#:include "common.fypp" submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_eigv_tridiag implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_slaebz( ijob, nitmax, n, mmax, minp, nbmin, abstol,reltol, pivmin, d, & !! SLAEBZ contains the iteration loops which compute and use the !! function N(w), which is the count of eigenvalues of a symmetric !! tridiagonal matrix T less than or equal to its argument w. It !! performs a choice of two types of loops: !! IJOB=1, followed by !! IJOB=2: It takes as input a list of intervals and returns a list of !! sufficiently small intervals whose union contains the same !! eigenvalues as the union of the original intervals. !! The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. !! The output interval (AB(j,1),AB(j,2)] will contain !! eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. !! IJOB=3: It performs a binary search in each input interval !! (AB(j,1),AB(j,2)] for a point w(j) such that !! N(w(j))=NVAL(j), and uses C(j) as the starting point of !! the search. If such a w(j) is found, then on output !! AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output !! (AB(j,1),AB(j,2)] will be a small interval containing the !! point where N(w) jumps through NVAL(j), unless that point !! lies outside the initial interval. !! Note that the intervals are in all cases half-open intervals, !! i.e., of the form (a,b] , which includes b but not a . !! To avoid underflow, the matrix should be scaled so that its largest !! element is no greater than overflow**(1/2) * underflow**(1/4) !! in absolute value. To assure the most accurate computation !! of small eigenvalues, the matrix should be scaled to be !! not much smaller than that, either. !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal !! Matrix", Report CS41, Computer Science Dept., Stanford !! University, July 21, 1966 !! Note: the arguments are, in general, *not* checked for unreasonable !! values. e, e2, nval, ab, c, mout,nab, work, iwork, info ) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: ijob, minp, mmax, n, nbmin, nitmax integer(${ik}$), intent(out) :: info, mout real(sp), intent(in) :: abstol, pivmin, reltol ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) integer(${ik}$), intent(inout) :: nab(mmax,*), nval(*) real(sp), intent(inout) :: ab(mmax,*), c(*) real(sp), intent(in) :: d(*), e(*), e2(*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl, klnew real(sp) :: tmp1, tmp2 ! Intrinsic Functions ! Executable Statements ! check for errors info = 0_${ik}$ if( ijob<1_${ik}$ .or. ijob>3_${ik}$ ) then info = -1_${ik}$ return end if ! initialize nab if( ijob==1_${ik}$ ) then ! compute the number of eigenvalues in the initial intervals. mout = 0_${ik}$ do ji = 1, minp do jp = 1, 2 tmp1 = d( 1_${ik}$ ) - ab( ji, jp ) if( abs( tmp1 )<pivmin )tmp1 = -pivmin nab( ji, jp ) = 0_${ik}$ if( tmp1<=zero )nab( ji, jp ) = 1_${ik}$ do j = 2, n tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp ) if( abs( tmp1 )<pivmin )tmp1 = -pivmin if( tmp1<=zero )nab( ji, jp ) = nab( ji, jp ) + 1_${ik}$ end do end do mout = mout + nab( ji, 2_${ik}$ ) - nab( ji, 1_${ik}$ ) end do return end if ! initialize for loop ! kf and kl have the following meaning: ! intervals 1,...,kf-1 have converged. ! intervals kf,...,kl still need to be refined. kf = 1_${ik}$ kl = minp ! if ijob=2, initialize c. ! if ijob=3, use the user-supplied starting point. if( ijob==2_${ik}$ ) then do ji = 1, minp c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) ) end do end if ! iteration loop loop_130: do jit = 1, nitmax ! loop over intervals if( kl-kf+1>=nbmin .and. nbmin>0_${ik}$ ) then ! begin of parallel version of the loop do ji = kf, kl ! compute n(c), the number of eigenvalues less than c work( ji ) = d( 1_${ik}$ ) - c( ji ) iwork( ji ) = 0_${ik}$ if( work( ji )<=pivmin ) then iwork( ji ) = 1_${ik}$ work( ji ) = min( work( ji ), -pivmin ) end if do j = 2, n work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji ) if( work( ji )<=pivmin ) then iwork( ji ) = iwork( ji ) + 1_${ik}$ work( ji ) = min( work( ji ), -pivmin ) end if end do end do if( ijob<=2_${ik}$ ) then ! ijob=2: choose all intervals containing eigenvalues. klnew = kl loop_70: do ji = kf, kl ! insure that n(w) is monotone iwork( ji ) = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), iwork( ji ) ) ) ! update the queue -- add intervals if both halves ! contain eigenvalues. if( iwork( ji )==nab( ji, 2_${ik}$ ) ) then ! no eigenvalue in the upper interval: ! just use the lower interval. ab( ji, 2_${ik}$ ) = c( ji ) else if( iwork( ji )==nab( ji, 1_${ik}$ ) ) then ! no eigenvalue in the lower interval: ! just use the upper interval. ab( ji, 1_${ik}$ ) = c( ji ) else klnew = klnew + 1_${ik}$ if( klnew<=mmax ) then ! eigenvalue in both intervals -- add upper to ! queue. ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ ) nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ ) ab( klnew, 1_${ik}$ ) = c( ji ) nab( klnew, 1_${ik}$ ) = iwork( ji ) ab( ji, 2_${ik}$ ) = c( ji ) nab( ji, 2_${ik}$ ) = iwork( ji ) else info = mmax + 1_${ik}$ end if end if end do loop_70 if( info/=0 )return kl = klnew else ! ijob=3: binary search. keep only the interval containing ! w s.t. n(w) = nval do ji = kf, kl if( iwork( ji )<=nval( ji ) ) then ab( ji, 1_${ik}$ ) = c( ji ) nab( ji, 1_${ik}$ ) = iwork( ji ) end if if( iwork( ji )>=nval( ji ) ) then ab( ji, 2_${ik}$ ) = c( ji ) nab( ji, 2_${ik}$ ) = iwork( ji ) end if end do end if else ! end of parallel version of the loop ! begin of serial version of the loop klnew = kl loop_100: do ji = kf, kl ! compute n(w), the number of eigenvalues less than w tmp1 = c( ji ) tmp2 = d( 1_${ik}$ ) - tmp1 itmp1 = 0_${ik}$ if( tmp2<=pivmin ) then itmp1 = 1_${ik}$ tmp2 = min( tmp2, -pivmin ) end if do j = 2, n tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1 if( tmp2<=pivmin ) then itmp1 = itmp1 + 1_${ik}$ tmp2 = min( tmp2, -pivmin ) end if end do if( ijob<=2_${ik}$ ) then ! ijob=2: choose all intervals containing eigenvalues. ! insure that n(w) is monotone itmp1 = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), itmp1 ) ) ! update the queue -- add intervals if both halves ! contain eigenvalues. if( itmp1==nab( ji, 2_${ik}$ ) ) then ! no eigenvalue in the upper interval: ! just use the lower interval. ab( ji, 2_${ik}$ ) = tmp1 else if( itmp1==nab( ji, 1_${ik}$ ) ) then ! no eigenvalue in the lower interval: ! just use the upper interval. ab( ji, 1_${ik}$ ) = tmp1 else if( klnew<mmax ) then ! eigenvalue in both intervals -- add upper to queue. klnew = klnew + 1_${ik}$ ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ ) nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ ) ab( klnew, 1_${ik}$ ) = tmp1 nab( klnew, 1_${ik}$ ) = itmp1 ab( ji, 2_${ik}$ ) = tmp1 nab( ji, 2_${ik}$ ) = itmp1 else info = mmax + 1_${ik}$ return end if else ! ijob=3: binary search. keep only the interval ! containing w s.t. n(w) = nval if( itmp1<=nval( ji ) ) then ab( ji, 1_${ik}$ ) = tmp1 nab( ji, 1_${ik}$ ) = itmp1 end if if( itmp1>=nval( ji ) ) then ab( ji, 2_${ik}$ ) = tmp1 nab( ji, 2_${ik}$ ) = itmp1 end if end if end do loop_100 kl = klnew end if ! check for convergence kfnew = kf loop_110: do ji = kf, kl tmp1 = abs( ab( ji, 2_${ik}$ )-ab( ji, 1_${ik}$ ) ) tmp2 = max( abs( ab( ji, 2_${ik}$ ) ), abs( ab( ji, 1_${ik}$ ) ) ) if( tmp1<max( abstol, pivmin, reltol*tmp2 ) .or.nab( ji, 1_${ik}$ )>=nab( ji, 2_${ik}$ ) ) & then ! converged -- swap with position kfnew, ! then increment kfnew if( ji>kfnew ) then tmp1 = ab( ji, 1_${ik}$ ) tmp2 = ab( ji, 2_${ik}$ ) itmp1 = nab( ji, 1_${ik}$ ) itmp2 = nab( ji, 2_${ik}$ ) ab( ji, 1_${ik}$ ) = ab( kfnew, 1_${ik}$ ) ab( ji, 2_${ik}$ ) = ab( kfnew, 2_${ik}$ ) nab( ji, 1_${ik}$ ) = nab( kfnew, 1_${ik}$ ) nab( ji, 2_${ik}$ ) = nab( kfnew, 2_${ik}$ ) ab( kfnew, 1_${ik}$ ) = tmp1 ab( kfnew, 2_${ik}$ ) = tmp2 nab( kfnew, 1_${ik}$ ) = itmp1 nab( kfnew, 2_${ik}$ ) = itmp2 if( ijob==3_${ik}$ ) then itmp1 = nval( ji ) nval( ji ) = nval( kfnew ) nval( kfnew ) = itmp1 end if end if kfnew = kfnew + 1_${ik}$ end if end do loop_110 kf = kfnew ! choose midpoints do ji = kf, kl c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) ) end do ! if no more intervals to refine, quit. if( kf>kl )go to 140 end do loop_130 ! converged 140 continue info = max( kl+1-kf, 0_${ik}$ ) mout = kl return end subroutine stdlib${ii}$_slaebz pure module subroutine stdlib${ii}$_dlaebz( ijob, nitmax, n, mmax, minp, nbmin, abstol,reltol, pivmin, d, & !! DLAEBZ contains the iteration loops which compute and use the !! function N(w), which is the count of eigenvalues of a symmetric !! tridiagonal matrix T less than or equal to its argument w. It !! performs a choice of two types of loops: !! IJOB=1, followed by !! IJOB=2: It takes as input a list of intervals and returns a list of !! sufficiently small intervals whose union contains the same !! eigenvalues as the union of the original intervals. !! The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. !! The output interval (AB(j,1),AB(j,2)] will contain !! eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. !! IJOB=3: It performs a binary search in each input interval !! (AB(j,1),AB(j,2)] for a point w(j) such that !! N(w(j))=NVAL(j), and uses C(j) as the starting point of !! the search. If such a w(j) is found, then on output !! AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output !! (AB(j,1),AB(j,2)] will be a small interval containing the !! point where N(w) jumps through NVAL(j), unless that point !! lies outside the initial interval. !! Note that the intervals are in all cases half-open intervals, !! i.e., of the form (a,b] , which includes b but not a . !! To avoid underflow, the matrix should be scaled so that its largest !! element is no greater than overflow**(1/2) * underflow**(1/4) !! in absolute value. To assure the most accurate computation !! of small eigenvalues, the matrix should be scaled to be !! not much smaller than that, either. !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal !! Matrix", Report CS41, Computer Science Dept., Stanford !! University, July 21, 1966 !! Note: the arguments are, in general, *not* checked for unreasonable !! values. e, e2, nval, ab, c, mout,nab, work, iwork, info ) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: ijob, minp, mmax, n, nbmin, nitmax integer(${ik}$), intent(out) :: info, mout real(dp), intent(in) :: abstol, pivmin, reltol ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) integer(${ik}$), intent(inout) :: nab(mmax,*), nval(*) real(dp), intent(inout) :: ab(mmax,*), c(*) real(dp), intent(in) :: d(*), e(*), e2(*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl, klnew real(dp) :: tmp1, tmp2 ! Intrinsic Functions ! Executable Statements ! check for errors info = 0_${ik}$ if( ijob<1_${ik}$ .or. ijob>3_${ik}$ ) then info = -1_${ik}$ return end if ! initialize nab if( ijob==1_${ik}$ ) then ! compute the number of eigenvalues in the initial intervals. mout = 0_${ik}$ do ji = 1, minp do jp = 1, 2 tmp1 = d( 1_${ik}$ ) - ab( ji, jp ) if( abs( tmp1 )<pivmin )tmp1 = -pivmin nab( ji, jp ) = 0_${ik}$ if( tmp1<=zero )nab( ji, jp ) = 1_${ik}$ do j = 2, n tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp ) if( abs( tmp1 )<pivmin )tmp1 = -pivmin if( tmp1<=zero )nab( ji, jp ) = nab( ji, jp ) + 1_${ik}$ end do end do mout = mout + nab( ji, 2_${ik}$ ) - nab( ji, 1_${ik}$ ) end do return end if ! initialize for loop ! kf and kl have the following meaning: ! intervals 1,...,kf-1 have converged. ! intervals kf,...,kl still need to be refined. kf = 1_${ik}$ kl = minp ! if ijob=2, initialize c. ! if ijob=3, use the user-supplied starting point. if( ijob==2_${ik}$ ) then do ji = 1, minp c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) ) end do end if ! iteration loop loop_130: do jit = 1, nitmax ! loop over intervals if( kl-kf+1>=nbmin .and. nbmin>0_${ik}$ ) then ! begin of parallel version of the loop do ji = kf, kl ! compute n(c), the number of eigenvalues less than c work( ji ) = d( 1_${ik}$ ) - c( ji ) iwork( ji ) = 0_${ik}$ if( work( ji )<=pivmin ) then iwork( ji ) = 1_${ik}$ work( ji ) = min( work( ji ), -pivmin ) end if do j = 2, n work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji ) if( work( ji )<=pivmin ) then iwork( ji ) = iwork( ji ) + 1_${ik}$ work( ji ) = min( work( ji ), -pivmin ) end if end do end do if( ijob<=2_${ik}$ ) then ! ijob=2: choose all intervals containing eigenvalues. klnew = kl loop_70: do ji = kf, kl ! insure that n(w) is monotone iwork( ji ) = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), iwork( ji ) ) ) ! update the queue -- add intervals if both halves ! contain eigenvalues. if( iwork( ji )==nab( ji, 2_${ik}$ ) ) then ! no eigenvalue in the upper interval: ! just use the lower interval. ab( ji, 2_${ik}$ ) = c( ji ) else if( iwork( ji )==nab( ji, 1_${ik}$ ) ) then ! no eigenvalue in the lower interval: ! just use the upper interval. ab( ji, 1_${ik}$ ) = c( ji ) else klnew = klnew + 1_${ik}$ if( klnew<=mmax ) then ! eigenvalue in both intervals -- add upper to ! queue. ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ ) nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ ) ab( klnew, 1_${ik}$ ) = c( ji ) nab( klnew, 1_${ik}$ ) = iwork( ji ) ab( ji, 2_${ik}$ ) = c( ji ) nab( ji, 2_${ik}$ ) = iwork( ji ) else info = mmax + 1_${ik}$ end if end if end do loop_70 if( info/=0 )return kl = klnew else ! ijob=3: binary search. keep only the interval containing ! w s.t. n(w) = nval do ji = kf, kl if( iwork( ji )<=nval( ji ) ) then ab( ji, 1_${ik}$ ) = c( ji ) nab( ji, 1_${ik}$ ) = iwork( ji ) end if if( iwork( ji )>=nval( ji ) ) then ab( ji, 2_${ik}$ ) = c( ji ) nab( ji, 2_${ik}$ ) = iwork( ji ) end if end do end if else ! end of parallel version of the loop ! begin of serial version of the loop klnew = kl loop_100: do ji = kf, kl ! compute n(w), the number of eigenvalues less than w tmp1 = c( ji ) tmp2 = d( 1_${ik}$ ) - tmp1 itmp1 = 0_${ik}$ if( tmp2<=pivmin ) then itmp1 = 1_${ik}$ tmp2 = min( tmp2, -pivmin ) end if do j = 2, n tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1 if( tmp2<=pivmin ) then itmp1 = itmp1 + 1_${ik}$ tmp2 = min( tmp2, -pivmin ) end if end do if( ijob<=2_${ik}$ ) then ! ijob=2: choose all intervals containing eigenvalues. ! insure that n(w) is monotone itmp1 = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), itmp1 ) ) ! update the queue -- add intervals if both halves ! contain eigenvalues. if( itmp1==nab( ji, 2_${ik}$ ) ) then ! no eigenvalue in the upper interval: ! just use the lower interval. ab( ji, 2_${ik}$ ) = tmp1 else if( itmp1==nab( ji, 1_${ik}$ ) ) then ! no eigenvalue in the lower interval: ! just use the upper interval. ab( ji, 1_${ik}$ ) = tmp1 else if( klnew<mmax ) then ! eigenvalue in both intervals -- add upper to queue. klnew = klnew + 1_${ik}$ ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ ) nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ ) ab( klnew, 1_${ik}$ ) = tmp1 nab( klnew, 1_${ik}$ ) = itmp1 ab( ji, 2_${ik}$ ) = tmp1 nab( ji, 2_${ik}$ ) = itmp1 else info = mmax + 1_${ik}$ return end if else ! ijob=3: binary search. keep only the interval ! containing w s.t. n(w) = nval if( itmp1<=nval( ji ) ) then ab( ji, 1_${ik}$ ) = tmp1 nab( ji, 1_${ik}$ ) = itmp1 end if if( itmp1>=nval( ji ) ) then ab( ji, 2_${ik}$ ) = tmp1 nab( ji, 2_${ik}$ ) = itmp1 end if end if end do loop_100 kl = klnew end if ! check for convergence kfnew = kf loop_110: do ji = kf, kl tmp1 = abs( ab( ji, 2_${ik}$ )-ab( ji, 1_${ik}$ ) ) tmp2 = max( abs( ab( ji, 2_${ik}$ ) ), abs( ab( ji, 1_${ik}$ ) ) ) if( tmp1<max( abstol, pivmin, reltol*tmp2 ) .or.nab( ji, 1_${ik}$ )>=nab( ji, 2_${ik}$ ) ) & then ! converged -- swap with position kfnew, ! then increment kfnew if( ji>kfnew ) then tmp1 = ab( ji, 1_${ik}$ ) tmp2 = ab( ji, 2_${ik}$ ) itmp1 = nab( ji, 1_${ik}$ ) itmp2 = nab( ji, 2_${ik}$ ) ab( ji, 1_${ik}$ ) = ab( kfnew, 1_${ik}$ ) ab( ji, 2_${ik}$ ) = ab( kfnew, 2_${ik}$ ) nab( ji, 1_${ik}$ ) = nab( kfnew, 1_${ik}$ ) nab( ji, 2_${ik}$ ) = nab( kfnew, 2_${ik}$ ) ab( kfnew, 1_${ik}$ ) = tmp1 ab( kfnew, 2_${ik}$ ) = tmp2 nab( kfnew, 1_${ik}$ ) = itmp1 nab( kfnew, 2_${ik}$ ) = itmp2 if( ijob==3_${ik}$ ) then itmp1 = nval( ji ) nval( ji ) = nval( kfnew ) nval( kfnew ) = itmp1 end if end if kfnew = kfnew + 1_${ik}$ end if end do loop_110 kf = kfnew ! choose midpoints do ji = kf, kl c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) ) end do ! if no more intervals to refine, quit. if( kf>kl )go to 140 end do loop_130 ! converged 140 continue info = max( kl+1-kf, 0_${ik}$ ) mout = kl return end subroutine stdlib${ii}$_dlaebz #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$laebz( ijob, nitmax, n, mmax, minp, nbmin, abstol,reltol, pivmin, d, & !! DLAEBZ: contains the iteration loops which compute and use the !! function N(w), which is the count of eigenvalues of a symmetric !! tridiagonal matrix T less than or equal to its argument w. It !! performs a choice of two types of loops: !! IJOB=1, followed by !! IJOB=2: It takes as input a list of intervals and returns a list of !! sufficiently small intervals whose union contains the same !! eigenvalues as the union of the original intervals. !! The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. !! The output interval (AB(j,1),AB(j,2)] will contain !! eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. !! IJOB=3: It performs a binary search in each input interval !! (AB(j,1),AB(j,2)] for a point w(j) such that !! N(w(j))=NVAL(j), and uses C(j) as the starting point of !! the search. If such a w(j) is found, then on output !! AB(j,1)=AB(j,2)=w. If no such w(j) is found, then on output !! (AB(j,1),AB(j,2)] will be a small interval containing the !! point where N(w) jumps through NVAL(j), unless that point !! lies outside the initial interval. !! Note that the intervals are in all cases half-open intervals, !! i.e., of the form (a,b] , which includes b but not a . !! To avoid underflow, the matrix should be scaled so that its largest !! element is no greater than overflow**(1/2) * underflow**(1/4) !! in absolute value. To assure the most accurate computation !! of small eigenvalues, the matrix should be scaled to be !! not much smaller than that, either. !! See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal !! Matrix", Report CS41, Computer Science Dept., Stanford !! University, July 21, 1966 !! Note: the arguments are, in general, *not* checked for unreasonable !! values. e, e2, nval, ab, c, mout,nab, work, iwork, info ) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: ijob, minp, mmax, n, nbmin, nitmax integer(${ik}$), intent(out) :: info, mout real(${rk}$), intent(in) :: abstol, pivmin, reltol ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) integer(${ik}$), intent(inout) :: nab(mmax,*), nval(*) real(${rk}$), intent(inout) :: ab(mmax,*), c(*) real(${rk}$), intent(in) :: d(*), e(*), e2(*) real(${rk}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: itmp1, itmp2, j, ji, jit, jp, kf, kfnew, kl, klnew real(${rk}$) :: tmp1, tmp2 ! Intrinsic Functions ! Executable Statements ! check for errors info = 0_${ik}$ if( ijob<1_${ik}$ .or. ijob>3_${ik}$ ) then info = -1_${ik}$ return end if ! initialize nab if( ijob==1_${ik}$ ) then ! compute the number of eigenvalues in the initial intervals. mout = 0_${ik}$ do ji = 1, minp do jp = 1, 2 tmp1 = d( 1_${ik}$ ) - ab( ji, jp ) if( abs( tmp1 )<pivmin )tmp1 = -pivmin nab( ji, jp ) = 0_${ik}$ if( tmp1<=zero )nab( ji, jp ) = 1_${ik}$ do j = 2, n tmp1 = d( j ) - e2( j-1 ) / tmp1 - ab( ji, jp ) if( abs( tmp1 )<pivmin )tmp1 = -pivmin if( tmp1<=zero )nab( ji, jp ) = nab( ji, jp ) + 1_${ik}$ end do end do mout = mout + nab( ji, 2_${ik}$ ) - nab( ji, 1_${ik}$ ) end do return end if ! initialize for loop ! kf and kl have the following meaning: ! intervals 1,...,kf-1 have converged. ! intervals kf,...,kl still need to be refined. kf = 1_${ik}$ kl = minp ! if ijob=2, initialize c. ! if ijob=3, use the user-supplied starting point. if( ijob==2_${ik}$ ) then do ji = 1, minp c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) ) end do end if ! iteration loop loop_130: do jit = 1, nitmax ! loop over intervals if( kl-kf+1>=nbmin .and. nbmin>0_${ik}$ ) then ! begin of parallel version of the loop do ji = kf, kl ! compute n(c), the number of eigenvalues less than c work( ji ) = d( 1_${ik}$ ) - c( ji ) iwork( ji ) = 0_${ik}$ if( work( ji )<=pivmin ) then iwork( ji ) = 1_${ik}$ work( ji ) = min( work( ji ), -pivmin ) end if do j = 2, n work( ji ) = d( j ) - e2( j-1 ) / work( ji ) - c( ji ) if( work( ji )<=pivmin ) then iwork( ji ) = iwork( ji ) + 1_${ik}$ work( ji ) = min( work( ji ), -pivmin ) end if end do end do if( ijob<=2_${ik}$ ) then ! ijob=2: choose all intervals containing eigenvalues. klnew = kl loop_70: do ji = kf, kl ! insure that n(w) is monotone iwork( ji ) = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), iwork( ji ) ) ) ! update the queue -- add intervals if both halves ! contain eigenvalues. if( iwork( ji )==nab( ji, 2_${ik}$ ) ) then ! no eigenvalue in the upper interval: ! just use the lower interval. ab( ji, 2_${ik}$ ) = c( ji ) else if( iwork( ji )==nab( ji, 1_${ik}$ ) ) then ! no eigenvalue in the lower interval: ! just use the upper interval. ab( ji, 1_${ik}$ ) = c( ji ) else klnew = klnew + 1_${ik}$ if( klnew<=mmax ) then ! eigenvalue in both intervals -- add upper to ! queue. ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ ) nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ ) ab( klnew, 1_${ik}$ ) = c( ji ) nab( klnew, 1_${ik}$ ) = iwork( ji ) ab( ji, 2_${ik}$ ) = c( ji ) nab( ji, 2_${ik}$ ) = iwork( ji ) else info = mmax + 1_${ik}$ end if end if end do loop_70 if( info/=0 )return kl = klnew else ! ijob=3: binary search. keep only the interval containing ! w s.t. n(w) = nval do ji = kf, kl if( iwork( ji )<=nval( ji ) ) then ab( ji, 1_${ik}$ ) = c( ji ) nab( ji, 1_${ik}$ ) = iwork( ji ) end if if( iwork( ji )>=nval( ji ) ) then ab( ji, 2_${ik}$ ) = c( ji ) nab( ji, 2_${ik}$ ) = iwork( ji ) end if end do end if else ! end of parallel version of the loop ! begin of serial version of the loop klnew = kl loop_100: do ji = kf, kl ! compute n(w), the number of eigenvalues less than w tmp1 = c( ji ) tmp2 = d( 1_${ik}$ ) - tmp1 itmp1 = 0_${ik}$ if( tmp2<=pivmin ) then itmp1 = 1_${ik}$ tmp2 = min( tmp2, -pivmin ) end if do j = 2, n tmp2 = d( j ) - e2( j-1 ) / tmp2 - tmp1 if( tmp2<=pivmin ) then itmp1 = itmp1 + 1_${ik}$ tmp2 = min( tmp2, -pivmin ) end if end do if( ijob<=2_${ik}$ ) then ! ijob=2: choose all intervals containing eigenvalues. ! insure that n(w) is monotone itmp1 = min( nab( ji, 2_${ik}$ ),max( nab( ji, 1_${ik}$ ), itmp1 ) ) ! update the queue -- add intervals if both halves ! contain eigenvalues. if( itmp1==nab( ji, 2_${ik}$ ) ) then ! no eigenvalue in the upper interval: ! just use the lower interval. ab( ji, 2_${ik}$ ) = tmp1 else if( itmp1==nab( ji, 1_${ik}$ ) ) then ! no eigenvalue in the lower interval: ! just use the upper interval. ab( ji, 1_${ik}$ ) = tmp1 else if( klnew<mmax ) then ! eigenvalue in both intervals -- add upper to queue. klnew = klnew + 1_${ik}$ ab( klnew, 2_${ik}$ ) = ab( ji, 2_${ik}$ ) nab( klnew, 2_${ik}$ ) = nab( ji, 2_${ik}$ ) ab( klnew, 1_${ik}$ ) = tmp1 nab( klnew, 1_${ik}$ ) = itmp1 ab( ji, 2_${ik}$ ) = tmp1 nab( ji, 2_${ik}$ ) = itmp1 else info = mmax + 1_${ik}$ return end if else ! ijob=3: binary search. keep only the interval ! containing w s.t. n(w) = nval if( itmp1<=nval( ji ) ) then ab( ji, 1_${ik}$ ) = tmp1 nab( ji, 1_${ik}$ ) = itmp1 end if if( itmp1>=nval( ji ) ) then ab( ji, 2_${ik}$ ) = tmp1 nab( ji, 2_${ik}$ ) = itmp1 end if end if end do loop_100 kl = klnew end if ! check for convergence kfnew = kf loop_110: do ji = kf, kl tmp1 = abs( ab( ji, 2_${ik}$ )-ab( ji, 1_${ik}$ ) ) tmp2 = max( abs( ab( ji, 2_${ik}$ ) ), abs( ab( ji, 1_${ik}$ ) ) ) if( tmp1<max( abstol, pivmin, reltol*tmp2 ) .or.nab( ji, 1_${ik}$ )>=nab( ji, 2_${ik}$ ) ) & then ! converged -- swap with position kfnew, ! then increment kfnew if( ji>kfnew ) then tmp1 = ab( ji, 1_${ik}$ ) tmp2 = ab( ji, 2_${ik}$ ) itmp1 = nab( ji, 1_${ik}$ ) itmp2 = nab( ji, 2_${ik}$ ) ab( ji, 1_${ik}$ ) = ab( kfnew, 1_${ik}$ ) ab( ji, 2_${ik}$ ) = ab( kfnew, 2_${ik}$ ) nab( ji, 1_${ik}$ ) = nab( kfnew, 1_${ik}$ ) nab( ji, 2_${ik}$ ) = nab( kfnew, 2_${ik}$ ) ab( kfnew, 1_${ik}$ ) = tmp1 ab( kfnew, 2_${ik}$ ) = tmp2 nab( kfnew, 1_${ik}$ ) = itmp1 nab( kfnew, 2_${ik}$ ) = itmp2 if( ijob==3_${ik}$ ) then itmp1 = nval( ji ) nval( ji ) = nval( kfnew ) nval( kfnew ) = itmp1 end if end if kfnew = kfnew + 1_${ik}$ end if end do loop_110 kf = kfnew ! choose midpoints do ji = kf, kl c( ji ) = half*( ab( ji, 1_${ik}$ )+ab( ji, 2_${ik}$ ) ) end do ! if no more intervals to refine, quit. if( kf>kl )go to 140 end do loop_130 ! converged 140 continue info = max( kl+1-kf, 0_${ik}$ ) mout = kl return end subroutine stdlib${ii}$_${ri}$laebz #:endif #:endfor pure integer(${ik}$) module function stdlib${ii}$_slaneg( n, d, lld, sigma, pivmin, r ) !! SLANEG computes the Sturm count, the number of negative pivots !! encountered while factoring tridiagonal T - sigma I = L D L^T. !! This implementation works directly on the factors without forming !! the tridiagonal matrix T. The Sturm count is also the number of !! eigenvalues of T less than sigma. !! This routine is called from SLARRB. !! The current routine does not use the PIVMIN parameter but rather !! requires IEEE-754 propagation of Infinities and NaNs. This !! routine also has no input range restrictions but does require !! default exception handling such that x/0 produces Inf when x is !! non-zero, and Inf/Inf produces NaN. For more information, see: !! Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in !! Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on !! Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 !! (Tech report version in LAWN 172 with the same title.) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: n, r real(sp), intent(in) :: pivmin, sigma ! Array Arguments real(sp), intent(in) :: d(*), lld(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: blklen = 128_${ik}$ ! some architectures propagate infinities and nans very slowly, so ! the code computes counts in blklen chunks. then a nan can ! propagate at most blklen columns before being detected. this is ! not a general tuning parameter; it needs only to be just large ! enough that the overhead is tiny in common cases. ! Local Scalars integer(${ik}$) :: bj, j, neg1, neg2, negcnt real(sp) :: bsav, dminus, dplus, gamma, p, t, tmp logical(lk) :: sawnan ! Intrinsic Functions ! Executable Statements negcnt = 0_${ik}$ ! i) upper part: l d l^t - sigma i = l+ d+ l+^t t = -sigma loop_210: do bj = 1, r-1, blklen neg1 = 0_${ik}$ bsav = t do j = bj, min(bj+blklen-1, r-1) dplus = d( j ) + t if( dplus<zero ) neg1 = neg1 + 1_${ik}$ tmp = t / dplus t = tmp * lld( j ) - sigma end do sawnan = stdlib${ii}$_sisnan( t ) ! run a slower version of the above loop if a nan is detected. ! a nan should occur only with a zero pivot after an infinite ! pivot. in that case, substituting 1 for t/dplus is the ! correct limit. if( sawnan ) then neg1 = 0_${ik}$ t = bsav do j = bj, min(bj+blklen-1, r-1) dplus = d( j ) + t if( dplus<zero ) neg1 = neg1 + 1_${ik}$ tmp = t / dplus if (stdlib${ii}$_sisnan(tmp)) tmp = one t = tmp * lld(j) - sigma end do end if negcnt = negcnt + neg1 end do loop_210 ! ii) lower part: l d l^t - sigma i = u- d- u-^t p = d( n ) - sigma do bj = n-1, r, -blklen neg2 = 0_${ik}$ bsav = p do j = bj, max(bj-blklen+1, r), -1 dminus = lld( j ) + p if( dminus<zero ) neg2 = neg2 + 1_${ik}$ tmp = p / dminus p = tmp * d( j ) - sigma end do sawnan = stdlib${ii}$_sisnan( p ) ! as above, run a slower version that substitutes 1 for inf/inf. if( sawnan ) then neg2 = 0_${ik}$ p = bsav do j = bj, max(bj-blklen+1, r), -1 dminus = lld( j ) + p if( dminus<zero ) neg2 = neg2 + 1_${ik}$ tmp = p / dminus if (stdlib${ii}$_sisnan(tmp)) tmp = one p = tmp * d(j) - sigma end do end if negcnt = negcnt + neg2 end do ! iii) twist index ! t was shifted by sigma initially. gamma = (t + sigma) + p if( gamma<zero ) negcnt = negcnt+1 stdlib${ii}$_slaneg = negcnt end function stdlib${ii}$_slaneg pure integer(${ik}$) module function stdlib${ii}$_dlaneg( n, d, lld, sigma, pivmin, r ) !! DLANEG computes the Sturm count, the number of negative pivots !! encountered while factoring tridiagonal T - sigma I = L D L^T. !! This implementation works directly on the factors without forming !! the tridiagonal matrix T. The Sturm count is also the number of !! eigenvalues of T less than sigma. !! This routine is called from DLARRB. !! The current routine does not use the PIVMIN parameter but rather !! requires IEEE-754 propagation of Infinities and NaNs. This !! routine also has no input range restrictions but does require !! default exception handling such that x/0 produces Inf when x is !! non-zero, and Inf/Inf produces NaN. For more information, see: !! Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in !! Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on !! Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 !! (Tech report version in LAWN 172 with the same title.) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: n, r real(dp), intent(in) :: pivmin, sigma ! Array Arguments real(dp), intent(in) :: d(*), lld(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: blklen = 128_${ik}$ ! some architectures propagate infinities and nans very slowly, so ! the code computes counts in blklen chunks. then a nan can ! propagate at most blklen columns before being detected. this is ! not a general tuning parameter; it needs only to be just large ! enough that the overhead is tiny in common cases. ! Local Scalars integer(${ik}$) :: bj, j, neg1, neg2, negcnt real(dp) :: bsav, dminus, dplus, gamma, p, t, tmp logical(lk) :: sawnan ! Intrinsic Functions ! Executable Statements negcnt = 0_${ik}$ ! i) upper part: l d l^t - sigma i = l+ d+ l+^t t = -sigma loop_210: do bj = 1, r-1, blklen neg1 = 0_${ik}$ bsav = t do j = bj, min(bj+blklen-1, r-1) dplus = d( j ) + t if( dplus<zero ) neg1 = neg1 + 1_${ik}$ tmp = t / dplus t = tmp * lld( j ) - sigma end do sawnan = stdlib${ii}$_disnan( t ) ! run a slower version of the above loop if a nan is detected. ! a nan should occur only with a zero pivot after an infinite ! pivot. in that case, substituting 1 for t/dplus is the ! correct limit. if( sawnan ) then neg1 = 0_${ik}$ t = bsav do j = bj, min(bj+blklen-1, r-1) dplus = d( j ) + t if( dplus<zero ) neg1 = neg1 + 1_${ik}$ tmp = t / dplus if (stdlib${ii}$_disnan(tmp)) tmp = one t = tmp * lld(j) - sigma end do end if negcnt = negcnt + neg1 end do loop_210 ! ii) lower part: l d l^t - sigma i = u- d- u-^t p = d( n ) - sigma do bj = n-1, r, -blklen neg2 = 0_${ik}$ bsav = p do j = bj, max(bj-blklen+1, r), -1 dminus = lld( j ) + p if( dminus<zero ) neg2 = neg2 + 1_${ik}$ tmp = p / dminus p = tmp * d( j ) - sigma end do sawnan = stdlib${ii}$_disnan( p ) ! as above, run a slower version that substitutes 1 for inf/inf. if( sawnan ) then neg2 = 0_${ik}$ p = bsav do j = bj, max(bj-blklen+1, r), -1 dminus = lld( j ) + p if( dminus<zero ) neg2 = neg2 + 1_${ik}$ tmp = p / dminus if (stdlib${ii}$_disnan(tmp)) tmp = one p = tmp * d(j) - sigma end do end if negcnt = negcnt + neg2 end do ! iii) twist index ! t was shifted by sigma initially. gamma = (t + sigma) + p if( gamma<zero ) negcnt = negcnt+1 stdlib${ii}$_dlaneg = negcnt end function stdlib${ii}$_dlaneg #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure integer(${ik}$) module function stdlib${ii}$_${ri}$laneg( n, d, lld, sigma, pivmin, r ) !! DLANEG: computes the Sturm count, the number of negative pivots !! encountered while factoring tridiagonal T - sigma I = L D L^T. !! This implementation works directly on the factors without forming !! the tridiagonal matrix T. The Sturm count is also the number of !! eigenvalues of T less than sigma. !! This routine is called from DLARRB. !! The current routine does not use the PIVMIN parameter but rather !! requires IEEE-754 propagation of Infinities and NaNs. This !! routine also has no input range restrictions but does require !! default exception handling such that x/0 produces Inf when x is !! non-zero, and Inf/Inf produces NaN. For more information, see: !! Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in !! Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on !! Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 !! (Tech report version in LAWN 172 with the same title.) ! -- lapack auxiliary routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments integer(${ik}$), intent(in) :: n, r real(${rk}$), intent(in) :: pivmin, sigma ! Array Arguments real(${rk}$), intent(in) :: d(*), lld(*) ! ===================================================================== ! Parameters integer(${ik}$), parameter :: blklen = 128_${ik}$ ! some architectures propagate infinities and nans very slowly, so ! the code computes counts in blklen chunks. then a nan can ! propagate at most blklen columns before being detected. this is ! not a general tuning parameter; it needs only to be just large ! enough that the overhead is tiny in common cases. ! Local Scalars integer(${ik}$) :: bj, j, neg1, neg2, negcnt real(${rk}$) :: bsav, dminus, dplus, gamma, p, t, tmp logical(lk) :: sawnan ! Intrinsic Functions ! Executable Statements negcnt = 0_${ik}$ ! i) upper part: l d l^t - sigma i = l+ d+ l+^t t = -sigma loop_210: do bj = 1, r-1, blklen neg1 = 0_${ik}$ bsav = t do j = bj, min(bj+blklen-1, r-1) dplus = d( j ) + t if( dplus<zero ) neg1 = neg1 + 1_${ik}$ tmp = t / dplus t = tmp * lld( j ) - sigma end do sawnan = stdlib${ii}$_${ri}$isnan( t ) ! run a slower version of the above loop if a nan is detected. ! a nan should occur only with a zero pivot after an infinite ! pivot. in that case, substituting 1 for t/dplus is the ! correct limit. if( sawnan ) then neg1 = 0_${ik}$ t = bsav do j = bj, min(bj+blklen-1, r-1) dplus = d( j ) + t if( dplus<zero ) neg1 = neg1 + 1_${ik}$ tmp = t / dplus if (stdlib${ii}$_${ri}$isnan(tmp)) tmp = one t = tmp * lld(j) - sigma end do end if negcnt = negcnt + neg1 end do loop_210 ! ii) lower part: l d l^t - sigma i = u- d- u-^t p = d( n ) - sigma do bj = n-1, r, -blklen neg2 = 0_${ik}$ bsav = p do j = bj, max(bj-blklen+1, r), -1 dminus = lld( j ) + p if( dminus<zero ) neg2 = neg2 + 1_${ik}$ tmp = p / dminus p = tmp * d( j ) - sigma end do sawnan = stdlib${ii}$_${ri}$isnan( p ) ! as above, run a slower version that substitutes 1 for inf/inf. if( sawnan ) then neg2 = 0_${ik}$ p = bsav do j = bj, max(bj-blklen+1, r), -1 dminus = lld( j ) + p if( dminus<zero ) neg2 = neg2 + 1_${ik}$ tmp = p / dminus if (stdlib${ii}$_${ri}$isnan(tmp)) tmp = one p = tmp * d(j) - sigma end do end if negcnt = negcnt + neg2 end do ! iii) twist index ! t was shifted by sigma initially. gamma = (t + sigma) + p if( gamma<zero ) negcnt = negcnt+1 stdlib${ii}$_${ri}$laneg = negcnt end function stdlib${ii}$_${ri}$laneg #:endif #:endfor pure module subroutine stdlib${ii}$_slaed0( icompq, qsiz, n, d, e, q, ldq, qstore, ldqs,work, iwork, info & !! SLAED0 computes all eigenvalues and corresponding eigenvectors of a !! symmetric tridiagonal matrix using the divide and conquer method. ) ! -- 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 integer(${ik}$), intent(in) :: icompq, ldq, ldqs, n, qsiz integer(${ik}$), intent(out) :: info ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(inout) :: d(*), e(*), q(ldq,*) real(sp), intent(out) :: qstore(ldqs,*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, & iq, iqptr, iwrem, j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, subpbs, & tlvls real(sp) :: temp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( icompq<0_${ik}$ .or. icompq>2_${ik}$ ) then info = -1_${ik}$ else if( ( icompq==1_${ik}$ ) .and. ( qsiz<max( 0_${ik}$, n ) ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldqs<max( 1_${ik}$, n ) ) then info = -9_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SLAED0', -info ) return end if ! quick return if possible if( n==0 )return smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'SLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ ) ! determine the size and placement of the submatrices, and save in ! the leading elements of iwork. iwork( 1_${ik}$ ) = n subpbs = 1_${ik}$ tlvls = 0_${ik}$ 10 continue if( iwork( subpbs )>smlsiz ) then do j = subpbs, 1, -1 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$ iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$ end do tlvls = tlvls + 1_${ik}$ subpbs = 2_${ik}$*subpbs go to 10 end if do j = 2, subpbs iwork( j ) = iwork( j ) + iwork( j-1 ) end do ! divide the matrix into subpbs submatrices of size at most smlsiz+1 ! using rank-1 modifications (cuts). spm1 = subpbs - 1_${ik}$ do i = 1, spm1 submat = iwork( i ) + 1_${ik}$ smm1 = submat - 1_${ik}$ d( smm1 ) = d( smm1 ) - abs( e( smm1 ) ) d( submat ) = d( submat ) - abs( e( smm1 ) ) end do indxq = 4_${ik}$*n + 3_${ik}$ if( icompq/=2_${ik}$ ) then ! set up workspaces for eigenvalues only/accumulate new vectors ! routine temp = log( real( n,KIND=sp) ) / log( two ) lgn = int( temp,KIND=${ik}$) if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ iprmpt = indxq + n + 1_${ik}$ iperm = iprmpt + n*lgn iqptr = iperm + n*lgn igivpt = iqptr + n + 2_${ik}$ igivcl = igivpt + n*lgn igivnm = 1_${ik}$ iq = igivnm + 2_${ik}$*n*lgn iwrem = iq + n**2_${ik}$ + 1_${ik}$ ! initialize pointers do i = 0, subpbs iwork( iprmpt+i ) = 1_${ik}$ iwork( igivpt+i ) = 1_${ik}$ end do iwork( iqptr ) = 1_${ik}$ end if ! solve each submatrix eigenproblem at the bottom of the divide and ! conquer tree. curr = 0_${ik}$ loop_70: do i = 0, spm1 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 1_${ik}$ ) else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+1 ) - iwork( i ) end if if( icompq==2_${ik}$ ) then call stdlib${ii}$_ssteqr( 'I', matsiz, d( submat ), e( submat ),q( submat, submat ), & ldq, work, info ) if( info/=0 )go to 130 else call stdlib${ii}$_ssteqr( 'I', matsiz, d( submat ), e( submat ),work( iq-1+iwork( & iqptr+curr ) ), matsiz, work,info ) if( info/=0 )go to 130 if( icompq==1_${ik}$ ) then call stdlib${ii}$_sgemm( 'N', 'N', qsiz, matsiz, matsiz, one,q( 1_${ik}$, submat ), ldq, & work( iq-1+iwork( iqptr+curr ) ), matsiz, zero, qstore( 1_${ik}$, submat ),ldqs ) end if iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$ curr = curr + 1_${ik}$ end if k = 1_${ik}$ do j = submat, iwork( i+1 ) iwork( indxq+j ) = k k = k + 1_${ik}$ end do end do loop_70 ! successively merge eigensystems of adjacent submatrices ! into eigensystem for the corresponding larger matrix. ! while ( subpbs > 1 ) curlvl = 1_${ik}$ 80 continue if( subpbs>1_${ik}$ ) then spm2 = subpbs - 2_${ik}$ loop_90: do i = 0, spm2, 2 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 2_${ik}$ ) msd2 = iwork( 1_${ik}$ ) curprb = 0_${ik}$ else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+2 ) - iwork( i ) msd2 = matsiz / 2_${ik}$ curprb = curprb + 1_${ik}$ end if ! merge lower order eigensystems (of size msd2 and matsiz - msd2) ! into an eigensystem of size matsiz. ! stdlib${ii}$_slaed1 is used only for the full eigensystem of a tridiagonal ! matrix. ! stdlib${ii}$_slaed7 handles the cases in which eigenvalues only or eigenvalues ! and eigenvectors of a full symmetric matrix (which was reduced to ! tridiagonal form) are desired. if( icompq==2_${ik}$ ) then call stdlib${ii}$_slaed1( matsiz, d( submat ), q( submat, submat ),ldq, iwork( & indxq+submat ),e( submat+msd2-1 ), msd2, work,iwork( subpbs+1 ), info ) else call stdlib${ii}$_slaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,d( submat ), & qstore( 1_${ik}$, submat ), ldqs,iwork( indxq+submat ), e( submat+msd2-1 ),msd2, & work( iq ), iwork( iqptr ),iwork( iprmpt ), iwork( iperm ),iwork( igivpt ), & iwork( igivcl ),work( igivnm ), work( iwrem ),iwork( subpbs+1 ), info ) end if if( info/=0 )go to 130 iwork( i / 2_${ik}$+1 ) = iwork( i+2 ) end do loop_90 subpbs = subpbs / 2_${ik}$ curlvl = curlvl + 1_${ik}$ go to 80 end if ! end while ! re-merge the eigenvalues/vectors which were deflated at the final ! merge step. if( icompq==1_${ik}$ ) then do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) call stdlib${ii}$_scopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ ) end do call stdlib${ii}$_scopy( n, work, 1_${ik}$, d, 1_${ik}$ ) else if( icompq==2_${ik}$ ) then do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) call stdlib${ii}$_scopy( n, q( 1_${ik}$, j ), 1_${ik}$, work( n*i+1 ), 1_${ik}$ ) end do call stdlib${ii}$_scopy( n, work, 1_${ik}$, d, 1_${ik}$ ) call stdlib${ii}$_slacpy( 'A', n, n, work( n+1 ), n, q, ldq ) else do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) end do call stdlib${ii}$_scopy( n, work, 1_${ik}$, d, 1_${ik}$ ) end if go to 140 130 continue info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ 140 continue return end subroutine stdlib${ii}$_slaed0 pure module subroutine stdlib${ii}$_dlaed0( icompq, qsiz, n, d, e, q, ldq, qstore, ldqs,work, iwork, info & !! DLAED0 computes all eigenvalues and corresponding eigenvectors of a !! symmetric tridiagonal matrix using the divide and conquer method. ) ! -- 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 integer(${ik}$), intent(in) :: icompq, ldq, ldqs, n, qsiz integer(${ik}$), intent(out) :: info ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(inout) :: d(*), e(*), q(ldq,*) real(dp), intent(out) :: qstore(ldqs,*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, & iq, iqptr, iwrem, j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, subpbs, & tlvls real(dp) :: temp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( icompq<0_${ik}$ .or. icompq>2_${ik}$ ) then info = -1_${ik}$ else if( ( icompq==1_${ik}$ ) .and. ( qsiz<max( 0_${ik}$, n ) ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldqs<max( 1_${ik}$, n ) ) then info = -9_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DLAED0', -info ) return end if ! quick return if possible if( n==0 )return smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'DLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ ) ! determine the size and placement of the submatrices, and save in ! the leading elements of iwork. iwork( 1_${ik}$ ) = n subpbs = 1_${ik}$ tlvls = 0_${ik}$ 10 continue if( iwork( subpbs )>smlsiz ) then do j = subpbs, 1, -1 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$ iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$ end do tlvls = tlvls + 1_${ik}$ subpbs = 2_${ik}$*subpbs go to 10 end if do j = 2, subpbs iwork( j ) = iwork( j ) + iwork( j-1 ) end do ! divide the matrix into subpbs submatrices of size at most smlsiz+1 ! using rank-1 modifications (cuts). spm1 = subpbs - 1_${ik}$ do i = 1, spm1 submat = iwork( i ) + 1_${ik}$ smm1 = submat - 1_${ik}$ d( smm1 ) = d( smm1 ) - abs( e( smm1 ) ) d( submat ) = d( submat ) - abs( e( smm1 ) ) end do indxq = 4_${ik}$*n + 3_${ik}$ if( icompq/=2_${ik}$ ) then ! set up workspaces for eigenvalues only/accumulate new vectors ! routine temp = log( real( n,KIND=dp) ) / log( two ) lgn = int( temp,KIND=${ik}$) if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ iprmpt = indxq + n + 1_${ik}$ iperm = iprmpt + n*lgn iqptr = iperm + n*lgn igivpt = iqptr + n + 2_${ik}$ igivcl = igivpt + n*lgn igivnm = 1_${ik}$ iq = igivnm + 2_${ik}$*n*lgn iwrem = iq + n**2_${ik}$ + 1_${ik}$ ! initialize pointers do i = 0, subpbs iwork( iprmpt+i ) = 1_${ik}$ iwork( igivpt+i ) = 1_${ik}$ end do iwork( iqptr ) = 1_${ik}$ end if ! solve each submatrix eigenproblem at the bottom of the divide and ! conquer tree. curr = 0_${ik}$ loop_70: do i = 0, spm1 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 1_${ik}$ ) else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+1 ) - iwork( i ) end if if( icompq==2_${ik}$ ) then call stdlib${ii}$_dsteqr( 'I', matsiz, d( submat ), e( submat ),q( submat, submat ), & ldq, work, info ) if( info/=0 )go to 130 else call stdlib${ii}$_dsteqr( 'I', matsiz, d( submat ), e( submat ),work( iq-1+iwork( & iqptr+curr ) ), matsiz, work,info ) if( info/=0 )go to 130 if( icompq==1_${ik}$ ) then call stdlib${ii}$_dgemm( 'N', 'N', qsiz, matsiz, matsiz, one,q( 1_${ik}$, submat ), ldq, & work( iq-1+iwork( iqptr+curr ) ), matsiz, zero, qstore( 1_${ik}$, submat ),ldqs ) end if iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$ curr = curr + 1_${ik}$ end if k = 1_${ik}$ do j = submat, iwork( i+1 ) iwork( indxq+j ) = k k = k + 1_${ik}$ end do end do loop_70 ! successively merge eigensystems of adjacent submatrices ! into eigensystem for the corresponding larger matrix. ! while ( subpbs > 1 ) curlvl = 1_${ik}$ 80 continue if( subpbs>1_${ik}$ ) then spm2 = subpbs - 2_${ik}$ loop_90: do i = 0, spm2, 2 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 2_${ik}$ ) msd2 = iwork( 1_${ik}$ ) curprb = 0_${ik}$ else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+2 ) - iwork( i ) msd2 = matsiz / 2_${ik}$ curprb = curprb + 1_${ik}$ end if ! merge lower order eigensystems (of size msd2 and matsiz - msd2) ! into an eigensystem of size matsiz. ! stdlib${ii}$_dlaed1 is used only for the full eigensystem of a tridiagonal ! matrix. ! stdlib${ii}$_dlaed7 handles the cases in which eigenvalues only or eigenvalues ! and eigenvectors of a full symmetric matrix (which was reduced to ! tridiagonal form) are desired. if( icompq==2_${ik}$ ) then call stdlib${ii}$_dlaed1( matsiz, d( submat ), q( submat, submat ),ldq, iwork( & indxq+submat ),e( submat+msd2-1 ), msd2, work,iwork( subpbs+1 ), info ) else call stdlib${ii}$_dlaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,d( submat ), & qstore( 1_${ik}$, submat ), ldqs,iwork( indxq+submat ), e( submat+msd2-1 ),msd2, & work( iq ), iwork( iqptr ),iwork( iprmpt ), iwork( iperm ),iwork( igivpt ), & iwork( igivcl ),work( igivnm ), work( iwrem ),iwork( subpbs+1 ), info ) end if if( info/=0 )go to 130 iwork( i / 2_${ik}$+1 ) = iwork( i+2 ) end do loop_90 subpbs = subpbs / 2_${ik}$ curlvl = curlvl + 1_${ik}$ go to 80 end if ! end while ! re-merge the eigenvalues/vectors which were deflated at the final ! merge step. if( icompq==1_${ik}$ ) then do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) call stdlib${ii}$_dcopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ ) end do call stdlib${ii}$_dcopy( n, work, 1_${ik}$, d, 1_${ik}$ ) else if( icompq==2_${ik}$ ) then do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) call stdlib${ii}$_dcopy( n, q( 1_${ik}$, j ), 1_${ik}$, work( n*i+1 ), 1_${ik}$ ) end do call stdlib${ii}$_dcopy( n, work, 1_${ik}$, d, 1_${ik}$ ) call stdlib${ii}$_dlacpy( 'A', n, n, work( n+1 ), n, q, ldq ) else do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) end do call stdlib${ii}$_dcopy( n, work, 1_${ik}$, d, 1_${ik}$ ) end if go to 140 130 continue info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ 140 continue return end subroutine stdlib${ii}$_dlaed0 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$laed0( icompq, qsiz, n, d, e, q, ldq, qstore, ldqs,work, iwork, info & !! DLAED0: computes all eigenvalues and corresponding eigenvectors of a !! symmetric tridiagonal matrix using the divide and conquer method. ) ! -- 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 integer(${ik}$), intent(in) :: icompq, ldq, ldqs, n, qsiz integer(${ik}$), intent(out) :: info ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(${rk}$), intent(inout) :: d(*), e(*), q(ldq,*) real(${rk}$), intent(out) :: qstore(ldqs,*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, & iq, iqptr, iwrem, j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, subpbs, & tlvls real(${rk}$) :: temp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( icompq<0_${ik}$ .or. icompq>2_${ik}$ ) then info = -1_${ik}$ else if( ( icompq==1_${ik}$ ) .and. ( qsiz<max( 0_${ik}$, n ) ) ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -7_${ik}$ else if( ldqs<max( 1_${ik}$, n ) ) then info = -9_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DLAED0', -info ) return end if ! quick return if possible if( n==0 )return smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'DLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ ) ! determine the size and placement of the submatrices, and save in ! the leading elements of iwork. iwork( 1_${ik}$ ) = n subpbs = 1_${ik}$ tlvls = 0_${ik}$ 10 continue if( iwork( subpbs )>smlsiz ) then do j = subpbs, 1, -1 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$ iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$ end do tlvls = tlvls + 1_${ik}$ subpbs = 2_${ik}$*subpbs go to 10 end if do j = 2, subpbs iwork( j ) = iwork( j ) + iwork( j-1 ) end do ! divide the matrix into subpbs submatrices of size at most smlsiz+1 ! using rank-1 modifications (cuts). spm1 = subpbs - 1_${ik}$ do i = 1, spm1 submat = iwork( i ) + 1_${ik}$ smm1 = submat - 1_${ik}$ d( smm1 ) = d( smm1 ) - abs( e( smm1 ) ) d( submat ) = d( submat ) - abs( e( smm1 ) ) end do indxq = 4_${ik}$*n + 3_${ik}$ if( icompq/=2_${ik}$ ) then ! set up workspaces for eigenvalues only/accumulate new vectors ! routine temp = log( real( n,KIND=${rk}$) ) / log( two ) lgn = int( temp,KIND=${ik}$) if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ iprmpt = indxq + n + 1_${ik}$ iperm = iprmpt + n*lgn iqptr = iperm + n*lgn igivpt = iqptr + n + 2_${ik}$ igivcl = igivpt + n*lgn igivnm = 1_${ik}$ iq = igivnm + 2_${ik}$*n*lgn iwrem = iq + n**2_${ik}$ + 1_${ik}$ ! initialize pointers do i = 0, subpbs iwork( iprmpt+i ) = 1_${ik}$ iwork( igivpt+i ) = 1_${ik}$ end do iwork( iqptr ) = 1_${ik}$ end if ! solve each submatrix eigenproblem at the bottom of the divide and ! conquer tree. curr = 0_${ik}$ loop_70: do i = 0, spm1 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 1_${ik}$ ) else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+1 ) - iwork( i ) end if if( icompq==2_${ik}$ ) then call stdlib${ii}$_${ri}$steqr( 'I', matsiz, d( submat ), e( submat ),q( submat, submat ), & ldq, work, info ) if( info/=0 )go to 130 else call stdlib${ii}$_${ri}$steqr( 'I', matsiz, d( submat ), e( submat ),work( iq-1+iwork( & iqptr+curr ) ), matsiz, work,info ) if( info/=0 )go to 130 if( icompq==1_${ik}$ ) then call stdlib${ii}$_${ri}$gemm( 'N', 'N', qsiz, matsiz, matsiz, one,q( 1_${ik}$, submat ), ldq, & work( iq-1+iwork( iqptr+curr ) ), matsiz, zero, qstore( 1_${ik}$, submat ),ldqs ) end if iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$ curr = curr + 1_${ik}$ end if k = 1_${ik}$ do j = submat, iwork( i+1 ) iwork( indxq+j ) = k k = k + 1_${ik}$ end do end do loop_70 ! successively merge eigensystems of adjacent submatrices ! into eigensystem for the corresponding larger matrix. ! while ( subpbs > 1 ) curlvl = 1_${ik}$ 80 continue if( subpbs>1_${ik}$ ) then spm2 = subpbs - 2_${ik}$ loop_90: do i = 0, spm2, 2 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 2_${ik}$ ) msd2 = iwork( 1_${ik}$ ) curprb = 0_${ik}$ else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+2 ) - iwork( i ) msd2 = matsiz / 2_${ik}$ curprb = curprb + 1_${ik}$ end if ! merge lower order eigensystems (of size msd2 and matsiz - msd2) ! into an eigensystem of size matsiz. ! stdlib${ii}$_${ri}$laed1 is used only for the full eigensystem of a tridiagonal ! matrix. ! stdlib${ii}$_${ri}$laed7 handles the cases in which eigenvalues only or eigenvalues ! and eigenvectors of a full symmetric matrix (which was reduced to ! tridiagonal form) are desired. if( icompq==2_${ik}$ ) then call stdlib${ii}$_${ri}$laed1( matsiz, d( submat ), q( submat, submat ),ldq, iwork( & indxq+submat ),e( submat+msd2-1 ), msd2, work,iwork( subpbs+1 ), info ) else call stdlib${ii}$_${ri}$laed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,d( submat ), & qstore( 1_${ik}$, submat ), ldqs,iwork( indxq+submat ), e( submat+msd2-1 ),msd2, & work( iq ), iwork( iqptr ),iwork( iprmpt ), iwork( iperm ),iwork( igivpt ), & iwork( igivcl ),work( igivnm ), work( iwrem ),iwork( subpbs+1 ), info ) end if if( info/=0 )go to 130 iwork( i / 2_${ik}$+1 ) = iwork( i+2 ) end do loop_90 subpbs = subpbs / 2_${ik}$ curlvl = curlvl + 1_${ik}$ go to 80 end if ! end while ! re-merge the eigenvalues/vectors which were deflated at the final ! merge step. if( icompq==1_${ik}$ ) then do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) call stdlib${ii}$_${ri}$copy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ ) end do call stdlib${ii}$_${ri}$copy( n, work, 1_${ik}$, d, 1_${ik}$ ) else if( icompq==2_${ik}$ ) then do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) call stdlib${ii}$_${ri}$copy( n, q( 1_${ik}$, j ), 1_${ik}$, work( n*i+1 ), 1_${ik}$ ) end do call stdlib${ii}$_${ri}$copy( n, work, 1_${ik}$, d, 1_${ik}$ ) call stdlib${ii}$_${ri}$lacpy( 'A', n, n, work( n+1 ), n, q, ldq ) else do i = 1, n j = iwork( indxq+i ) work( i ) = d( j ) end do call stdlib${ii}$_${ri}$copy( n, work, 1_${ik}$, d, 1_${ik}$ ) end if go to 140 130 continue info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ 140 continue return end subroutine stdlib${ii}$_${ri}$laed0 #:endif #:endfor pure module subroutine stdlib${ii}$_claed0( qsiz, n, d, e, q, ldq, qstore, ldqs, rwork,iwork, info ) !! Using the divide and conquer method, CLAED0: computes all eigenvalues !! of a symmetric tridiagonal matrix which is one diagonal block of !! those from reducing a dense or band Hermitian matrix and !! corresponding eigenvectors of the dense or band matrix. ! -- 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 integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldq, ldqs, n, qsiz ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(inout) :: d(*), e(*) real(sp), intent(out) :: rwork(*) complex(sp), intent(inout) :: q(ldq,*) complex(sp), intent(out) :: qstore(ldqs,*) ! ===================================================================== ! warning: n could be as big as qsiz! ! Local Scalars integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, & iq, iqptr, iwrem, j, k, lgn, ll, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, & subpbs, tlvls real(sp) :: temp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ ! if( icompq < 0 .or. icompq > 2 ) then ! info = -1 ! else if( ( icompq == 1 ) .and. ( qsiz < max( 0, n ) ) ) ! $ then if( qsiz<max( 0_${ik}$, n ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldqs<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CLAED0', -info ) return end if ! quick return if possible if( n==0 )return smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'CLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ ) ! determine the size and placement of the submatrices, and save in ! the leading elements of iwork. iwork( 1_${ik}$ ) = n subpbs = 1_${ik}$ tlvls = 0_${ik}$ 10 continue if( iwork( subpbs )>smlsiz ) then do j = subpbs, 1, -1 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$ iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$ end do tlvls = tlvls + 1_${ik}$ subpbs = 2_${ik}$*subpbs go to 10 end if do j = 2, subpbs iwork( j ) = iwork( j ) + iwork( j-1 ) end do ! divide the matrix into subpbs submatrices of size at most smlsiz+1 ! using rank-1 modifications (cuts). spm1 = subpbs - 1_${ik}$ do i = 1, spm1 submat = iwork( i ) + 1_${ik}$ smm1 = submat - 1_${ik}$ d( smm1 ) = d( smm1 ) - abs( e( smm1 ) ) d( submat ) = d( submat ) - abs( e( smm1 ) ) end do indxq = 4_${ik}$*n + 3_${ik}$ ! set up workspaces for eigenvalues only/accumulate new vectors ! routine temp = log( real( n,KIND=sp) ) / log( two ) lgn = int( temp,KIND=${ik}$) if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ iprmpt = indxq + n + 1_${ik}$ iperm = iprmpt + n*lgn iqptr = iperm + n*lgn igivpt = iqptr + n + 2_${ik}$ igivcl = igivpt + n*lgn igivnm = 1_${ik}$ iq = igivnm + 2_${ik}$*n*lgn iwrem = iq + n**2_${ik}$ + 1_${ik}$ ! initialize pointers do i = 0, subpbs iwork( iprmpt+i ) = 1_${ik}$ iwork( igivpt+i ) = 1_${ik}$ end do iwork( iqptr ) = 1_${ik}$ ! solve each submatrix eigenproblem at the bottom of the divide and ! conquer tree. curr = 0_${ik}$ do i = 0, spm1 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 1_${ik}$ ) else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+1 ) - iwork( i ) end if ll = iq - 1_${ik}$ + iwork( iqptr+curr ) call stdlib${ii}$_ssteqr( 'I', matsiz, d( submat ), e( submat ),rwork( ll ), matsiz, & rwork, info ) call stdlib${ii}$_clacrm( qsiz, matsiz, q( 1_${ik}$, submat ), ldq, rwork( ll ),matsiz, qstore( & 1_${ik}$, submat ), ldqs,rwork( iwrem ) ) iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$ curr = curr + 1_${ik}$ if( info>0_${ik}$ ) then info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ return end if k = 1_${ik}$ do j = submat, iwork( i+1 ) iwork( indxq+j ) = k k = k + 1_${ik}$ end do end do ! successively merge eigensystems of adjacent submatrices ! into eigensystem for the corresponding larger matrix. ! while ( subpbs > 1 ) curlvl = 1_${ik}$ 80 continue if( subpbs>1_${ik}$ ) then spm2 = subpbs - 2_${ik}$ do i = 0, spm2, 2 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 2_${ik}$ ) msd2 = iwork( 1_${ik}$ ) curprb = 0_${ik}$ else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+2 ) - iwork( i ) msd2 = matsiz / 2_${ik}$ curprb = curprb + 1_${ik}$ end if ! merge lower order eigensystems (of size msd2 and matsiz - msd2) ! into an eigensystem of size matsiz. stdlib${ii}$_claed7 handles the case ! when the eigenvectors of a full or band hermitian matrix (which ! was reduced to tridiagonal form) are desired. ! i am free to use q as a valuable working space until loop 150. call stdlib${ii}$_claed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,d( submat ), & qstore( 1_${ik}$, submat ), ldqs,e( submat+msd2-1 ), iwork( indxq+submat ),rwork( iq ), & iwork( iqptr ), iwork( iprmpt ),iwork( iperm ), iwork( igivpt ),iwork( igivcl ), & rwork( igivnm ),q( 1_${ik}$, submat ), rwork( iwrem ),iwork( subpbs+1 ), info ) if( info>0_${ik}$ ) then info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ return end if iwork( i / 2_${ik}$+1 ) = iwork( i+2 ) end do subpbs = subpbs / 2_${ik}$ curlvl = curlvl + 1_${ik}$ go to 80 end if ! end while ! re-merge the eigenvalues/vectors which were deflated at the final ! merge step. do i = 1, n j = iwork( indxq+i ) rwork( i ) = d( j ) call stdlib${ii}$_ccopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ ) end do call stdlib${ii}$_scopy( n, rwork, 1_${ik}$, d, 1_${ik}$ ) return end subroutine stdlib${ii}$_claed0 pure module subroutine stdlib${ii}$_zlaed0( qsiz, n, d, e, q, ldq, qstore, ldqs, rwork,iwork, info ) !! Using the divide and conquer method, ZLAED0: computes all eigenvalues !! of a symmetric tridiagonal matrix which is one diagonal block of !! those from reducing a dense or band Hermitian matrix and !! corresponding eigenvectors of the dense or band matrix. ! -- 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 integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldq, ldqs, n, qsiz ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(inout) :: d(*), e(*) real(dp), intent(out) :: rwork(*) complex(dp), intent(inout) :: q(ldq,*) complex(dp), intent(out) :: qstore(ldqs,*) ! ===================================================================== ! warning: n could be as big as qsiz! ! Local Scalars integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, & iq, iqptr, iwrem, j, k, lgn, ll, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, & subpbs, tlvls real(dp) :: temp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ ! if( icompq < 0 .or. icompq > 2 ) then ! info = -1 ! else if( ( icompq == 1 ) .and. ( qsiz < max( 0, n ) ) ) ! $ then if( qsiz<max( 0_${ik}$, n ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldqs<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZLAED0', -info ) return end if ! quick return if possible if( n==0 )return smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'ZLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ ) ! determine the size and placement of the submatrices, and save in ! the leading elements of iwork. iwork( 1_${ik}$ ) = n subpbs = 1_${ik}$ tlvls = 0_${ik}$ 10 continue if( iwork( subpbs )>smlsiz ) then do j = subpbs, 1, -1 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$ iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$ end do tlvls = tlvls + 1_${ik}$ subpbs = 2_${ik}$*subpbs go to 10 end if do j = 2, subpbs iwork( j ) = iwork( j ) + iwork( j-1 ) end do ! divide the matrix into subpbs submatrices of size at most smlsiz+1 ! using rank-1 modifications (cuts). spm1 = subpbs - 1_${ik}$ do i = 1, spm1 submat = iwork( i ) + 1_${ik}$ smm1 = submat - 1_${ik}$ d( smm1 ) = d( smm1 ) - abs( e( smm1 ) ) d( submat ) = d( submat ) - abs( e( smm1 ) ) end do indxq = 4_${ik}$*n + 3_${ik}$ ! set up workspaces for eigenvalues only/accumulate new vectors ! routine temp = log( real( n,KIND=dp) ) / log( two ) lgn = int( temp,KIND=${ik}$) if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ iprmpt = indxq + n + 1_${ik}$ iperm = iprmpt + n*lgn iqptr = iperm + n*lgn igivpt = iqptr + n + 2_${ik}$ igivcl = igivpt + n*lgn igivnm = 1_${ik}$ iq = igivnm + 2_${ik}$*n*lgn iwrem = iq + n**2_${ik}$ + 1_${ik}$ ! initialize pointers do i = 0, subpbs iwork( iprmpt+i ) = 1_${ik}$ iwork( igivpt+i ) = 1_${ik}$ end do iwork( iqptr ) = 1_${ik}$ ! solve each submatrix eigenproblem at the bottom of the divide and ! conquer tree. curr = 0_${ik}$ do i = 0, spm1 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 1_${ik}$ ) else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+1 ) - iwork( i ) end if ll = iq - 1_${ik}$ + iwork( iqptr+curr ) call stdlib${ii}$_dsteqr( 'I', matsiz, d( submat ), e( submat ),rwork( ll ), matsiz, & rwork, info ) call stdlib${ii}$_zlacrm( qsiz, matsiz, q( 1_${ik}$, submat ), ldq, rwork( ll ),matsiz, qstore( & 1_${ik}$, submat ), ldqs,rwork( iwrem ) ) iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$ curr = curr + 1_${ik}$ if( info>0_${ik}$ ) then info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ return end if k = 1_${ik}$ do j = submat, iwork( i+1 ) iwork( indxq+j ) = k k = k + 1_${ik}$ end do end do ! successively merge eigensystems of adjacent submatrices ! into eigensystem for the corresponding larger matrix. ! while ( subpbs > 1 ) curlvl = 1_${ik}$ 80 continue if( subpbs>1_${ik}$ ) then spm2 = subpbs - 2_${ik}$ do i = 0, spm2, 2 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 2_${ik}$ ) msd2 = iwork( 1_${ik}$ ) curprb = 0_${ik}$ else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+2 ) - iwork( i ) msd2 = matsiz / 2_${ik}$ curprb = curprb + 1_${ik}$ end if ! merge lower order eigensystems (of size msd2 and matsiz - msd2) ! into an eigensystem of size matsiz. stdlib${ii}$_zlaed7 handles the case ! when the eigenvectors of a full or band hermitian matrix (which ! was reduced to tridiagonal form) are desired. ! i am free to use q as a valuable working space until loop 150. call stdlib${ii}$_zlaed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,d( submat ), & qstore( 1_${ik}$, submat ), ldqs,e( submat+msd2-1 ), iwork( indxq+submat ),rwork( iq ), & iwork( iqptr ), iwork( iprmpt ),iwork( iperm ), iwork( igivpt ),iwork( igivcl ), & rwork( igivnm ),q( 1_${ik}$, submat ), rwork( iwrem ),iwork( subpbs+1 ), info ) if( info>0_${ik}$ ) then info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ return end if iwork( i / 2_${ik}$+1 ) = iwork( i+2 ) end do subpbs = subpbs / 2_${ik}$ curlvl = curlvl + 1_${ik}$ go to 80 end if ! end while ! re-merge the eigenvalues/vectors which were deflated at the final ! merge step. do i = 1, n j = iwork( indxq+i ) rwork( i ) = d( j ) call stdlib${ii}$_zcopy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ ) end do call stdlib${ii}$_dcopy( n, rwork, 1_${ik}$, d, 1_${ik}$ ) return end subroutine stdlib${ii}$_zlaed0 #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$laed0( qsiz, n, d, e, q, ldq, qstore, ldqs, rwork,iwork, info ) !! Using the divide and conquer method, ZLAED0: computes all eigenvalues !! of a symmetric tridiagonal matrix which is one diagonal block of !! those from reducing a dense or band Hermitian matrix and !! corresponding eigenvectors of the dense or band matrix. ! -- 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 integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: ldq, ldqs, n, qsiz ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(${ck}$), intent(inout) :: d(*), e(*) real(${ck}$), intent(out) :: rwork(*) complex(${ck}$), intent(inout) :: q(ldq,*) complex(${ck}$), intent(out) :: qstore(ldqs,*) ! ===================================================================== ! warning: n could be as big as qsiz! ! Local Scalars integer(${ik}$) :: curlvl, curprb, curr, i, igivcl, igivnm, igivpt, indxq, iperm, iprmpt, & iq, iqptr, iwrem, j, k, lgn, ll, matsiz, msd2, smlsiz, smm1, spm1, spm2, submat, & subpbs, tlvls real(${ck}$) :: temp ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ ! if( icompq < 0 .or. icompq > 2 ) then ! info = -1 ! else if( ( icompq == 1 ) .and. ( qsiz < max( 0, n ) ) ) ! $ then if( qsiz<max( 0_${ik}$, n ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( ldqs<max( 1_${ik}$, n ) ) then info = -8_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZLAED0', -info ) return end if ! quick return if possible if( n==0 )return smlsiz = stdlib${ii}$_ilaenv( 9_${ik}$, 'ZLAED0', ' ', 0_${ik}$, 0_${ik}$, 0_${ik}$, 0_${ik}$ ) ! determine the size and placement of the submatrices, and save in ! the leading elements of iwork. iwork( 1_${ik}$ ) = n subpbs = 1_${ik}$ tlvls = 0_${ik}$ 10 continue if( iwork( subpbs )>smlsiz ) then do j = subpbs, 1, -1 iwork( 2_${ik}$*j ) = ( iwork( j )+1_${ik}$ ) / 2_${ik}$ iwork( 2_${ik}$*j-1 ) = iwork( j ) / 2_${ik}$ end do tlvls = tlvls + 1_${ik}$ subpbs = 2_${ik}$*subpbs go to 10 end if do j = 2, subpbs iwork( j ) = iwork( j ) + iwork( j-1 ) end do ! divide the matrix into subpbs submatrices of size at most smlsiz+1 ! using rank-1 modifications (cuts). spm1 = subpbs - 1_${ik}$ do i = 1, spm1 submat = iwork( i ) + 1_${ik}$ smm1 = submat - 1_${ik}$ d( smm1 ) = d( smm1 ) - abs( e( smm1 ) ) d( submat ) = d( submat ) - abs( e( smm1 ) ) end do indxq = 4_${ik}$*n + 3_${ik}$ ! set up workspaces for eigenvalues only/accumulate new vectors ! routine temp = log( real( n,KIND=${ck}$) ) / log( two ) lgn = int( temp,KIND=${ik}$) if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ if( 2_${ik}$**lgn<n )lgn = lgn + 1_${ik}$ iprmpt = indxq + n + 1_${ik}$ iperm = iprmpt + n*lgn iqptr = iperm + n*lgn igivpt = iqptr + n + 2_${ik}$ igivcl = igivpt + n*lgn igivnm = 1_${ik}$ iq = igivnm + 2_${ik}$*n*lgn iwrem = iq + n**2_${ik}$ + 1_${ik}$ ! initialize pointers do i = 0, subpbs iwork( iprmpt+i ) = 1_${ik}$ iwork( igivpt+i ) = 1_${ik}$ end do iwork( iqptr ) = 1_${ik}$ ! solve each submatrix eigenproblem at the bottom of the divide and ! conquer tree. curr = 0_${ik}$ do i = 0, spm1 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 1_${ik}$ ) else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+1 ) - iwork( i ) end if ll = iq - 1_${ik}$ + iwork( iqptr+curr ) call stdlib${ii}$_${c2ri(ci)}$steqr( 'I', matsiz, d( submat ), e( submat ),rwork( ll ), matsiz, & rwork, info ) call stdlib${ii}$_${ci}$lacrm( qsiz, matsiz, q( 1_${ik}$, submat ), ldq, rwork( ll ),matsiz, qstore( & 1_${ik}$, submat ), ldqs,rwork( iwrem ) ) iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2_${ik}$ curr = curr + 1_${ik}$ if( info>0_${ik}$ ) then info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ return end if k = 1_${ik}$ do j = submat, iwork( i+1 ) iwork( indxq+j ) = k k = k + 1_${ik}$ end do end do ! successively merge eigensystems of adjacent submatrices ! into eigensystem for the corresponding larger matrix. ! while ( subpbs > 1 ) curlvl = 1_${ik}$ 80 continue if( subpbs>1_${ik}$ ) then spm2 = subpbs - 2_${ik}$ do i = 0, spm2, 2 if( i==0_${ik}$ ) then submat = 1_${ik}$ matsiz = iwork( 2_${ik}$ ) msd2 = iwork( 1_${ik}$ ) curprb = 0_${ik}$ else submat = iwork( i ) + 1_${ik}$ matsiz = iwork( i+2 ) - iwork( i ) msd2 = matsiz / 2_${ik}$ curprb = curprb + 1_${ik}$ end if ! merge lower order eigensystems (of size msd2 and matsiz - msd2) ! into an eigensystem of size matsiz. stdlib${ii}$_${ci}$laed7 handles the case ! when the eigenvectors of a full or band hermitian matrix (which ! was reduced to tridiagonal form) are desired. ! i am free to use q as a valuable working space until loop 150. call stdlib${ii}$_${ci}$laed7( matsiz, msd2, qsiz, tlvls, curlvl, curprb,d( submat ), & qstore( 1_${ik}$, submat ), ldqs,e( submat+msd2-1 ), iwork( indxq+submat ),rwork( iq ), & iwork( iqptr ), iwork( iprmpt ),iwork( iperm ), iwork( igivpt ),iwork( igivcl ), & rwork( igivnm ),q( 1_${ik}$, submat ), rwork( iwrem ),iwork( subpbs+1 ), info ) if( info>0_${ik}$ ) then info = submat*( n+1 ) + submat + matsiz - 1_${ik}$ return end if iwork( i / 2_${ik}$+1 ) = iwork( i+2 ) end do subpbs = subpbs / 2_${ik}$ curlvl = curlvl + 1_${ik}$ go to 80 end if ! end while ! re-merge the eigenvalues/vectors which were deflated at the final ! merge step. do i = 1, n j = iwork( indxq+i ) rwork( i ) = d( j ) call stdlib${ii}$_${ci}$copy( qsiz, qstore( 1_${ik}$, j ), 1_${ik}$, q( 1_${ik}$, i ), 1_${ik}$ ) end do call stdlib${ii}$_${c2ri(ci)}$copy( n, rwork, 1_${ik}$, d, 1_${ik}$ ) return end subroutine stdlib${ii}$_${ci}$laed0 #:endif #:endfor pure module subroutine stdlib${ii}$_slaed1( n, d, q, ldq, indxq, rho, cutpnt, work, iwork,info ) !! SLAED1 computes the updated eigensystem of a diagonal !! matrix after modification by a rank-one symmetric matrix. This !! routine is used only for the eigenproblem which requires all !! eigenvalues and eigenvectors of a tridiagonal matrix. SLAED7 handles !! the case in which eigenvalues only or eigenvalues and eigenvectors !! of a full symmetric matrix (which was reduced to tridiagonal form) !! are desired. !! T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out) !! where Z = Q**T*u, u is a vector of length N with ones in the !! CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. !! The eigenvectors of the original matrix are stored in Q, and the !! eigenvalues are in D. The algorithm consists of three stages: !! The first stage consists of deflating the size of the problem !! when there are multiple eigenvalues or if there is a zero in !! the Z vector. For each such occurrence the dimension of the !! secular equation problem is reduced by one. This stage is !! performed by the routine SLAED2. !! The second stage consists of calculating the updated !! eigenvalues. This is done by finding the roots of the secular !! equation via the routine SLAED4 (as called by SLAED3). !! This routine also calculates the eigenvectors of the current !! problem. !! The final stage consists of computing the updated eigenvectors !! directly using the updated eigenvalues. The eigenvectors for !! the current problem are multiplied with the eigenvectors from !! the overall problem. ! -- 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 integer(${ik}$), intent(in) :: cutpnt, ldq, n integer(${ik}$), intent(out) :: info real(sp), intent(inout) :: rho ! Array Arguments integer(${ik}$), intent(inout) :: indxq(*) integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(inout) :: d(*), q(ldq,*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: coltyp, cpp1, i, idlmda, indx, indxc, indxp, iq2, is, iw, iz, k, n1, & n2 ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( n<0_${ik}$ ) then info = -1_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( min( 1_${ik}$, n / 2_${ik}$ )>cutpnt .or. ( n / 2_${ik}$ )<cutpnt ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SLAED1', -info ) return end if ! quick return if possible if( n==0 )return ! the following values are integer pointers which indicate ! the portion of the workspace ! used by a particular array in stdlib${ii}$_slaed2 and stdlib${ii}$_slaed3. iz = 1_${ik}$ idlmda = iz + n iw = idlmda + n iq2 = iw + n indx = 1_${ik}$ indxc = indx + n coltyp = indxc + n indxp = coltyp + n ! form the z-vector which consists of the last row of q_1 and the ! first row of q_2. call stdlib${ii}$_scopy( cutpnt, q( cutpnt, 1_${ik}$ ), ldq, work( iz ), 1_${ik}$ ) cpp1 = cutpnt + 1_${ik}$ call stdlib${ii}$_scopy( n-cutpnt, q( cpp1, cpp1 ), ldq, work( iz+cutpnt ), 1_${ik}$ ) ! deflate eigenvalues. call stdlib${ii}$_slaed2( k, n, cutpnt, d, q, ldq, indxq, rho, work( iz ),work( idlmda ), & work( iw ), work( iq2 ),iwork( indx ), iwork( indxc ), iwork( indxp ),iwork( coltyp ), & info ) if( info/=0 )go to 20 ! solve secular equation. if( k/=0_${ik}$ ) then is = ( iwork( coltyp )+iwork( coltyp+1 ) )*cutpnt +( iwork( coltyp+1 )+iwork( & coltyp+2 ) )*( n-cutpnt ) + iq2 call stdlib${ii}$_slaed3( k, n, cutpnt, d, q, ldq, rho, work( idlmda ),work( iq2 ), iwork(& indxc ), iwork( coltyp ),work( iw ), work( is ), info ) if( info/=0 )go to 20 ! prepare the indxq sorting permutation. n1 = k n2 = n - k call stdlib${ii}$_slamrg( n1, n2, d, 1_${ik}$, -1_${ik}$, indxq ) else do i = 1, n indxq( i ) = i end do end if 20 continue return end subroutine stdlib${ii}$_slaed1 pure module subroutine stdlib${ii}$_dlaed1( n, d, q, ldq, indxq, rho, cutpnt, work, iwork,info ) !! DLAED1 computes the updated eigensystem of a diagonal !! matrix after modification by a rank-one symmetric matrix. This !! routine is used only for the eigenproblem which requires all !! eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles !! the case in which eigenvalues only or eigenvalues and eigenvectors !! of a full symmetric matrix (which was reduced to tridiagonal form) !! are desired. !! T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out) !! where Z = Q**T*u, u is a vector of length N with ones in the !! CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. !! The eigenvectors of the original matrix are stored in Q, and the !! eigenvalues are in D. The algorithm consists of three stages: !! The first stage consists of deflating the size of the problem !! when there are multiple eigenvalues or if there is a zero in !! the Z vector. For each such occurrence the dimension of the !! secular equation problem is reduced by one. This stage is !! performed by the routine DLAED2. !! The second stage consists of calculating the updated !! eigenvalues. This is done by finding the roots of the secular !! equation via the routine DLAED4 (as called by DLAED3). !! This routine also calculates the eigenvectors of the current !! problem. !! The final stage consists of computing the updated eigenvectors !! directly using the updated eigenvalues. The eigenvectors for !! the current problem are multiplied with the eigenvectors from !! the overall problem. ! -- 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 integer(${ik}$), intent(in) :: cutpnt, ldq, n integer(${ik}$), intent(out) :: info real(dp), intent(inout) :: rho ! Array Arguments integer(${ik}$), intent(inout) :: indxq(*) integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(inout) :: d(*), q(ldq,*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: coltyp, i, idlmda, indx, indxc, indxp, iq2, is, iw, iz, k, n1, n2, & zpp1 ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( n<0_${ik}$ ) then info = -1_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( min( 1_${ik}$, n / 2_${ik}$ )>cutpnt .or. ( n / 2_${ik}$ )<cutpnt ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DLAED1', -info ) return end if ! quick return if possible if( n==0 )return ! the following values are integer pointers which indicate ! the portion of the workspace ! used by a particular array in stdlib${ii}$_dlaed2 and stdlib${ii}$_dlaed3. iz = 1_${ik}$ idlmda = iz + n iw = idlmda + n iq2 = iw + n indx = 1_${ik}$ indxc = indx + n coltyp = indxc + n indxp = coltyp + n ! form the z-vector which consists of the last row of q_1 and the ! first row of q_2. call stdlib${ii}$_dcopy( cutpnt, q( cutpnt, 1_${ik}$ ), ldq, work( iz ), 1_${ik}$ ) zpp1 = cutpnt + 1_${ik}$ call stdlib${ii}$_dcopy( n-cutpnt, q( zpp1, zpp1 ), ldq, work( iz+cutpnt ), 1_${ik}$ ) ! deflate eigenvalues. call stdlib${ii}$_dlaed2( k, n, cutpnt, d, q, ldq, indxq, rho, work( iz ),work( idlmda ), & work( iw ), work( iq2 ),iwork( indx ), iwork( indxc ), iwork( indxp ),iwork( coltyp ), & info ) if( info/=0 )go to 20 ! solve secular equation. if( k/=0_${ik}$ ) then is = ( iwork( coltyp )+iwork( coltyp+1 ) )*cutpnt +( iwork( coltyp+1 )+iwork( & coltyp+2 ) )*( n-cutpnt ) + iq2 call stdlib${ii}$_dlaed3( k, n, cutpnt, d, q, ldq, rho, work( idlmda ),work( iq2 ), iwork(& indxc ), iwork( coltyp ),work( iw ), work( is ), info ) if( info/=0 )go to 20 ! prepare the indxq sorting permutation. n1 = k n2 = n - k call stdlib${ii}$_dlamrg( n1, n2, d, 1_${ik}$, -1_${ik}$, indxq ) else do i = 1, n indxq( i ) = i end do end if 20 continue return end subroutine stdlib${ii}$_dlaed1 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$laed1( n, d, q, ldq, indxq, rho, cutpnt, work, iwork,info ) !! DLAED1: computes the updated eigensystem of a diagonal !! matrix after modification by a rank-one symmetric matrix. This !! routine is used only for the eigenproblem which requires all !! eigenvalues and eigenvectors of a tridiagonal matrix. DLAED7 handles !! the case in which eigenvalues only or eigenvalues and eigenvectors !! of a full symmetric matrix (which was reduced to tridiagonal form) !! are desired. !! T = Q(in) ( D(in) + RHO * Z*Z**T ) Q**T(in) = Q(out) * D(out) * Q**T(out) !! where Z = Q**T*u, u is a vector of length N with ones in the !! CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. !! The eigenvectors of the original matrix are stored in Q, and the !! eigenvalues are in D. The algorithm consists of three stages: !! The first stage consists of deflating the size of the problem !! when there are multiple eigenvalues or if there is a zero in !! the Z vector. For each such occurrence the dimension of the !! secular equation problem is reduced by one. This stage is !! performed by the routine DLAED2. !! The second stage consists of calculating the updated !! eigenvalues. This is done by finding the roots of the secular !! equation via the routine DLAED4 (as called by DLAED3). !! This routine also calculates the eigenvectors of the current !! problem. !! The final stage consists of computing the updated eigenvectors !! directly using the updated eigenvalues. The eigenvectors for !! the current problem are multiplied with the eigenvectors from !! the overall problem. ! -- 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 integer(${ik}$), intent(in) :: cutpnt, ldq, n integer(${ik}$), intent(out) :: info real(${rk}$), intent(inout) :: rho ! Array Arguments integer(${ik}$), intent(inout) :: indxq(*) integer(${ik}$), intent(out) :: iwork(*) real(${rk}$), intent(inout) :: d(*), q(ldq,*) real(${rk}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: coltyp, i, idlmda, indx, indxc, indxp, iq2, is, iw, iz, k, n1, n2, & zpp1 ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( n<0_${ik}$ ) then info = -1_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( min( 1_${ik}$, n / 2_${ik}$ )>cutpnt .or. ( n / 2_${ik}$ )<cutpnt ) then info = -7_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DLAED1', -info ) return end if ! quick return if possible if( n==0 )return ! the following values are integer pointers which indicate ! the portion of the workspace ! used by a particular array in stdlib${ii}$_${ri}$laed2 and stdlib${ii}$_${ri}$laed3. iz = 1_${ik}$ idlmda = iz + n iw = idlmda + n iq2 = iw + n indx = 1_${ik}$ indxc = indx + n coltyp = indxc + n indxp = coltyp + n ! form the z-vector which consists of the last row of q_1 and the ! first row of q_2. call stdlib${ii}$_${ri}$copy( cutpnt, q( cutpnt, 1_${ik}$ ), ldq, work( iz ), 1_${ik}$ ) zpp1 = cutpnt + 1_${ik}$ call stdlib${ii}$_${ri}$copy( n-cutpnt, q( zpp1, zpp1 ), ldq, work( iz+cutpnt ), 1_${ik}$ ) ! deflate eigenvalues. call stdlib${ii}$_${ri}$laed2( k, n, cutpnt, d, q, ldq, indxq, rho, work( iz ),work( idlmda ), & work( iw ), work( iq2 ),iwork( indx ), iwork( indxc ), iwork( indxp ),iwork( coltyp ), & info ) if( info/=0 )go to 20 ! solve secular equation. if( k/=0_${ik}$ ) then is = ( iwork( coltyp )+iwork( coltyp+1 ) )*cutpnt +( iwork( coltyp+1 )+iwork( & coltyp+2 ) )*( n-cutpnt ) + iq2 call stdlib${ii}$_${ri}$laed3( k, n, cutpnt, d, q, ldq, rho, work( idlmda ),work( iq2 ), iwork(& indxc ), iwork( coltyp ),work( iw ), work( is ), info ) if( info/=0 )go to 20 ! prepare the indxq sorting permutation. n1 = k n2 = n - k call stdlib${ii}$_${ri}$lamrg( n1, n2, d, 1_${ik}$, -1_${ik}$, indxq ) else do i = 1, n indxq( i ) = i end do end if 20 continue return end subroutine stdlib${ii}$_${ri}$laed1 #:endif #:endfor pure module subroutine stdlib${ii}$_slaed2( k, n, n1, d, q, ldq, indxq, rho, z, dlamda, w,q2, indx, indxc,& !! SLAED2 merges the two sets of eigenvalues together into a single !! sorted set. Then it tries to deflate the size of the problem. !! There are two ways in which deflation can occur: when two or more !! eigenvalues are close together or if there is a tiny entry in the !! Z vector. For each such occurrence the order of the related secular !! equation problem is reduced by one. indxp, coltyp, 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 integer(${ik}$), intent(out) :: info, k integer(${ik}$), intent(in) :: ldq, n, n1 real(sp), intent(inout) :: rho ! Array Arguments integer(${ik}$), intent(out) :: coltyp(*), indx(*), indxc(*), indxp(*) integer(${ik}$), intent(inout) :: indxq(*) real(sp), intent(inout) :: d(*), q(ldq,*), z(*) real(sp), intent(out) :: dlamda(*), q2(*), w(*) ! ===================================================================== ! Parameters real(sp), parameter :: mone = -1.0_sp ! Local Arrays integer(${ik}$) :: ctot(4_${ik}$), psm(4_${ik}$) ! Local Scalars integer(${ik}$) :: ct, i, imax, iq1, iq2, j, jmax, js, k2, n1p1, n2, nj, pj real(sp) :: c, eps, s, t, tau, tol ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ if( n<0_${ik}$ ) then info = -2_${ik}$ else if( ldq<max( 1_${ik}$, n ) ) then info = -6_${ik}$ else if( min( 1_${ik}$, ( n / 2_${ik}$ ) )>n1 .or. ( n / 2_${ik}$ )<n1 ) then info = -3_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SLAED2', -info ) return end if ! quick return if possible if( n==0 )return n2 = n - n1 n1p1 = n1 + 1_${ik}$ if( rho<zero ) then call stdlib${ii}$_sscal( n2, mone, z( n1p1 ), 1_${ik}$ ) end if ! normalize z so that norm(z) = 1. since z is the concatenation of ! two normalized vectors, norm2(z) = sqrt(2). t = one / sqrt( two ) call stdlib${ii}$_sscal( n, t, z, 1_${ik}$ ) ! rho = abs( norm(z)**2 * rho ) rho = abs( two*rho ) ! sort the eigenvalues into increasing order do i = n1p1, n indxq( i ) = indxq( i ) + n1 end do ! re-integrate the deflated parts from the last pass do i = 1, n dlamda( i ) = d( indxq( i ) ) end do call stdlib${ii}$_slamrg( n1, n2, dlamda, 1_${ik}$, 1_${ik}$, indxc ) do i = 1, n indx( i ) = indxq( indxc( i ) ) end do ! calculate the allowable deflation tolerance imax = stdlib${ii}$_isamax( n, z, 1_${ik}$ ) jmax = stdlib${ii}$_isamax( n, d, 1_${ik}$ ) eps = stdlib${ii}$_slamch( 'EPSILON' ) tol = eight*eps*max( abs( d( jmax ) ), abs( z( imax ) ) ) ! if the rank-1 modifier is small enough, no more needs to be done ! except to reorganize q so that its columns correspond with the ! elements in d. if( rho*abs( z( imax ) )<=tol ) then k = 0_${ik}$ iq2 = 1_${ik}$ do j = 1, n i = indx( j ) call stdlib${ii}$_scopy( n, q( 1_${ik}$, i ), 1_${ik}$, q2( iq2 ), 1_${ik}$ ) dlamda( j ) = d( i ) iq2 = iq2 + n end do call stdlib${ii}$_slacpy( 'A', n, n, q2, n, q, ldq ) call stdlib${ii}$_scopy( n, dlamda, 1_${ik}$, d, 1_${ik}$ ) go to 190 end if ! if there are multiple eigenvalues then the problem deflates. here ! the number of equal eigenvalues are found. as each equal ! eigenvalue is found, an elementary reflector is computed to rotate ! the corresponding eigensubspace so that the corresponding ! components of z are zero in this new basis. do i = 1, n1 coltyp( i ) = 1_${ik}$ end do do i = n1p1, n coltyp( i ) = 3_${ik}$ end do k = 0_${ik}$ k2 = n + 1_${ik}$ do j = 1, n nj = indx( j ) if( rho*abs( z( nj ) )<=tol ) then ! deflate due to small z component. k2 = k2 - 1_${ik}$ coltyp( nj ) = 4_${ik}$ indxp( k2 ) = nj if( j==n )go to 100 else pj = nj go to 80 end if end do 80 continue j = j + 1_${ik}$ nj = indx( j ) if( j>n )go to 100 if( rho*abs( z( nj ) )<=tol ) then ! deflate due to small z component. k2 = k2 - 1_${ik}$ coltyp( nj ) = 4_${ik}$ indxp( k2 ) = nj else ! check if eigenvalues are close enough to allow deflation. s = z( pj ) c = z( nj ) ! find sqrt(a**2+b**2) without overflow or ! destructive underflow. tau = stdlib${ii}$_slapy2( c, s ) t = d( nj ) - d( pj ) c = c / tau s = -s / tau if( abs( t*c*s )<=tol ) then ! deflation is possible. z( nj ) = tau z( pj ) = zero if( coltyp( nj )/=coltyp( pj ) )coltyp( nj ) = 2_${ik}$ coltyp( pj ) = 4_${ik}$ call stdlib${ii}$_srot( n, q( 1_${ik}$, pj ), 1_${ik}$, q( 1_${ik}$, nj ), 1_${ik}$, c, s ) t = d( pj )*c**2_${ik}$ + d( nj )*s**2_${ik}$ d( nj ) = d( pj )*s**2_${ik}$ + d( nj )*c**2_${ik}$ d( pj ) = t k2 = k2 - 1_${ik}$ i = 1_${ik}$ 90 continue if( k2+i<=n ) then if( d( pj )<d( indxp( k2+i ) ) ) then indxp( k2+i-1 ) = indxp( k2+i ) indxp( k2+i ) = pj i = i + 1_${ik}$ go to 90 else indxp( k2+i-1 ) = pj end if else indxp( k2+i-1 ) = pj end if pj = nj else k = k + 1_${ik}$ dlamda( k ) = d( pj ) w( k ) = z( pj ) indxp( k ) = pj pj = nj end if end if go to 80 100 continue ! record the last eigenvalue. k = k + 1_${ik}$ dlamda( k ) = d( pj ) w( k ) = z( pj ) indxp( k ) = pj ! count up the total number of the various types of columns, then ! form a permutation which positions the four column types into ! four uniform groups (although one or more of these groups may be ! empty). do j = 1, 4 ctot( j ) = 0_${ik}$ end do do j = 1, n ct = coltyp( j ) ctot( ct ) = ctot( ct ) + 1_${ik}$ end do ! psm(*) = position in submatrix (of types 1 through 4) psm( 1_${ik}$ ) = 1_${ik}$ psm( 2_${ik}$ ) = 1_${ik}$ + ctot( 1_${ik}$ ) psm( 3_${ik}$ ) = psm( 2_${ik}$ ) + ctot( 2_${ik}$ ) psm( 4_${ik}$ ) = psm( 3_${ik}$ ) + ctot( 3_${ik}$ ) k = n - ctot( 4_${ik}$ ) ! fill out the indxc array so that the permutation which it induces ! will place all type-1 columns first, all type-2 columns next, ! then all type-3's, and finally all type-4's. do j = 1, n js = indxp( j ) ct = coltyp( js ) indx( psm( ct ) ) = js indxc( psm( ct ) ) = j psm( ct ) = psm( ct ) + 1_${ik}$ end do ! sort the eigenvalues and corresponding eigenvectors into dlamda ! and q2 respectively. the eigenvalues/vectors which were not ! deflated go into the first k slots of dlamda and q2 respectively, ! while those which were deflated go into the last n - k slots. i = 1_${ik}$ iq1 = 1_${ik}$ iq2 = 1_${ik}$ + ( ctot( 1_${ik}$ )+ctot( 2_${ik}$ ) )*n1 do j = 1, ctot( 1 ) js = indx( i ) call stdlib${ii}$_scopy( n1, q( 1_${ik}$, js ), 1_${ik}$, q2( iq1 ), 1_${ik}$ ) z( i ) = d( js ) i = i + 1_${ik}$ iq1 = iq1 + n1 end do do j = 1, ctot( 2 ) js = indx( i ) call stdlib${ii}$_scopy( n1, q( 1_${ik}$, js ), 1_${ik}$, q2( iq1 ), 1_${ik}$ ) call stdlib${ii}$_scopy( n2, q( n1+1, js ), 1_${ik}$, q2( iq2 ), 1_${ik}$ ) z( i ) = d( js ) i = i + 1_${ik}$ iq1 = iq1 + n1 iq2 = iq2 + n2 end do do j = 1, ctot( 3 ) js = indx( i ) call stdlib${ii}$_scopy( n2, q( n1+1, js ), 1_${ik}$, q2( iq2 ), 1_${ik}$ ) z( i ) = d( js ) i = i + 1_${ik}$ iq2 = iq2 + n2 end do iq1 = iq2 do j = 1, ctot( 4 ) js = indx( i ) call stdlib${ii}$_scopy( n, q( 1_${ik}$, js ), 1_${ik}$, q2( iq2 ), 1_${ik}$ ) iq2 = iq2 + n z( i ) = d( js ) i = i + 1_${ik}$ end do ! the deflated eigenvalues and their corresponding vectors go back ! into the last n - k slots of d and q respectively. if( k<n ) then call stdlib${ii}$_slacpy( 'A', n, ctot( 4_${ik}$ ), q2( iq1 ), n,q( 1_${ik}$, k+1 ), ldq ) call stdlib${ii}$_scopy( n-k, z( k+1 ), 1_${ik}$, d( k+1 ), 1_${ik}$ ) end if ! copy ctot into coltyp for referencing in stdlib${ii}$_slaed3. do j = 1, 4 coltyp( j ) = ctot( j ) end do 190 continue return end subroutine stdlib${ii}$_slaed2 pure module subroutine stdlib${ii}$_dlaed2( k, n, n1, d, q, ldq, indxq, rho, z, dlamda, w,q2, indx, indxc,& !! DLAED2 merges the two sets of eigenvalues together into a single !! sorted set. Then it tries to deflate the size of the problem. !! There are two ways in which deflation can occur: when two or more !! eigenvalues are close together or if there is a tiny entry in the !! Z vector. For each such occurrence the order of the related secular !! equation problem is reduced by one. indxp, coltyp, 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 integer(${ik}$), intent(out) :: info, k integer(${ik}$), intent(in) :: ldq, n, n1 real(dp), intent(inout) :: rho ! Array Arguments integer(${ik}$), intent(out) :: coltyp(*), indx(*), indxc(*), indxp(*) integer(${ik}$), intent(inout) :: indxq(*) real(dp), intent(inout) :: d(*), q(ldq,*), z(*) real(dp), intent(out) :: dlamda(*), q2(*), w(*) ! ===================================================================== ! Parameters real(dp), parameter :: mone = -1.0_dp ! Local Arrays integer(${ik}$) :: ctot(4_${ik}$