#: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