#:include "common.fypp" submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_svd_comp implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_sgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info ) !! SGEBRD reduces a general real M-by-N matrix A to upper or lower !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, lwork, m, n ! Array Arguments real(sp), intent(inout) :: a(lda,*) real(sp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'SGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) lwkopt = ( m+n )*nb work( 1_${ik}$ ) = real( lwkopt,KIND=sp) lquery = ( lwork==-1_${ik}$ ) if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then info = -10_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGEBRD', -info ) return else if( lquery ) then return end if ! quick return if possible minmn = min( m, n ) if( minmn==0_${ik}$ ) then work( 1_${ik}$ ) = 1_${ik}$ return end if ws = max( m, n ) ldwrkx = m ldwrky = n if( nb>1_${ik}$ .and. nb<minmn ) then ! set the crossover point nx. nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'SGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) ! determine when to switch from blocked to unblocked code. if( nx<minmn ) then ws = ( m+n )*nb if( lwork<ws ) then ! not enough work space for the optimal nb, consider using ! a smaller block size. nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'SGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( lwork>=( m+n )*nbmin ) then nb = lwork / ( m+n ) else nb = 1_${ik}$ nx = minmn end if end if end if else nx = minmn end if do i = 1, minmn - nx, nb ! reduce rows and columns i:i+nb-1 to bidiagonal form and return ! the matrices x and y which are needed to update the unreduced ! part of the matrix call stdlib${ii}$_slabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), & taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky ) ! update the trailing submatrix a(i+nb:m,i+nb:n), using an update ! of the form a := a - v*y**t - x*u**t call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, a( i+& nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, one,a( i+nb, i+nb ), lda ) call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, & work( nb+1 ), ldwrkx, a( i, i+nb ), lda,one, a( i+nb, i+nb ), lda ) ! copy diagonal and off-diagonal elements of b back into a if( m>=n ) then do j = i, i + nb - 1 a( j, j ) = d( j ) a( j, j+1 ) = e( j ) end do else do j = i, i + nb - 1 a( j, j ) = d( j ) a( j+1, j ) = e( j ) end do end if end do ! use unblocked code to reduce the remainder of the matrix call stdlib${ii}$_sgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), & work, iinfo ) work( 1_${ik}$ ) = ws return end subroutine stdlib${ii}$_sgebrd pure module subroutine stdlib${ii}$_dgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info ) !! DGEBRD reduces a general real M-by-N matrix A to upper or lower !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, lwork, m, n ! Array Arguments real(dp), intent(inout) :: a(lda,*) real(dp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) lwkopt = ( m+n )*nb work( 1_${ik}$ ) = real( lwkopt,KIND=dp) lquery = ( lwork==-1_${ik}$ ) if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then info = -10_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGEBRD', -info ) return else if( lquery ) then return end if ! quick return if possible minmn = min( m, n ) if( minmn==0_${ik}$ ) then work( 1_${ik}$ ) = 1_${ik}$ return end if ws = max( m, n ) ldwrkx = m ldwrky = n if( nb>1_${ik}$ .and. nb<minmn ) then ! set the crossover point nx. nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) ! determine when to switch from blocked to unblocked code. if( nx<minmn ) then ws = ( m+n )*nb if( lwork<ws ) then ! not enough work space for the optimal nb, consider using ! a smaller block size. nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( lwork>=( m+n )*nbmin ) then nb = lwork / ( m+n ) else nb = 1_${ik}$ nx = minmn end if end if end if else nx = minmn end if do i = 1, minmn - nx, nb ! reduce rows and columns i:i+nb-1 to bidiagonal form and return ! the matrices x and y which are needed to update the unreduced ! part of the matrix call stdlib${ii}$_dlabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), & taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky ) ! update the trailing submatrix a(i+nb:m,i+nb:n), using an update ! of the form a := a - v*y**t - x*u**t call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, a( i+& nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, one,a( i+nb, i+nb ), lda ) call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, & work( nb+1 ), ldwrkx, a( i, i+nb ), lda,one, a( i+nb, i+nb ), lda ) ! copy diagonal and off-diagonal elements of b back into a if( m>=n ) then do j = i, i + nb - 1 a( j, j ) = d( j ) a( j, j+1 ) = e( j ) end do else do j = i, i + nb - 1 a( j, j ) = d( j ) a( j+1, j ) = e( j ) end do end if end do ! use unblocked code to reduce the remainder of the matrix call stdlib${ii}$_dgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), & work, iinfo ) work( 1_${ik}$ ) = ws return end subroutine stdlib${ii}$_dgebrd #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$gebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info ) !! DGEBRD: reduces a general real M-by-N matrix A to upper or lower !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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(out) :: info integer(${ik}$), intent(in) :: lda, lwork, m, n ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) lwkopt = ( m+n )*nb work( 1_${ik}$ ) = real( lwkopt,KIND=${rk}$) lquery = ( lwork==-1_${ik}$ ) if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then info = -10_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGEBRD', -info ) return else if( lquery ) then return end if ! quick return if possible minmn = min( m, n ) if( minmn==0_${ik}$ ) then work( 1_${ik}$ ) = 1_${ik}$ return end if ws = max( m, n ) ldwrkx = m ldwrky = n if( nb>1_${ik}$ .and. nb<minmn ) then ! set the crossover point nx. nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) ! determine when to switch from blocked to unblocked code. if( nx<minmn ) then ws = ( m+n )*nb if( lwork<ws ) then ! not enough work space for the optimal nb, consider using ! a smaller block size. nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'DGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( lwork>=( m+n )*nbmin ) then nb = lwork / ( m+n ) else nb = 1_${ik}$ nx = minmn end if end if end if else nx = minmn end if do i = 1, minmn - nx, nb ! reduce rows and columns i:i+nb-1 to bidiagonal form and return ! the matrices x and y which are needed to update the unreduced ! part of the matrix call stdlib${ii}$_${ri}$labrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), & taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky ) ! update the trailing submatrix a(i+nb:m,i+nb:n), using an update ! of the form a := a - v*y**t - x*u**t call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, a( i+& nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, one,a( i+nb, i+nb ), lda ) call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -one, & work( nb+1 ), ldwrkx, a( i, i+nb ), lda,one, a( i+nb, i+nb ), lda ) ! copy diagonal and off-diagonal elements of b back into a if( m>=n ) then do j = i, i + nb - 1 a( j, j ) = d( j ) a( j, j+1 ) = e( j ) end do else do j = i, i + nb - 1 a( j, j ) = d( j ) a( j+1, j ) = e( j ) end do end if end do ! use unblocked code to reduce the remainder of the matrix call stdlib${ii}$_${ri}$gebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), & work, iinfo ) work( 1_${ik}$ ) = ws return end subroutine stdlib${ii}$_${ri}$gebrd #:endif #:endfor pure module subroutine stdlib${ii}$_cgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info ) !! CGEBRD reduces a general complex M-by-N matrix A to upper or lower !! bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, lwork, m, n ! Array Arguments real(sp), intent(out) :: d(*), e(*) complex(sp), intent(inout) :: a(lda,*) complex(sp), intent(out) :: taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'CGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) lwkopt = ( m+n )*nb work( 1_${ik}$ ) = real( lwkopt,KIND=sp) lquery = ( lwork==-1_${ik}$ ) if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then info = -10_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGEBRD', -info ) return else if( lquery ) then return end if ! quick return if possible minmn = min( m, n ) if( minmn==0_${ik}$ ) then work( 1_${ik}$ ) = 1_${ik}$ return end if ws = max( m, n ) ldwrkx = m ldwrky = n if( nb>1_${ik}$ .and. nb<minmn ) then ! set the crossover point nx. nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'CGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) ! determine when to switch from blocked to unblocked code. if( nx<minmn ) then ws = ( m+n )*nb if( lwork<ws ) then ! not enough work space for the optimal nb, consider using ! a smaller block size. nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'CGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( lwork>=( m+n )*nbmin ) then nb = lwork / ( m+n ) else nb = 1_${ik}$ nx = minmn end if end if end if else nx = minmn end if do i = 1, minmn - nx, nb ! reduce rows and columns i:i+ib-1 to bidiagonal form and return ! the matrices x and y which are needed to update the unreduced ! part of the matrix call stdlib${ii}$_clabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), & taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky ) ! update the trailing submatrix a(i+ib:m,i+ib:n), using ! an update of the form a := a - v*y**h - x*u**h call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE', m-i-nb+1,n-i-nb+1, nb, -& cone, a( i+nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, cone,a( i+nb, i+nb ), lda ) call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -cone, & work( nb+1 ), ldwrkx, a( i, i+nb ), lda,cone, a( i+nb, i+nb ), lda ) ! copy diagonal and off-diagonal elements of b back into a if( m>=n ) then do j = i, i + nb - 1 a( j, j ) = d( j ) a( j, j+1 ) = e( j ) end do else do j = i, i + nb - 1 a( j, j ) = d( j ) a( j+1, j ) = e( j ) end do end if end do ! use unblocked code to reduce the remainder of the matrix call stdlib${ii}$_cgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), & work, iinfo ) work( 1_${ik}$ ) = ws return end subroutine stdlib${ii}$_cgebrd pure module subroutine stdlib${ii}$_zgebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info ) !! ZGEBRD reduces a general complex M-by-N matrix A to upper or lower !! bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, lwork, m, n ! Array Arguments real(dp), intent(out) :: d(*), e(*) complex(dp), intent(inout) :: a(lda,*) complex(dp), intent(out) :: taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) lwkopt = ( m+n )*nb work( 1_${ik}$ ) = real( lwkopt,KIND=dp) lquery = ( lwork==-1_${ik}$ ) if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then info = -10_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGEBRD', -info ) return else if( lquery ) then return end if ! quick return if possible minmn = min( m, n ) if( minmn==0_${ik}$ ) then work( 1_${ik}$ ) = 1_${ik}$ return end if ws = max( m, n ) ldwrkx = m ldwrky = n if( nb>1_${ik}$ .and. nb<minmn ) then ! set the crossover point nx. nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) ! determine when to switch from blocked to unblocked code. if( nx<minmn ) then ws = ( m+n )*nb if( lwork<ws ) then ! not enough work space for the optimal nb, consider using ! a smaller block size. nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( lwork>=( m+n )*nbmin ) then nb = lwork / ( m+n ) else nb = 1_${ik}$ nx = minmn end if end if end if else nx = minmn end if do i = 1, minmn - nx, nb ! reduce rows and columns i:i+ib-1 to bidiagonal form and return ! the matrices x and y which are needed to update the unreduced ! part of the matrix call stdlib${ii}$_zlabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), & taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky ) ! update the trailing submatrix a(i+ib:m,i+ib:n), using ! an update of the form a := a - v*y**h - x*u**h call stdlib${ii}$_zgemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE', m-i-nb+1,n-i-nb+1, nb, -& cone, a( i+nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, cone,a( i+nb, i+nb ), lda ) call stdlib${ii}$_zgemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -cone, & work( nb+1 ), ldwrkx, a( i, i+nb ), lda,cone, a( i+nb, i+nb ), lda ) ! copy diagonal and off-diagonal elements of b back into a if( m>=n ) then do j = i, i + nb - 1 a( j, j ) = d( j ) a( j, j+1 ) = e( j ) end do else do j = i, i + nb - 1 a( j, j ) = d( j ) a( j+1, j ) = e( j ) end do end if end do ! use unblocked code to reduce the remainder of the matrix call stdlib${ii}$_zgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), & work, iinfo ) work( 1_${ik}$ ) = ws return end subroutine stdlib${ii}$_zgebrd #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$gebrd( m, n, a, lda, d, e, tauq, taup, work, lwork,info ) !! ZGEBRD: reduces a general complex M-by-N matrix A to upper or lower !! bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, lwork, m, n ! Array Arguments real(${ck}$), intent(out) :: d(*), e(*) complex(${ck}$), intent(inout) :: a(lda,*) complex(${ck}$), intent(out) :: taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: lquery integer(${ik}$) :: i, iinfo, j, ldwrkx, ldwrky, lwkopt, minmn, nb, nbmin, nx, ws ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ nb = max( 1_${ik}$, stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) lwkopt = ( m+n )*nb work( 1_${ik}$ ) = real( lwkopt,KIND=${ck}$) lquery = ( lwork==-1_${ik}$ ) if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ else if( lwork<max( 1_${ik}$, m, n ) .and. .not.lquery ) then info = -10_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGEBRD', -info ) return else if( lquery ) then return end if ! quick return if possible minmn = min( m, n ) if( minmn==0_${ik}$ ) then work( 1_${ik}$ ) = 1_${ik}$ return end if ws = max( m, n ) ldwrkx = m ldwrky = n if( nb>1_${ik}$ .and. nb<minmn ) then ! set the crossover point nx. nx = max( nb, stdlib${ii}$_ilaenv( 3_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) ) ! determine when to switch from blocked to unblocked code. if( nx<minmn ) then ws = ( m+n )*nb if( lwork<ws ) then ! not enough work space for the optimal nb, consider using ! a smaller block size. nbmin = stdlib${ii}$_ilaenv( 2_${ik}$, 'ZGEBRD', ' ', m, n, -1_${ik}$, -1_${ik}$ ) if( lwork>=( m+n )*nbmin ) then nb = lwork / ( m+n ) else nb = 1_${ik}$ nx = minmn end if end if end if else nx = minmn end if do i = 1, minmn - nx, nb ! reduce rows and columns i:i+ib-1 to bidiagonal form and return ! the matrices x and y which are needed to update the unreduced ! part of the matrix call stdlib${ii}$_${ci}$labrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),tauq( i ), & taup( i ), work, ldwrkx,work( ldwrkx*nb+1 ), ldwrky ) ! update the trailing submatrix a(i+ib:m,i+ib:n), using ! an update of the form a := a - v*y**h - x*u**h call stdlib${ii}$_${ci}$gemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE', m-i-nb+1,n-i-nb+1, nb, -& cone, a( i+nb, i ), lda,work( ldwrkx*nb+nb+1 ), ldwrky, cone,a( i+nb, i+nb ), lda ) call stdlib${ii}$_${ci}$gemm( 'NO TRANSPOSE', 'NO TRANSPOSE', m-i-nb+1, n-i-nb+1,nb, -cone, & work( nb+1 ), ldwrkx, a( i, i+nb ), lda,cone, a( i+nb, i+nb ), lda ) ! copy diagonal and off-diagonal elements of b back into a if( m>=n ) then do j = i, i + nb - 1 a( j, j ) = d( j ) a( j, j+1 ) = e( j ) end do else do j = i, i + nb - 1 a( j, j ) = d( j ) a( j+1, j ) = e( j ) end do end if end do ! use unblocked code to reduce the remainder of the matrix call stdlib${ii}$_${ci}$gebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),tauq( i ), taup( i ), & work, iinfo ) work( 1_${ik}$ ) = ws return end subroutine stdlib${ii}$_${ci}$gebrd #:endif #:endfor pure module subroutine stdlib${ii}$_sgebd2( m, n, a, lda, d, e, tauq, taup, work, info ) !! SGEBD2 reduces a real general m by n matrix A to upper or lower !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, m, n ! Array Arguments real(sp), intent(inout) :: a(lda,*) real(sp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGEBD2', -info ) return end if if( m>=n ) then ! reduce to upper bidiagonal form do i = 1, n ! generate elementary reflector h(i) to annihilate a(i+1:m,i) call stdlib${ii}$_slarfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) ) d( i ) = a( i, i ) a( i, i ) = one ! apply h(i) to a(i:m,i+1:n) from the left if( i<n )call stdlib${ii}$_slarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$, tauq( i ),a( i, i+& 1_${ik}$ ), lda, work ) a( i, i ) = d( i ) if( i<n ) then ! generate elementary reflector g(i) to annihilate ! a(i,i+2:n) call stdlib${ii}$_slarfg( n-i, a( i, i+1 ), a( i, min( i+2, n ) ),lda, taup( i ) ) e( i ) = a( i, i+1 ) a( i, i+1 ) = one ! apply g(i) to a(i+1:m,i+1:n) from the right call stdlib${ii}$_slarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 & ), lda, work ) a( i, i+1 ) = e( i ) else taup( i ) = zero end if end do else ! reduce to lower bidiagonal form do i = 1, m ! generate elementary reflector g(i) to annihilate a(i,i+1:n) call stdlib${ii}$_slarfg( n-i+1, a( i, i ), a( i, min( i+1, n ) ), lda,taup( i ) ) d( i ) = a( i, i ) a( i, i ) = one ! apply g(i) to a(i+1:m,i:n) from the right if( i<m )call stdlib${ii}$_slarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+& 1_${ik}$, i ), lda, work ) a( i, i ) = d( i ) if( i<m ) then ! generate elementary reflector h(i) to annihilate ! a(i+2:m,i) call stdlib${ii}$_slarfg( m-i, a( i+1, i ), a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) ) e( i ) = a( i+1, i ) a( i+1, i ) = one ! apply h(i) to a(i+1:m,i+1:n) from the left call stdlib${ii}$_slarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$, tauq( i ),a( i+1, i+1 ), & lda, work ) a( i+1, i ) = e( i ) else tauq( i ) = zero end if end do end if return end subroutine stdlib${ii}$_sgebd2 pure module subroutine stdlib${ii}$_dgebd2( m, n, a, lda, d, e, tauq, taup, work, info ) !! DGEBD2 reduces a real general m by n matrix A to upper or lower !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, m, n ! Array Arguments real(dp), intent(inout) :: a(lda,*) real(dp), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGEBD2', -info ) return end if if( m>=n ) then ! reduce to upper bidiagonal form do i = 1, n ! generate elementary reflector h(i) to annihilate a(i+1:m,i) call stdlib${ii}$_dlarfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) ) d( i ) = a( i, i ) a( i, i ) = one ! apply h(i) to a(i:m,i+1:n) from the left if( i<n )call stdlib${ii}$_dlarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$, tauq( i ),a( i, i+& 1_${ik}$ ), lda, work ) a( i, i ) = d( i ) if( i<n ) then ! generate elementary reflector g(i) to annihilate ! a(i,i+2:n) call stdlib${ii}$_dlarfg( n-i, a( i, i+1 ), a( i, min( i+2, n ) ),lda, taup( i ) ) e( i ) = a( i, i+1 ) a( i, i+1 ) = one ! apply g(i) to a(i+1:m,i+1:n) from the right call stdlib${ii}$_dlarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 & ), lda, work ) a( i, i+1 ) = e( i ) else taup( i ) = zero end if end do else ! reduce to lower bidiagonal form do i = 1, m ! generate elementary reflector g(i) to annihilate a(i,i+1:n) call stdlib${ii}$_dlarfg( n-i+1, a( i, i ), a( i, min( i+1, n ) ), lda,taup( i ) ) d( i ) = a( i, i ) a( i, i ) = one ! apply g(i) to a(i+1:m,i:n) from the right if( i<m )call stdlib${ii}$_dlarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+& 1_${ik}$, i ), lda, work ) a( i, i ) = d( i ) if( i<m ) then ! generate elementary reflector h(i) to annihilate ! a(i+2:m,i) call stdlib${ii}$_dlarfg( m-i, a( i+1, i ), a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) ) e( i ) = a( i+1, i ) a( i+1, i ) = one ! apply h(i) to a(i+1:m,i+1:n) from the left call stdlib${ii}$_dlarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$, tauq( i ),a( i+1, i+1 ), & lda, work ) a( i+1, i ) = e( i ) else tauq( i ) = zero end if end do end if return end subroutine stdlib${ii}$_dgebd2 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$gebd2( m, n, a, lda, d, e, tauq, taup, work, info ) !! DGEBD2: reduces a real general m by n matrix A to upper or lower !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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(out) :: info integer(${ik}$), intent(in) :: lda, m, n ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(out) :: d(*), e(*), taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGEBD2', -info ) return end if if( m>=n ) then ! reduce to upper bidiagonal form do i = 1, n ! generate elementary reflector h(i) to annihilate a(i+1:m,i) call stdlib${ii}$_${ri}$larfg( m-i+1, a( i, i ), a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) ) d( i ) = a( i, i ) a( i, i ) = one ! apply h(i) to a(i:m,i+1:n) from the left if( i<n )call stdlib${ii}$_${ri}$larf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$, tauq( i ),a( i, i+& 1_${ik}$ ), lda, work ) a( i, i ) = d( i ) if( i<n ) then ! generate elementary reflector g(i) to annihilate ! a(i,i+2:n) call stdlib${ii}$_${ri}$larfg( n-i, a( i, i+1 ), a( i, min( i+2, n ) ),lda, taup( i ) ) e( i ) = a( i, i+1 ) a( i, i+1 ) = one ! apply g(i) to a(i+1:m,i+1:n) from the right call stdlib${ii}$_${ri}$larf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 & ), lda, work ) a( i, i+1 ) = e( i ) else taup( i ) = zero end if end do else ! reduce to lower bidiagonal form do i = 1, m ! generate elementary reflector g(i) to annihilate a(i,i+1:n) call stdlib${ii}$_${ri}$larfg( n-i+1, a( i, i ), a( i, min( i+1, n ) ), lda,taup( i ) ) d( i ) = a( i, i ) a( i, i ) = one ! apply g(i) to a(i+1:m,i:n) from the right if( i<m )call stdlib${ii}$_${ri}$larf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+& 1_${ik}$, i ), lda, work ) a( i, i ) = d( i ) if( i<m ) then ! generate elementary reflector h(i) to annihilate ! a(i+2:m,i) call stdlib${ii}$_${ri}$larfg( m-i, a( i+1, i ), a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) ) e( i ) = a( i+1, i ) a( i+1, i ) = one ! apply h(i) to a(i+1:m,i+1:n) from the left call stdlib${ii}$_${ri}$larf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$, tauq( i ),a( i+1, i+1 ), & lda, work ) a( i+1, i ) = e( i ) else tauq( i ) = zero end if end do end if return end subroutine stdlib${ii}$_${ri}$gebd2 #:endif #:endfor pure module subroutine stdlib${ii}$_cgebd2( m, n, a, lda, d, e, tauq, taup, work, info ) !! CGEBD2 reduces a complex general m by n matrix A to upper or lower !! real bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, m, n ! Array Arguments real(sp), intent(out) :: d(*), e(*) complex(sp), intent(inout) :: a(lda,*) complex(sp), intent(out) :: taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i complex(sp) :: alpha ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGEBD2', -info ) return end if if( m>=n ) then ! reduce to upper bidiagonal form do i = 1, n ! generate elementary reflector h(i) to annihilate a(i+1:m,i) alpha = a( i, i ) call stdlib${ii}$_clarfg( m-i+1, alpha, a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) ) d( i ) = real( alpha,KIND=sp) a( i, i ) = cone ! apply h(i)**h to a(i:m,i+1:n) from the left if( i<n )call stdlib${ii}$_clarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$,conjg( tauq( i ) ), & a( i, i+1 ), lda, work ) a( i, i ) = d( i ) if( i<n ) then ! generate elementary reflector g(i) to annihilate ! a(i,i+2:n) call stdlib${ii}$_clacgv( n-i, a( i, i+1 ), lda ) alpha = a( i, i+1 ) call stdlib${ii}$_clarfg( n-i, alpha, a( i, min( i+2, n ) ),lda, taup( i ) ) e( i ) = real( alpha,KIND=sp) a( i, i+1 ) = cone ! apply g(i) to a(i+1:m,i+1:n) from the right call stdlib${ii}$_clarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 & ), lda, work ) call stdlib${ii}$_clacgv( n-i, a( i, i+1 ), lda ) a( i, i+1 ) = e( i ) else taup( i ) = czero end if end do else ! reduce to lower bidiagonal form do i = 1, m ! generate elementary reflector g(i) to annihilate a(i,i+1:n) call stdlib${ii}$_clacgv( n-i+1, a( i, i ), lda ) alpha = a( i, i ) call stdlib${ii}$_clarfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,taup( i ) ) d( i ) = real( alpha,KIND=sp) a( i, i ) = cone ! apply g(i) to a(i+1:m,i:n) from the right if( i<m )call stdlib${ii}$_clarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+& 1_${ik}$, i ), lda, work ) call stdlib${ii}$_clacgv( n-i+1, a( i, i ), lda ) a( i, i ) = d( i ) if( i<m ) then ! generate elementary reflector h(i) to annihilate ! a(i+2:m,i) alpha = a( i+1, i ) call stdlib${ii}$_clarfg( m-i, alpha, a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) ) e( i ) = real( alpha,KIND=sp) a( i+1, i ) = cone ! apply h(i)**h to a(i+1:m,i+1:n) from the left call stdlib${ii}$_clarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$,conjg( tauq( i ) ), a( i+& 1_${ik}$, i+1 ), lda,work ) a( i+1, i ) = e( i ) else tauq( i ) = czero end if end do end if return end subroutine stdlib${ii}$_cgebd2 pure module subroutine stdlib${ii}$_zgebd2( m, n, a, lda, d, e, tauq, taup, work, info ) !! ZGEBD2 reduces a complex general m by n matrix A to upper or lower !! real bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, m, n ! Array Arguments real(dp), intent(out) :: d(*), e(*) complex(dp), intent(inout) :: a(lda,*) complex(dp), intent(out) :: taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i complex(dp) :: alpha ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGEBD2', -info ) return end if if( m>=n ) then ! reduce to upper bidiagonal form do i = 1, n ! generate elementary reflector h(i) to annihilate a(i+1:m,i) alpha = a( i, i ) call stdlib${ii}$_zlarfg( m-i+1, alpha, a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) ) d( i ) = real( alpha,KIND=dp) a( i, i ) = cone ! apply h(i)**h to a(i:m,i+1:n) from the left if( i<n )call stdlib${ii}$_zlarf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$,conjg( tauq( i ) ), & a( i, i+1 ), lda, work ) a( i, i ) = d( i ) if( i<n ) then ! generate elementary reflector g(i) to annihilate ! a(i,i+2:n) call stdlib${ii}$_zlacgv( n-i, a( i, i+1 ), lda ) alpha = a( i, i+1 ) call stdlib${ii}$_zlarfg( n-i, alpha, a( i, min( i+2, n ) ), lda,taup( i ) ) e( i ) = real( alpha,KIND=dp) a( i, i+1 ) = cone ! apply g(i) to a(i+1:m,i+1:n) from the right call stdlib${ii}$_zlarf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 & ), lda, work ) call stdlib${ii}$_zlacgv( n-i, a( i, i+1 ), lda ) a( i, i+1 ) = e( i ) else taup( i ) = czero end if end do else ! reduce to lower bidiagonal form do i = 1, m ! generate elementary reflector g(i) to annihilate a(i,i+1:n) call stdlib${ii}$_zlacgv( n-i+1, a( i, i ), lda ) alpha = a( i, i ) call stdlib${ii}$_zlarfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,taup( i ) ) d( i ) = real( alpha,KIND=dp) a( i, i ) = cone ! apply g(i) to a(i+1:m,i:n) from the right if( i<m )call stdlib${ii}$_zlarf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+& 1_${ik}$, i ), lda, work ) call stdlib${ii}$_zlacgv( n-i+1, a( i, i ), lda ) a( i, i ) = d( i ) if( i<m ) then ! generate elementary reflector h(i) to annihilate ! a(i+2:m,i) alpha = a( i+1, i ) call stdlib${ii}$_zlarfg( m-i, alpha, a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) ) e( i ) = real( alpha,KIND=dp) a( i+1, i ) = cone ! apply h(i)**h to a(i+1:m,i+1:n) from the left call stdlib${ii}$_zlarf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$,conjg( tauq( i ) ), a( i+& 1_${ik}$, i+1 ), lda,work ) a( i+1, i ) = e( i ) else tauq( i ) = czero end if end do end if return end subroutine stdlib${ii}$_zgebd2 #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$gebd2( m, n, a, lda, d, e, tauq, taup, work, info ) !! ZGEBD2: reduces a complex general m by n matrix A to upper or lower !! real bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. ! -- 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) :: lda, m, n ! Array Arguments real(${ck}$), intent(out) :: d(*), e(*) complex(${ck}$), intent(inout) :: a(lda,*) complex(${ck}$), intent(out) :: taup(*), tauq(*), work(*) ! ===================================================================== ! Local Scalars integer(${ik}$) :: i complex(${ck}$) :: alpha ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ if( m<0_${ik}$ ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, m ) ) then info = -4_${ik}$ end if if( info<0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGEBD2', -info ) return end if if( m>=n ) then ! reduce to upper bidiagonal form do i = 1, n ! generate elementary reflector h(i) to annihilate a(i+1:m,i) alpha = a( i, i ) call stdlib${ii}$_${ci}$larfg( m-i+1, alpha, a( min( i+1, m ), i ), 1_${ik}$,tauq( i ) ) d( i ) = real( alpha,KIND=${ck}$) a( i, i ) = cone ! apply h(i)**h to a(i:m,i+1:n) from the left if( i<n )call stdlib${ii}$_${ci}$larf( 'LEFT', m-i+1, n-i, a( i, i ), 1_${ik}$,conjg( tauq( i ) ), & a( i, i+1 ), lda, work ) a( i, i ) = d( i ) if( i<n ) then ! generate elementary reflector g(i) to annihilate ! a(i,i+2:n) call stdlib${ii}$_${ci}$lacgv( n-i, a( i, i+1 ), lda ) alpha = a( i, i+1 ) call stdlib${ii}$_${ci}$larfg( n-i, alpha, a( i, min( i+2, n ) ), lda,taup( i ) ) e( i ) = real( alpha,KIND=${ck}$) a( i, i+1 ) = cone ! apply g(i) to a(i+1:m,i+1:n) from the right call stdlib${ii}$_${ci}$larf( 'RIGHT', m-i, n-i, a( i, i+1 ), lda,taup( i ), a( i+1, i+1 & ), lda, work ) call stdlib${ii}$_${ci}$lacgv( n-i, a( i, i+1 ), lda ) a( i, i+1 ) = e( i ) else taup( i ) = czero end if end do else ! reduce to lower bidiagonal form do i = 1, m ! generate elementary reflector g(i) to annihilate a(i,i+1:n) call stdlib${ii}$_${ci}$lacgv( n-i+1, a( i, i ), lda ) alpha = a( i, i ) call stdlib${ii}$_${ci}$larfg( n-i+1, alpha, a( i, min( i+1, n ) ), lda,taup( i ) ) d( i ) = real( alpha,KIND=${ck}$) a( i, i ) = cone ! apply g(i) to a(i+1:m,i:n) from the right if( i<m )call stdlib${ii}$_${ci}$larf( 'RIGHT', m-i, n-i+1, a( i, i ), lda,taup( i ), a( i+& 1_${ik}$, i ), lda, work ) call stdlib${ii}$_${ci}$lacgv( n-i+1, a( i, i ), lda ) a( i, i ) = d( i ) if( i<m ) then ! generate elementary reflector h(i) to annihilate ! a(i+2:m,i) alpha = a( i+1, i ) call stdlib${ii}$_${ci}$larfg( m-i, alpha, a( min( i+2, m ), i ), 1_${ik}$,tauq( i ) ) e( i ) = real( alpha,KIND=${ck}$) a( i+1, i ) = cone ! apply h(i)**h to a(i+1:m,i+1:n) from the left call stdlib${ii}$_${ci}$larf( 'LEFT', m-i, n-i, a( i+1, i ), 1_${ik}$,conjg( tauq( i ) ), a( i+& 1_${ik}$, i+1 ), lda,work ) a( i+1, i ) = e( i ) else tauq( i ) = czero end if end do end if return end subroutine stdlib${ii}$_${ci}$gebd2 #:endif #:endfor pure module subroutine stdlib${ii}$_sgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, & !! SGBBRD reduces a real general m-by-n band matrix A to upper !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! The routine computes B, and optionally forms Q or P**T, or computes !! Q**T*C for a given matrix C. ldc, work, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc ! Array Arguments real(sp), intent(inout) :: ab(ldab,*), c(ldc,*) real(sp), intent(out) :: d(*), e(*), pt(ldpt,*), q(ldq,*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: wantb, wantc, wantpt, wantq integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mn,& mu, mu0, nr, nrt real(sp) :: ra, rb, rc, rs ! Intrinsic Functions ! Executable Statements ! test the input parameters wantb = stdlib_lsame( vect, 'B' ) wantq = stdlib_lsame( vect, 'Q' ) .or. wantb wantpt = stdlib_lsame( vect, 'P' ) .or. wantb wantc = ncc>0_${ik}$ klu1 = kl + ku + 1_${ik}$ info = 0_${ik}$ if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then info = -1_${ik}$ else if( m<0_${ik}$ ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ncc<0_${ik}$ ) then info = -4_${ik}$ else if( kl<0_${ik}$ ) then info = -5_${ik}$ else if( ku<0_${ik}$ ) then info = -6_${ik}$ else if( ldab<klu1 ) then info = -8_${ik}$ else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then info = -12_${ik}$ else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then info = -14_${ik}$ else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SGBBRD', -info ) return end if ! initialize q and p**t to the unit matrix, if needed if( wantq )call stdlib${ii}$_slaset( 'FULL', m, m, zero, one, q, ldq ) if( wantpt )call stdlib${ii}$_slaset( 'FULL', n, n, zero, one, pt, ldpt ) ! quick return if possible. if( m==0 .or. n==0 )return minmn = min( m, n ) if( kl+ku>1_${ik}$ ) then ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce ! first to lower bidiagonal form and then transform to upper ! bidiagonal if( ku>0_${ik}$ ) then ml0 = 1_${ik}$ mu0 = 2_${ik}$ else ml0 = 2_${ik}$ mu0 = 1_${ik}$ end if ! wherever possible, plane rotations are generated and applied in ! vector operations of length nr over the index set j1:j2:klu1. ! the sines of the plane rotations are stored in work(1:max(m,n)) ! and the cosines in work(max(m,n)+1:2*max(m,n)). mn = max( m, n ) klm = min( m-1, kl ) kun = min( n-1, ku ) kb = klm + kun kb1 = kb + 1_${ik}$ inca = kb1*ldab nr = 0_${ik}$ j1 = klm + 2_${ik}$ j2 = 1_${ik}$ - kun loop_90: do i = 1, minmn ! reduce i-th column and i-th row of matrix to bidiagonal form ml = klm + 1_${ik}$ mu = kun + 1_${ik}$ loop_80: do kk = 1, kb j1 = j1 + kb j2 = j2 + kb ! generate plane rotations to annihilate nonzero elements ! which have been created below the band if( nr>0_${ik}$ )call stdlib${ii}$_slargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, & work( mn+j1 ), kb1 ) ! apply plane rotations from the left do l = 1, kb if( j2-klm+l-1>n ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( & klu1-l+1, j1-klm+l-1 ), inca,work( mn+j1 ), work( j1 ), kb1 ) end do if( ml>ml0 ) then if( ml<=m-i+1 ) then ! generate plane rotation to annihilate a(i+ml-1,i) ! within the band, and apply rotation from the left call stdlib${ii}$_slartg( ab( ku+ml-1, i ), ab( ku+ml, i ),work( mn+i+ml-1 ), & work( i+ml-1 ),ra ) ab( ku+ml-1, i ) = ra if( i<n )call stdlib${ii}$_srot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-& 1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,work( mn+i+ml-1 ), work( i+ml-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantq ) then ! accumulate product of plane rotations in q do j = j1, j2, kb1 call stdlib${ii}$_srot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,work( mn+j ), work( j & ) ) end do end if if( wantc ) then ! apply plane rotations to c do j = j1, j2, kb1 call stdlib${ii}$_srot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,work( mn+j ), & work( j ) ) end do end if if( j2+kun>n ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j-1,j+ku) above the band ! and store it in work(n+1:2*n) work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun ) ab( 1_${ik}$, j+kun ) = work( mn+j )*ab( 1_${ik}$, j+kun ) end do ! generate plane rotations to annihilate nonzero elements ! which have been generated above the band if( nr>0_${ik}$ )call stdlib${ii}$_slargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,& work( mn+j1+kun ),kb1 ) ! apply plane rotations from the right do l = 1, kb if( j2+l-1>m ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_slartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+& kun ), inca,work( mn+j1+kun ), work( j1+kun ),kb1 ) end do if( ml==ml0 .and. mu>mu0 ) then if( mu<=n-i+1 ) then ! generate plane rotation to annihilate a(i,i+mu-1) ! within the band, and apply rotation from the right call stdlib${ii}$_slartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),work( & mn+i+mu-1 ), work( i+mu-1 ),ra ) ab( ku-mu+3, i+mu-2 ) = ra call stdlib${ii}$_srot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-& mu+3, i+mu-1 ), 1_${ik}$,work( mn+i+mu-1 ), work( i+mu-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantpt ) then ! accumulate product of plane rotations in p**t do j = j1, j2, kb1 call stdlib${ii}$_srot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, work( & mn+j+kun ),work( j+kun ) ) end do end if if( j2+kb>m ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j+kl+ku,j+ku-1) below the ! band and store it in work(1:n) work( j+kb ) = work( j+kun )*ab( klu1, j+kun ) ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun ) end do if( ml>ml0 ) then ml = ml - 1_${ik}$ else mu = mu - 1_${ik}$ end if end do loop_80 end do loop_90 end if if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then ! a has been reduced to lower bidiagonal form ! transform lower bidiagonal form to upper bidiagonal by applying ! plane rotations from the left, storing diagonal elements in d ! and off-diagonal elements in e do i = 1, min( m-1, n ) call stdlib${ii}$_slartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra ) d( i ) = ra if( i<n ) then e( i ) = rs*ab( 1_${ik}$, i+1 ) ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 ) end if if( wantq )call stdlib${ii}$_srot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc, rs ) if( wantc )call stdlib${ii}$_srot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs ) end do if( m<=n )d( m ) = ab( 1_${ik}$, m ) else if( ku>0_${ik}$ ) then ! a has been reduced to upper bidiagonal form if( m<n ) then ! annihilate a(m,m+1) by applying plane rotations from the ! right, storing diagonal elements in d and off-diagonal ! elements in e rb = ab( ku, m+1 ) do i = m, 1, -1 call stdlib${ii}$_slartg( ab( ku+1, i ), rb, rc, rs, ra ) d( i ) = ra if( i>1_${ik}$ ) then rb = -rs*ab( ku, i ) e( i-1 ) = rc*ab( ku, i ) end if if( wantpt )call stdlib${ii}$_srot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, rs ) end do else ! copy off-diagonal elements to e and diagonal elements to d do i = 1, minmn - 1 e( i ) = ab( ku, i+1 ) end do do i = 1, minmn d( i ) = ab( ku+1, i ) end do end if else ! a is diagonal. set elements of e to zero and copy diagonal ! elements to d. do i = 1, minmn - 1 e( i ) = zero end do do i = 1, minmn d( i ) = ab( 1_${ik}$, i ) end do end if return end subroutine stdlib${ii}$_sgbbrd pure module subroutine stdlib${ii}$_dgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, & !! DGBBRD reduces a real general m-by-n band matrix A to upper !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! The routine computes B, and optionally forms Q or P**T, or computes !! Q**T*C for a given matrix C. ldc, work, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc ! Array Arguments real(dp), intent(inout) :: ab(ldab,*), c(ldc,*) real(dp), intent(out) :: d(*), e(*), pt(ldpt,*), q(ldq,*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: wantb, wantc, wantpt, wantq integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mn,& mu, mu0, nr, nrt real(dp) :: ra, rb, rc, rs ! Intrinsic Functions ! Executable Statements ! test the input parameters wantb = stdlib_lsame( vect, 'B' ) wantq = stdlib_lsame( vect, 'Q' ) .or. wantb wantpt = stdlib_lsame( vect, 'P' ) .or. wantb wantc = ncc>0_${ik}$ klu1 = kl + ku + 1_${ik}$ info = 0_${ik}$ if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then info = -1_${ik}$ else if( m<0_${ik}$ ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ncc<0_${ik}$ ) then info = -4_${ik}$ else if( kl<0_${ik}$ ) then info = -5_${ik}$ else if( ku<0_${ik}$ ) then info = -6_${ik}$ else if( ldab<klu1 ) then info = -8_${ik}$ else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then info = -12_${ik}$ else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then info = -14_${ik}$ else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGBBRD', -info ) return end if ! initialize q and p**t to the unit matrix, if needed if( wantq )call stdlib${ii}$_dlaset( 'FULL', m, m, zero, one, q, ldq ) if( wantpt )call stdlib${ii}$_dlaset( 'FULL', n, n, zero, one, pt, ldpt ) ! quick return if possible. if( m==0 .or. n==0 )return minmn = min( m, n ) if( kl+ku>1_${ik}$ ) then ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce ! first to lower bidiagonal form and then transform to upper ! bidiagonal if( ku>0_${ik}$ ) then ml0 = 1_${ik}$ mu0 = 2_${ik}$ else ml0 = 2_${ik}$ mu0 = 1_${ik}$ end if ! wherever possible, plane rotations are generated and applied in ! vector operations of length nr over the index set j1:j2:klu1. ! the sines of the plane rotations are stored in work(1:max(m,n)) ! and the cosines in work(max(m,n)+1:2*max(m,n)). mn = max( m, n ) klm = min( m-1, kl ) kun = min( n-1, ku ) kb = klm + kun kb1 = kb + 1_${ik}$ inca = kb1*ldab nr = 0_${ik}$ j1 = klm + 2_${ik}$ j2 = 1_${ik}$ - kun loop_90: do i = 1, minmn ! reduce i-th column and i-th row of matrix to bidiagonal form ml = klm + 1_${ik}$ mu = kun + 1_${ik}$ loop_80: do kk = 1, kb j1 = j1 + kb j2 = j2 + kb ! generate plane rotations to annihilate nonzero elements ! which have been created below the band if( nr>0_${ik}$ )call stdlib${ii}$_dlargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, & work( mn+j1 ), kb1 ) ! apply plane rotations from the left do l = 1, kb if( j2-klm+l-1>n ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( & klu1-l+1, j1-klm+l-1 ), inca,work( mn+j1 ), work( j1 ), kb1 ) end do if( ml>ml0 ) then if( ml<=m-i+1 ) then ! generate plane rotation to annihilate a(i+ml-1,i) ! within the band, and apply rotation from the left call stdlib${ii}$_dlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),work( mn+i+ml-1 ), & work( i+ml-1 ),ra ) ab( ku+ml-1, i ) = ra if( i<n )call stdlib${ii}$_drot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-& 1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,work( mn+i+ml-1 ), work( i+ml-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantq ) then ! accumulate product of plane rotations in q do j = j1, j2, kb1 call stdlib${ii}$_drot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,work( mn+j ), work( j & ) ) end do end if if( wantc ) then ! apply plane rotations to c do j = j1, j2, kb1 call stdlib${ii}$_drot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,work( mn+j ), & work( j ) ) end do end if if( j2+kun>n ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j-1,j+ku) above the band ! and store it in work(n+1:2*n) work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun ) ab( 1_${ik}$, j+kun ) = work( mn+j )*ab( 1_${ik}$, j+kun ) end do ! generate plane rotations to annihilate nonzero elements ! which have been generated above the band if( nr>0_${ik}$ )call stdlib${ii}$_dlargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,& work( mn+j1+kun ),kb1 ) ! apply plane rotations from the right do l = 1, kb if( j2+l-1>m ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_dlartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+& kun ), inca,work( mn+j1+kun ), work( j1+kun ),kb1 ) end do if( ml==ml0 .and. mu>mu0 ) then if( mu<=n-i+1 ) then ! generate plane rotation to annihilate a(i,i+mu-1) ! within the band, and apply rotation from the right call stdlib${ii}$_dlartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),work( & mn+i+mu-1 ), work( i+mu-1 ),ra ) ab( ku-mu+3, i+mu-2 ) = ra call stdlib${ii}$_drot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-& mu+3, i+mu-1 ), 1_${ik}$,work( mn+i+mu-1 ), work( i+mu-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantpt ) then ! accumulate product of plane rotations in p**t do j = j1, j2, kb1 call stdlib${ii}$_drot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, work( & mn+j+kun ),work( j+kun ) ) end do end if if( j2+kb>m ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j+kl+ku,j+ku-1) below the ! band and store it in work(1:n) work( j+kb ) = work( j+kun )*ab( klu1, j+kun ) ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun ) end do if( ml>ml0 ) then ml = ml - 1_${ik}$ else mu = mu - 1_${ik}$ end if end do loop_80 end do loop_90 end if if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then ! a has been reduced to lower bidiagonal form ! transform lower bidiagonal form to upper bidiagonal by applying ! plane rotations from the left, storing diagonal elements in d ! and off-diagonal elements in e do i = 1, min( m-1, n ) call stdlib${ii}$_dlartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra ) d( i ) = ra if( i<n ) then e( i ) = rs*ab( 1_${ik}$, i+1 ) ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 ) end if if( wantq )call stdlib${ii}$_drot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc, rs ) if( wantc )call stdlib${ii}$_drot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs ) end do if( m<=n )d( m ) = ab( 1_${ik}$, m ) else if( ku>0_${ik}$ ) then ! a has been reduced to upper bidiagonal form if( m<n ) then ! annihilate a(m,m+1) by applying plane rotations from the ! right, storing diagonal elements in d and off-diagonal ! elements in e rb = ab( ku, m+1 ) do i = m, 1, -1 call stdlib${ii}$_dlartg( ab( ku+1, i ), rb, rc, rs, ra ) d( i ) = ra if( i>1_${ik}$ ) then rb = -rs*ab( ku, i ) e( i-1 ) = rc*ab( ku, i ) end if if( wantpt )call stdlib${ii}$_drot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, rs ) end do else ! copy off-diagonal elements to e and diagonal elements to d do i = 1, minmn - 1 e( i ) = ab( ku, i+1 ) end do do i = 1, minmn d( i ) = ab( ku+1, i ) end do end if else ! a is diagonal. set elements of e to zero and copy diagonal ! elements to d. do i = 1, minmn - 1 e( i ) = zero end do do i = 1, minmn d( i ) = ab( 1_${ik}$, i ) end do end if return end subroutine stdlib${ii}$_dgbbrd #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$gbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, & !! DGBBRD: reduces a real general m-by-n band matrix A to upper !! bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. !! The routine computes B, and optionally forms Q or P**T, or computes !! Q**T*C for a given matrix C. ldc, work, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc ! Array Arguments real(${rk}$), intent(inout) :: ab(ldab,*), c(ldc,*) real(${rk}$), intent(out) :: d(*), e(*), pt(ldpt,*), q(ldq,*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: wantb, wantc, wantpt, wantq integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mn,& mu, mu0, nr, nrt real(${rk}$) :: ra, rb, rc, rs ! Intrinsic Functions ! Executable Statements ! test the input parameters wantb = stdlib_lsame( vect, 'B' ) wantq = stdlib_lsame( vect, 'Q' ) .or. wantb wantpt = stdlib_lsame( vect, 'P' ) .or. wantb wantc = ncc>0_${ik}$ klu1 = kl + ku + 1_${ik}$ info = 0_${ik}$ if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then info = -1_${ik}$ else if( m<0_${ik}$ ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ncc<0_${ik}$ ) then info = -4_${ik}$ else if( kl<0_${ik}$ ) then info = -5_${ik}$ else if( ku<0_${ik}$ ) then info = -6_${ik}$ else if( ldab<klu1 ) then info = -8_${ik}$ else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then info = -12_${ik}$ else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then info = -14_${ik}$ else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DGBBRD', -info ) return end if ! initialize q and p**t to the unit matrix, if needed if( wantq )call stdlib${ii}$_${ri}$laset( 'FULL', m, m, zero, one, q, ldq ) if( wantpt )call stdlib${ii}$_${ri}$laset( 'FULL', n, n, zero, one, pt, ldpt ) ! quick return if possible. if( m==0 .or. n==0 )return minmn = min( m, n ) if( kl+ku>1_${ik}$ ) then ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce ! first to lower bidiagonal form and then transform to upper ! bidiagonal if( ku>0_${ik}$ ) then ml0 = 1_${ik}$ mu0 = 2_${ik}$ else ml0 = 2_${ik}$ mu0 = 1_${ik}$ end if ! wherever possible, plane rotations are generated and applied in ! vector operations of length nr over the index set j1:j2:klu1. ! the sines of the plane rotations are stored in work(1:max(m,n)) ! and the cosines in work(max(m,n)+1:2*max(m,n)). mn = max( m, n ) klm = min( m-1, kl ) kun = min( n-1, ku ) kb = klm + kun kb1 = kb + 1_${ik}$ inca = kb1*ldab nr = 0_${ik}$ j1 = klm + 2_${ik}$ j2 = 1_${ik}$ - kun loop_90: do i = 1, minmn ! reduce i-th column and i-th row of matrix to bidiagonal form ml = klm + 1_${ik}$ mu = kun + 1_${ik}$ loop_80: do kk = 1, kb j1 = j1 + kb j2 = j2 + kb ! generate plane rotations to annihilate nonzero elements ! which have been created below the band if( nr>0_${ik}$ )call stdlib${ii}$_${ri}$largv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, & work( mn+j1 ), kb1 ) ! apply plane rotations from the left do l = 1, kb if( j2-klm+l-1>n ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_${ri}$lartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( & klu1-l+1, j1-klm+l-1 ), inca,work( mn+j1 ), work( j1 ), kb1 ) end do if( ml>ml0 ) then if( ml<=m-i+1 ) then ! generate plane rotation to annihilate a(i+ml-1,i) ! within the band, and apply rotation from the left call stdlib${ii}$_${ri}$lartg( ab( ku+ml-1, i ), ab( ku+ml, i ),work( mn+i+ml-1 ), & work( i+ml-1 ),ra ) ab( ku+ml-1, i ) = ra if( i<n )call stdlib${ii}$_${ri}$rot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-& 1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,work( mn+i+ml-1 ), work( i+ml-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantq ) then ! accumulate product of plane rotations in q do j = j1, j2, kb1 call stdlib${ii}$_${ri}$rot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,work( mn+j ), work( j & ) ) end do end if if( wantc ) then ! apply plane rotations to c do j = j1, j2, kb1 call stdlib${ii}$_${ri}$rot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,work( mn+j ), & work( j ) ) end do end if if( j2+kun>n ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j-1,j+ku) above the band ! and store it in work(n+1:2*n) work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun ) ab( 1_${ik}$, j+kun ) = work( mn+j )*ab( 1_${ik}$, j+kun ) end do ! generate plane rotations to annihilate nonzero elements ! which have been generated above the band if( nr>0_${ik}$ )call stdlib${ii}$_${ri}$largv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,& work( mn+j1+kun ),kb1 ) ! apply plane rotations from the right do l = 1, kb if( j2+l-1>m ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_${ri}$lartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+& kun ), inca,work( mn+j1+kun ), work( j1+kun ),kb1 ) end do if( ml==ml0 .and. mu>mu0 ) then if( mu<=n-i+1 ) then ! generate plane rotation to annihilate a(i,i+mu-1) ! within the band, and apply rotation from the right call stdlib${ii}$_${ri}$lartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),work( & mn+i+mu-1 ), work( i+mu-1 ),ra ) ab( ku-mu+3, i+mu-2 ) = ra call stdlib${ii}$_${ri}$rot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-& mu+3, i+mu-1 ), 1_${ik}$,work( mn+i+mu-1 ), work( i+mu-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantpt ) then ! accumulate product of plane rotations in p**t do j = j1, j2, kb1 call stdlib${ii}$_${ri}$rot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, work( & mn+j+kun ),work( j+kun ) ) end do end if if( j2+kb>m ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j+kl+ku,j+ku-1) below the ! band and store it in work(1:n) work( j+kb ) = work( j+kun )*ab( klu1, j+kun ) ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun ) end do if( ml>ml0 ) then ml = ml - 1_${ik}$ else mu = mu - 1_${ik}$ end if end do loop_80 end do loop_90 end if if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then ! a has been reduced to lower bidiagonal form ! transform lower bidiagonal form to upper bidiagonal by applying ! plane rotations from the left, storing diagonal elements in d ! and off-diagonal elements in e do i = 1, min( m-1, n ) call stdlib${ii}$_${ri}$lartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra ) d( i ) = ra if( i<n ) then e( i ) = rs*ab( 1_${ik}$, i+1 ) ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 ) end if if( wantq )call stdlib${ii}$_${ri}$rot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc, rs ) if( wantc )call stdlib${ii}$_${ri}$rot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs ) end do if( m<=n )d( m ) = ab( 1_${ik}$, m ) else if( ku>0_${ik}$ ) then ! a has been reduced to upper bidiagonal form if( m<n ) then ! annihilate a(m,m+1) by applying plane rotations from the ! right, storing diagonal elements in d and off-diagonal ! elements in e rb = ab( ku, m+1 ) do i = m, 1, -1 call stdlib${ii}$_${ri}$lartg( ab( ku+1, i ), rb, rc, rs, ra ) d( i ) = ra if( i>1_${ik}$ ) then rb = -rs*ab( ku, i ) e( i-1 ) = rc*ab( ku, i ) end if if( wantpt )call stdlib${ii}$_${ri}$rot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, rs ) end do else ! copy off-diagonal elements to e and diagonal elements to d do i = 1, minmn - 1 e( i ) = ab( ku, i+1 ) end do do i = 1, minmn d( i ) = ab( ku+1, i ) end do end if else ! a is diagonal. set elements of e to zero and copy diagonal ! elements to d. do i = 1, minmn - 1 e( i ) = zero end do do i = 1, minmn d( i ) = ab( 1_${ik}$, i ) end do end if return end subroutine stdlib${ii}$_${ri}$gbbrd #:endif #:endfor pure module subroutine stdlib${ii}$_cgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, & !! CGBBRD reduces a complex general m-by-n band matrix A to real upper !! bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! The routine computes B, and optionally forms Q or P**H, or computes !! Q**H*C for a given matrix C. ldc, work, rwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc ! Array Arguments real(sp), intent(out) :: d(*), e(*), rwork(*) complex(sp), intent(inout) :: ab(ldab,*), c(ldc,*) complex(sp), intent(out) :: pt(ldpt,*), q(ldq,*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: wantb, wantc, wantpt, wantq integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mu,& mu0, nr, nrt real(sp) :: abst, rc complex(sp) :: ra, rb, rs, t ! Intrinsic Functions ! Executable Statements ! test the input parameters wantb = stdlib_lsame( vect, 'B' ) wantq = stdlib_lsame( vect, 'Q' ) .or. wantb wantpt = stdlib_lsame( vect, 'P' ) .or. wantb wantc = ncc>0_${ik}$ klu1 = kl + ku + 1_${ik}$ info = 0_${ik}$ if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then info = -1_${ik}$ else if( m<0_${ik}$ ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ncc<0_${ik}$ ) then info = -4_${ik}$ else if( kl<0_${ik}$ ) then info = -5_${ik}$ else if( ku<0_${ik}$ ) then info = -6_${ik}$ else if( ldab<klu1 ) then info = -8_${ik}$ else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then info = -12_${ik}$ else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then info = -14_${ik}$ else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CGBBRD', -info ) return end if ! initialize q and p**h to the unit matrix, if needed if( wantq )call stdlib${ii}$_claset( 'FULL', m, m, czero, cone, q, ldq ) if( wantpt )call stdlib${ii}$_claset( 'FULL', n, n, czero, cone, pt, ldpt ) ! quick return if possible. if( m==0 .or. n==0 )return minmn = min( m, n ) if( kl+ku>1_${ik}$ ) then ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce ! first to lower bidiagonal form and then transform to upper ! bidiagonal if( ku>0_${ik}$ ) then ml0 = 1_${ik}$ mu0 = 2_${ik}$ else ml0 = 2_${ik}$ mu0 = 1_${ik}$ end if ! wherever possible, plane rotations are generated and applied in ! vector operations of length nr over the index set j1:j2:klu1. ! the complex sines of the plane rotations are stored in work, ! and the real cosines in rwork. klm = min( m-1, kl ) kun = min( n-1, ku ) kb = klm + kun kb1 = kb + 1_${ik}$ inca = kb1*ldab nr = 0_${ik}$ j1 = klm + 2_${ik}$ j2 = 1_${ik}$ - kun loop_90: do i = 1, minmn ! reduce i-th column and i-th row of matrix to bidiagonal form ml = klm + 1_${ik}$ mu = kun + 1_${ik}$ loop_80: do kk = 1, kb j1 = j1 + kb j2 = j2 + kb ! generate plane rotations to annihilate nonzero elements ! which have been created below the band if( nr>0_${ik}$ )call stdlib${ii}$_clargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, & rwork( j1 ), kb1 ) ! apply plane rotations from the left do l = 1, kb if( j2-klm+l-1>n ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_clartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( & klu1-l+1, j1-klm+l-1 ), inca,rwork( j1 ), work( j1 ), kb1 ) end do if( ml>ml0 ) then if( ml<=m-i+1 ) then ! generate plane rotation to annihilate a(i+ml-1,i) ! within the band, and apply rotation from the left call stdlib${ii}$_clartg( ab( ku+ml-1, i ), ab( ku+ml, i ),rwork( i+ml-1 ), & work( i+ml-1 ), ra ) ab( ku+ml-1, i ) = ra if( i<n )call stdlib${ii}$_crot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-& 1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,rwork( i+ml-1 ), work( i+ml-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantq ) then ! accumulate product of plane rotations in q do j = j1, j2, kb1 call stdlib${ii}$_crot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,rwork( j ), conjg( & work( j ) ) ) end do end if if( wantc ) then ! apply plane rotations to c do j = j1, j2, kb1 call stdlib${ii}$_crot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,rwork( j ), & work( j ) ) end do end if if( j2+kun>n ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j-1,j+ku) above the band ! and store it in work(n+1:2*n) work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun ) ab( 1_${ik}$, j+kun ) = rwork( j )*ab( 1_${ik}$, j+kun ) end do ! generate plane rotations to annihilate nonzero elements ! which have been generated above the band if( nr>0_${ik}$ )call stdlib${ii}$_clargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,& rwork( j1+kun ),kb1 ) ! apply plane rotations from the right do l = 1, kb if( j2+l-1>m ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_clartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+& kun ), inca,rwork( j1+kun ), work( j1+kun ), kb1 ) end do if( ml==ml0 .and. mu>mu0 ) then if( mu<=n-i+1 ) then ! generate plane rotation to annihilate a(i,i+mu-1) ! within the band, and apply rotation from the right call stdlib${ii}$_clartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),rwork( & i+mu-1 ), work( i+mu-1 ), ra ) ab( ku-mu+3, i+mu-2 ) = ra call stdlib${ii}$_crot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-& mu+3, i+mu-1 ), 1_${ik}$,rwork( i+mu-1 ), work( i+mu-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantpt ) then ! accumulate product of plane rotations in p**h do j = j1, j2, kb1 call stdlib${ii}$_crot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, rwork(& j+kun ),conjg( work( j+kun ) ) ) end do end if if( j2+kb>m ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j+kl+ku,j+ku-1) below the ! band and store it in work(1:n) work( j+kb ) = work( j+kun )*ab( klu1, j+kun ) ab( klu1, j+kun ) = rwork( j+kun )*ab( klu1, j+kun ) end do if( ml>ml0 ) then ml = ml - 1_${ik}$ else mu = mu - 1_${ik}$ end if end do loop_80 end do loop_90 end if if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then ! a has been reduced to complex lower bidiagonal form ! transform lower bidiagonal form to upper bidiagonal by applying ! plane rotations from the left, overwriting superdiagonal ! elements on subdiagonal elements do i = 1, min( m-1, n ) call stdlib${ii}$_clartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra ) ab( 1_${ik}$, i ) = ra if( i<n ) then ab( 2_${ik}$, i ) = rs*ab( 1_${ik}$, i+1 ) ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 ) end if if( wantq )call stdlib${ii}$_crot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc,conjg( rs ) ) if( wantc )call stdlib${ii}$_crot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs ) end do else ! a has been reduced to complex upper bidiagonal form or is ! diagonal if( ku>0_${ik}$ .and. m<n ) then ! annihilate a(m,m+1) by applying plane rotations from the ! right rb = ab( ku, m+1 ) do i = m, 1, -1 call stdlib${ii}$_clartg( ab( ku+1, i ), rb, rc, rs, ra ) ab( ku+1, i ) = ra if( i>1_${ik}$ ) then rb = -conjg( rs )*ab( ku, i ) ab( ku, i ) = rc*ab( ku, i ) end if if( wantpt )call stdlib${ii}$_crot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, & conjg( rs ) ) end do end if end if ! make diagonal and superdiagonal elements real, storing them in d ! and e t = ab( ku+1, 1_${ik}$ ) loop_120: do i = 1, minmn abst = abs( t ) d( i ) = abst if( abst/=zero ) then t = t / abst else t = cone end if if( wantq )call stdlib${ii}$_cscal( m, t, q( 1_${ik}$, i ), 1_${ik}$ ) if( wantc )call stdlib${ii}$_cscal( ncc, conjg( t ), c( i, 1_${ik}$ ), ldc ) if( i<minmn ) then if( ku==0_${ik}$ .and. kl==0_${ik}$ ) then e( i ) = zero t = ab( 1_${ik}$, i+1 ) else if( ku==0_${ik}$ ) then t = ab( 2_${ik}$, i )*conjg( t ) else t = ab( ku, i+1 )*conjg( t ) end if abst = abs( t ) e( i ) = abst if( abst/=zero ) then t = t / abst else t = cone end if if( wantpt )call stdlib${ii}$_cscal( n, t, pt( i+1, 1_${ik}$ ), ldpt ) t = ab( ku+1, i+1 )*conjg( t ) end if end if end do loop_120 return end subroutine stdlib${ii}$_cgbbrd pure module subroutine stdlib${ii}$_zgbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, & !! ZGBBRD reduces a complex general m-by-n band matrix A to real upper !! bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! The routine computes B, and optionally forms Q or P**H, or computes !! Q**H*C for a given matrix C. ldc, work, rwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc ! Array Arguments real(dp), intent(out) :: d(*), e(*), rwork(*) complex(dp), intent(inout) :: ab(ldab,*), c(ldc,*) complex(dp), intent(out) :: pt(ldpt,*), q(ldq,*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: wantb, wantc, wantpt, wantq integer(${ik}$) :: i, inca, j, j1, j2, kb, kb1, kk, klm, klu1, kun, l, minmn, ml, ml0, mu,& mu0, nr, nrt real(dp) :: abst, rc complex(dp) :: ra, rb, rs, t ! Intrinsic Functions ! Executable Statements ! test the input parameters wantb = stdlib_lsame( vect, 'B' ) wantq = stdlib_lsame( vect, 'Q' ) .or. wantb wantpt = stdlib_lsame( vect, 'P' ) .or. wantb wantc = ncc>0_${ik}$ klu1 = kl + ku + 1_${ik}$ info = 0_${ik}$ if( .not.wantq .and. .not.wantpt .and. .not.stdlib_lsame( vect, 'N' ) )then info = -1_${ik}$ else if( m<0_${ik}$ ) then info = -2_${ik}$ else if( n<0_${ik}$ ) then info = -3_${ik}$ else if( ncc<0_${ik}$ ) then info = -4_${ik}$ else if( kl<0_${ik}$ ) then info = -5_${ik}$ else if( ku<0_${ik}$ ) then info = -6_${ik}$ else if( ldab<klu1 ) then info = -8_${ik}$ else if( ldq<1_${ik}$ .or. wantq .and. ldq<max( 1_${ik}$, m ) ) then info = -12_${ik}$ else if( ldpt<1_${ik}$ .or. wantpt .and. ldpt<max( 1_${ik}$, n ) ) then info = -14_${ik}$ else if( ldc<1_${ik}$ .or. wantc .and. ldc<max( 1_${ik}$, m ) ) then info = -16_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZGBBRD', -info ) return end if ! initialize q and p**h to the unit matrix, if needed if( wantq )call stdlib${ii}$_zlaset( 'FULL', m, m, czero, cone, q, ldq ) if( wantpt )call stdlib${ii}$_zlaset( 'FULL', n, n, czero, cone, pt, ldpt ) ! quick return if possible. if( m==0 .or. n==0 )return minmn = min( m, n ) if( kl+ku>1_${ik}$ ) then ! reduce to upper bidiagonal form if ku > 0; if ku = 0, reduce ! first to lower bidiagonal form and then transform to upper ! bidiagonal if( ku>0_${ik}$ ) then ml0 = 1_${ik}$ mu0 = 2_${ik}$ else ml0 = 2_${ik}$ mu0 = 1_${ik}$ end if ! wherever possible, plane rotations are generated and applied in ! vector operations of length nr over the index set j1:j2:klu1. ! the complex sines of the plane rotations are stored in work, ! and the real cosines in rwork. klm = min( m-1, kl ) kun = min( n-1, ku ) kb = klm + kun kb1 = kb + 1_${ik}$ inca = kb1*ldab nr = 0_${ik}$ j1 = klm + 2_${ik}$ j2 = 1_${ik}$ - kun loop_90: do i = 1, minmn ! reduce i-th column and i-th row of matrix to bidiagonal form ml = klm + 1_${ik}$ mu = kun + 1_${ik}$ loop_80: do kk = 1, kb j1 = j1 + kb j2 = j2 + kb ! generate plane rotations to annihilate nonzero elements ! which have been created below the band if( nr>0_${ik}$ )call stdlib${ii}$_zlargv( nr, ab( klu1, j1-klm-1 ), inca,work( j1 ), kb1, & rwork( j1 ), kb1 ) ! apply plane rotations from the left do l = 1, kb if( j2-klm+l-1>n ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_zlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,ab( & klu1-l+1, j1-klm+l-1 ), inca,rwork( j1 ), work( j1 ), kb1 ) end do if( ml>ml0 ) then if( ml<=m-i+1 ) then ! generate plane rotation to annihilate a(i+ml-1,i) ! within the band, and apply rotation from the left call stdlib${ii}$_zlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),rwork( i+ml-1 ), & work( i+ml-1 ), ra ) ab( ku+ml-1, i ) = ra if( i<n )call stdlib${ii}$_zrot( min( ku+ml-2, n-i ),ab( ku+ml-2, i+1 ), ldab-& 1_${ik}$,ab( ku+ml-1, i+1 ), ldab-1,rwork( i+ml-1 ), work( i+ml-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantq ) then ! accumulate product of plane rotations in q do j = j1, j2, kb1 call stdlib${ii}$_zrot( m, q( 1_${ik}$, j-1 ), 1_${ik}$, q( 1_${ik}$, j ), 1_${ik}$,rwork( j ), conjg( & work( j ) ) ) end do end if if( wantc ) then ! apply plane rotations to c do j = j1, j2, kb1 call stdlib${ii}$_zrot( ncc, c( j-1, 1_${ik}$ ), ldc, c( j, 1_${ik}$ ), ldc,rwork( j ), & work( j ) ) end do end if if( j2+kun>n ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j-1,j+ku) above the band ! and store it in work(n+1:2*n) work( j+kun ) = work( j )*ab( 1_${ik}$, j+kun ) ab( 1_${ik}$, j+kun ) = rwork( j )*ab( 1_${ik}$, j+kun ) end do ! generate plane rotations to annihilate nonzero elements ! which have been generated above the band if( nr>0_${ik}$ )call stdlib${ii}$_zlargv( nr, ab( 1_${ik}$, j1+kun-1 ), inca,work( j1+kun ), kb1,& rwork( j1+kun ),kb1 ) ! apply plane rotations from the right do l = 1, kb if( j2+l-1>m ) then nrt = nr - 1_${ik}$ else nrt = nr end if if( nrt>0_${ik}$ )call stdlib${ii}$_zlartv( nrt, ab( l+1, j1+kun-1 ), inca,ab( l, j1+& kun ), inca,rwork( j1+kun ), work( j1+kun ), kb1 ) end do if( ml==ml0 .and. mu>mu0 ) then if( mu<=n-i+1 ) then ! generate plane rotation to annihilate a(i,i+mu-1) ! within the band, and apply rotation from the right call stdlib${ii}$_zlartg( ab( ku-mu+3, i+mu-2 ),ab( ku-mu+2, i+mu-1 ),rwork( & i+mu-1 ), work( i+mu-1 ), ra ) ab( ku-mu+3, i+mu-2 ) = ra call stdlib${ii}$_zrot( min( kl+mu-2, m-i ),ab( ku-mu+4, i+mu-2 ), 1_${ik}$,ab( ku-& mu+3, i+mu-1 ), 1_${ik}$,rwork( i+mu-1 ), work( i+mu-1 ) ) end if nr = nr + 1_${ik}$ j1 = j1 - kb1 end if if( wantpt ) then ! accumulate product of plane rotations in p**h do j = j1, j2, kb1 call stdlib${ii}$_zrot( n, pt( j+kun-1, 1_${ik}$ ), ldpt,pt( j+kun, 1_${ik}$ ), ldpt, rwork(& j+kun ),conjg( work( j+kun ) ) ) end do end if if( j2+kb>m ) then ! adjust j2 to keep within the bounds of the matrix nr = nr - 1_${ik}$ j2 = j2 - kb1 end if do j = j1, j2, kb1 ! create nonzero element a(j+kl+ku,j+ku-1) below the ! band and store it in work(1:n) work( j+kb ) = work( j+kun )*ab( klu1, j+kun ) ab( klu1, j+kun ) = rwork( j+kun )*ab( klu1, j+kun ) end do if( ml>ml0 ) then ml = ml - 1_${ik}$ else mu = mu - 1_${ik}$ end if end do loop_80 end do loop_90 end if if( ku==0_${ik}$ .and. kl>0_${ik}$ ) then ! a has been reduced to complex lower bidiagonal form ! transform lower bidiagonal form to upper bidiagonal by applying ! plane rotations from the left, overwriting superdiagonal ! elements on subdiagonal elements do i = 1, min( m-1, n ) call stdlib${ii}$_zlartg( ab( 1_${ik}$, i ), ab( 2_${ik}$, i ), rc, rs, ra ) ab( 1_${ik}$, i ) = ra if( i<n ) then ab( 2_${ik}$, i ) = rs*ab( 1_${ik}$, i+1 ) ab( 1_${ik}$, i+1 ) = rc*ab( 1_${ik}$, i+1 ) end if if( wantq )call stdlib${ii}$_zrot( m, q( 1_${ik}$, i ), 1_${ik}$, q( 1_${ik}$, i+1 ), 1_${ik}$, rc,conjg( rs ) ) if( wantc )call stdlib${ii}$_zrot( ncc, c( i, 1_${ik}$ ), ldc, c( i+1, 1_${ik}$ ), ldc, rc,rs ) end do else ! a has been reduced to complex upper bidiagonal form or is ! diagonal if( ku>0_${ik}$ .and. m<n ) then ! annihilate a(m,m+1) by applying plane rotations from the ! right rb = ab( ku, m+1 ) do i = m, 1, -1 call stdlib${ii}$_zlartg( ab( ku+1, i ), rb, rc, rs, ra ) ab( ku+1, i ) = ra if( i>1_${ik}$ ) then rb = -conjg( rs )*ab( ku, i ) ab( ku, i ) = rc*ab( ku, i ) end if if( wantpt )call stdlib${ii}$_zrot( n, pt( i, 1_${ik}$ ), ldpt, pt( m+1, 1_${ik}$ ), ldpt,rc, & conjg( rs ) ) end do end if end if ! make diagonal and superdiagonal elements real, storing them in d ! and e t = ab( ku+1, 1_${ik}$ ) loop_120: do i = 1, minmn abst = abs( t ) d( i ) = abst if( abst/=zero ) then t = t / abst else t = cone end if if( wantq )call stdlib${ii}$_zscal( m, t, q( 1_${ik}$, i ), 1_${ik}$ ) if( wantc )call stdlib${ii}$_zscal( ncc, conjg( t ), c( i, 1_${ik}$ ), ldc ) if( i<minmn ) then if( ku==0_${ik}$ .and. kl==0_${ik}$ ) then e( i ) = zero t = ab( 1_${ik}$, i+1 ) else if( ku==0_${ik}$ ) then t = ab( 2_${ik}$, i )*conjg( t ) else t = ab( ku, i+1 )*conjg( t ) end if abst = abs( t ) e( i ) = abst if( abst/=zero ) then t = t / abst else t = cone end if if( wantpt )call stdlib${ii}$_zscal( n, t, pt( i+1, 1_${ik}$ ), ldpt ) t = ab( ku+1, i+1 )*conjg( t ) end if end if end do loop_120 return end subroutine stdlib${ii}$_zgbbrd #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$gbbrd( vect, m, n, ncc, kl, ku, ab, ldab, d, e, q,ldq, pt, ldpt, c, & !! ZGBBRD: reduces a complex general m-by-n band matrix A to real upper !! bidiagonal form B by a unitary transformation: Q**H * A * P = B. !! The routine computes B, and optionally forms Q or P**H, or computes !! Q**H*C for a given matrix C. ldc, work, rwork, info ) ! -- lapack computational routine -- ! -- lapack is a software package provided by univ. of tennessee, -- ! -- univ. of california berkeley, univ. of colorado denver and nag ltd..-- use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone ! Scalar Arguments character, intent(in) :: vect integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc ! Array Arguments real(${ck}$), intent(out) :: d(*), e(*), rwork(*) complex(${ck}$), intent(inout) :: ab(ldab,*), c(ldc,*) complex(${ck}$), intent(out) :: pt(ldpt,*), q(ldq,*), work(*) ! ===================================================================== ! Local Scalars logical(lk) :: wantb