#:include "common.fypp"
submodule(stdlib_lapack_eig_svd_lsq) stdlib_lapack_lsq_constrained
implicit none
contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
pure module subroutine stdlib${ii}$_sgglse( m, n, p, a, lda, b, ldb, c, d, x, work, lwork,info )
!! SGGLSE solves the linear equality-constrained least squares (LSE)
!! problem:
!! minimize || c - A*x ||_2 subject to B*x = d
!! where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!! M-vector, and d is a given P-vector. It is assumed that
!! P <= N <= M+P, and
!! rank(B) = P and rank( (A) ) = N.
!! ( (B) )
!! These conditions ensure that the LSE problem has a unique solution,
!! which is obtained using a generalized RQ factorization of the
!! matrices (B, A) given by
!! B = (0 R)*Q, A = Z*T*Q.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
real(sp), intent(inout) :: a(lda,*), b(ldb,*), c(*), d(*)
real(sp), intent(out) :: work(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: lopt, lwkmin, lwkopt, mn, nb, nb1, nb2, nb3, nb4, nr
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
mn = min( m, n )
lquery = ( lwork==-1_${ik}$ )
if( m<0_${ik}$ ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p>n .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, m ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, p ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SGEQRF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SGERQF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SORMQR', ' ', m, n, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SORMRQ', ' ', m, n, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = p + mn + max( m, n )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGGLSE', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! compute the grq factorization of matrices b and a:
! b*q**t = ( 0 t12 ) p z**t*a*q**t = ( r11 r12 ) n-p
! n-p p ( 0 r22 ) m+p-n
! n-p p
! where t12 and r11 are upper triangular, and q and z are
! orthogonal.
call stdlib${ii}$_sggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),work( p+mn+1 ), lwork-p-&
mn, info )
lopt = work( p+mn+1 )
! update c = z**t *c = ( c1 ) n-p
! ( c2 ) m+p-n
call stdlib${ii}$_sormqr( 'LEFT', 'TRANSPOSE', m, 1_${ik}$, mn, a, lda, work( p+1 ),c, max( 1_${ik}$, m ), &
work( p+mn+1 ), lwork-p-mn, info )
lopt = max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
! solve t12*x2 = d for x2
if( p>0_${ik}$ ) then
call stdlib${ii}$_strtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', p, 1_${ik}$,b( 1_${ik}$, n-p+1 ), ldb, d,&
p, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
! put the solution in x
call stdlib${ii}$_scopy( p, d, 1_${ik}$, x( n-p+1 ), 1_${ik}$ )
! update c1
call stdlib${ii}$_sgemv( 'NO TRANSPOSE', n-p, p, -one, a( 1_${ik}$, n-p+1 ), lda,d, 1_${ik}$, one, c, 1_${ik}$ &
)
end if
! solve r11*x1 = c1 for x1
if( n>p ) then
call stdlib${ii}$_strtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n-p, 1_${ik}$,a, lda, c, n-p, &
info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! put the solutions in x
call stdlib${ii}$_scopy( n-p, c, 1_${ik}$, x, 1_${ik}$ )
end if
! compute the residual vector:
if( m<n ) then
nr = m + p - n
if( nr>0_${ik}$ )call stdlib${ii}$_sgemv( 'NO TRANSPOSE', nr, n-m, -one, a( n-p+1, m+1 ),lda, d( &
nr+1 ), 1_${ik}$, one, c( n-p+1 ), 1_${ik}$ )
else
nr = p
end if
if( nr>0_${ik}$ ) then
call stdlib${ii}$_strmv( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', nr,a( n-p+1, n-p+1 ), lda, &
d, 1_${ik}$ )
call stdlib${ii}$_saxpy( nr, -one, d, 1_${ik}$, c( n-p+1 ), 1_${ik}$ )
end if
! backward transformation x = q**t*x
call stdlib${ii}$_sormrq( 'LEFT', 'TRANSPOSE', n, 1_${ik}$, p, b, ldb, work( 1_${ik}$ ), x,n, work( p+mn+1 &
), lwork-p-mn, info )
work( 1_${ik}$ ) = p + mn + max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_sgglse
pure module subroutine stdlib${ii}$_dgglse( m, n, p, a, lda, b, ldb, c, d, x, work, lwork,info )
!! DGGLSE solves the linear equality-constrained least squares (LSE)
!! problem:
!! minimize || c - A*x ||_2 subject to B*x = d
!! where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!! M-vector, and d is a given P-vector. It is assumed that
!! P <= N <= M+P, and
!! rank(B) = P and rank( (A) ) = N.
!! ( (B) )
!! These conditions ensure that the LSE problem has a unique solution,
!! which is obtained using a generalized RQ factorization of the
!! matrices (B, A) given by
!! B = (0 R)*Q, A = Z*T*Q.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
real(dp), intent(inout) :: a(lda,*), b(ldb,*), c(*), d(*)
real(dp), intent(out) :: work(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: lopt, lwkmin, lwkopt, mn, nb, nb1, nb2, nb3, nb4, nr
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
mn = min( m, n )
lquery = ( lwork==-1_${ik}$ )
if( m<0_${ik}$ ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p>n .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, m ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, p ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEQRF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGERQF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMQR', ' ', m, n, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMRQ', ' ', m, n, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = p + mn + max( m, n )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGGLSE', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! compute the grq factorization of matrices b and a:
! b*q**t = ( 0 t12 ) p z**t*a*q**t = ( r11 r12 ) n-p
! n-p p ( 0 r22 ) m+p-n
! n-p p
! where t12 and r11 are upper triangular, and q and z are
! orthogonal.
call stdlib${ii}$_dggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),work( p+mn+1 ), lwork-p-&
mn, info )
lopt = work( p+mn+1 )
! update c = z**t *c = ( c1 ) n-p
! ( c2 ) m+p-n
call stdlib${ii}$_dormqr( 'LEFT', 'TRANSPOSE', m, 1_${ik}$, mn, a, lda, work( p+1 ),c, max( 1_${ik}$, m ), &
work( p+mn+1 ), lwork-p-mn, info )
lopt = max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
! solve t12*x2 = d for x2
if( p>0_${ik}$ ) then
call stdlib${ii}$_dtrtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', p, 1_${ik}$,b( 1_${ik}$, n-p+1 ), ldb, d,&
p, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
! put the solution in x
call stdlib${ii}$_dcopy( p, d, 1_${ik}$, x( n-p+1 ), 1_${ik}$ )
! update c1
call stdlib${ii}$_dgemv( 'NO TRANSPOSE', n-p, p, -one, a( 1_${ik}$, n-p+1 ), lda,d, 1_${ik}$, one, c, 1_${ik}$ &
)
end if
! solve r11*x1 = c1 for x1
if( n>p ) then
call stdlib${ii}$_dtrtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n-p, 1_${ik}$,a, lda, c, n-p, &
info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! put the solutions in x
call stdlib${ii}$_dcopy( n-p, c, 1_${ik}$, x, 1_${ik}$ )
end if
! compute the residual vector:
if( m<n ) then
nr = m + p - n
if( nr>0_${ik}$ )call stdlib${ii}$_dgemv( 'NO TRANSPOSE', nr, n-m, -one, a( n-p+1, m+1 ),lda, d( &
nr+1 ), 1_${ik}$, one, c( n-p+1 ), 1_${ik}$ )
else
nr = p
end if
if( nr>0_${ik}$ ) then
call stdlib${ii}$_dtrmv( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', nr,a( n-p+1, n-p+1 ), lda, &
d, 1_${ik}$ )
call stdlib${ii}$_daxpy( nr, -one, d, 1_${ik}$, c( n-p+1 ), 1_${ik}$ )
end if
! backward transformation x = q**t*x
call stdlib${ii}$_dormrq( 'LEFT', 'TRANSPOSE', n, 1_${ik}$, p, b, ldb, work( 1_${ik}$ ), x,n, work( p+mn+1 &
), lwork-p-mn, info )
work( 1_${ik}$ ) = p + mn + max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_dgglse
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$gglse( m, n, p, a, lda, b, ldb, c, d, x, work, lwork,info )
!! DGGLSE: solves the linear equality-constrained least squares (LSE)
!! problem:
!! minimize || c - A*x ||_2 subject to B*x = d
!! where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!! M-vector, and d is a given P-vector. It is assumed that
!! P <= N <= M+P, and
!! rank(B) = P and rank( (A) ) = N.
!! ( (B) )
!! These conditions ensure that the LSE problem has a unique solution,
!! which is obtained using a generalized RQ factorization of the
!! matrices (B, A) given by
!! B = (0 R)*Q, A = Z*T*Q.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*), c(*), d(*)
real(${rk}$), intent(out) :: work(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: lopt, lwkmin, lwkopt, mn, nb, nb1, nb2, nb3, nb4, nr
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
mn = min( m, n )
lquery = ( lwork==-1_${ik}$ )
if( m<0_${ik}$ ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p>n .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, m ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, p ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEQRF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGERQF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMQR', ' ', m, n, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMRQ', ' ', m, n, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = p + mn + max( m, n )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGGLSE', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! compute the grq factorization of matrices b and a:
! b*q**t = ( 0 t12 ) p z**t*a*q**t = ( r11 r12 ) n-p
! n-p p ( 0 r22 ) m+p-n
! n-p p
! where t12 and r11 are upper triangular, and q and z are
! orthogonal.
call stdlib${ii}$_${ri}$ggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),work( p+mn+1 ), lwork-p-&
mn, info )
lopt = work( p+mn+1 )
! update c = z**t *c = ( c1 ) n-p
! ( c2 ) m+p-n
call stdlib${ii}$_${ri}$ormqr( 'LEFT', 'TRANSPOSE', m, 1_${ik}$, mn, a, lda, work( p+1 ),c, max( 1_${ik}$, m ), &
work( p+mn+1 ), lwork-p-mn, info )
lopt = max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
! solve t12*x2 = d for x2
if( p>0_${ik}$ ) then
call stdlib${ii}$_${ri}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', p, 1_${ik}$,b( 1_${ik}$, n-p+1 ), ldb, d,&
p, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
! put the solution in x
call stdlib${ii}$_${ri}$copy( p, d, 1_${ik}$, x( n-p+1 ), 1_${ik}$ )
! update c1
call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', n-p, p, -one, a( 1_${ik}$, n-p+1 ), lda,d, 1_${ik}$, one, c, 1_${ik}$ &
)
end if
! solve r11*x1 = c1 for x1
if( n>p ) then
call stdlib${ii}$_${ri}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n-p, 1_${ik}$,a, lda, c, n-p, &
info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! put the solutions in x
call stdlib${ii}$_${ri}$copy( n-p, c, 1_${ik}$, x, 1_${ik}$ )
end if
! compute the residual vector:
if( m<n ) then
nr = m + p - n
if( nr>0_${ik}$ )call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', nr, n-m, -one, a( n-p+1, m+1 ),lda, d( &
nr+1 ), 1_${ik}$, one, c( n-p+1 ), 1_${ik}$ )
else
nr = p
end if
if( nr>0_${ik}$ ) then
call stdlib${ii}$_${ri}$trmv( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', nr,a( n-p+1, n-p+1 ), lda, &
d, 1_${ik}$ )
call stdlib${ii}$_${ri}$axpy( nr, -one, d, 1_${ik}$, c( n-p+1 ), 1_${ik}$ )
end if
! backward transformation x = q**t*x
call stdlib${ii}$_${ri}$ormrq( 'LEFT', 'TRANSPOSE', n, 1_${ik}$, p, b, ldb, work( 1_${ik}$ ), x,n, work( p+mn+1 &
), lwork-p-mn, info )
work( 1_${ik}$ ) = p + mn + max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_${ri}$gglse
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cgglse( m, n, p, a, lda, b, ldb, c, d, x, work, lwork,info )
!! CGGLSE solves the linear equality-constrained least squares (LSE)
!! problem:
!! minimize || c - A*x ||_2 subject to B*x = d
!! where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!! M-vector, and d is a given P-vector. It is assumed that
!! P <= N <= M+P, and
!! rank(B) = P and rank( (A) ) = N.
!! ( (B) )
!! These conditions ensure that the LSE problem has a unique solution,
!! which is obtained using a generalized RQ factorization of the
!! matrices (B, A) given by
!! B = (0 R)*Q, A = Z*T*Q.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
complex(sp), intent(inout) :: a(lda,*), b(ldb,*), c(*), d(*)
complex(sp), intent(out) :: work(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: lopt, lwkmin, lwkopt, mn, nb, nb1, nb2, nb3, nb4, nr
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
mn = min( m, n )
lquery = ( lwork==-1_${ik}$ )
if( m<0_${ik}$ ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p>n .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, m ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, p ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CGEQRF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CGERQF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CUNMQR', ' ', m, n, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CUNMRQ', ' ', m, n, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = p + mn + max( m, n )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGGLSE', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! compute the grq factorization of matrices b and a:
! b*q**h = ( 0 t12 ) p z**h*a*q**h = ( r11 r12 ) n-p
! n-p p ( 0 r22 ) m+p-n
! n-p p
! where t12 and r11 are upper triangular, and q and z are
! unitary.
call stdlib${ii}$_cggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),work( p+mn+1 ), lwork-p-&
mn, info )
lopt = real( work( p+mn+1 ),KIND=sp)
! update c = z**h *c = ( c1 ) n-p
! ( c2 ) m+p-n
call stdlib${ii}$_cunmqr( 'LEFT', 'CONJUGATE TRANSPOSE', m, 1_${ik}$, mn, a, lda,work( p+1 ), c, &
max( 1_${ik}$, m ), work( p+mn+1 ),lwork-p-mn, info )
lopt = max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
! solve t12*x2 = d for x2
if( p>0_${ik}$ ) then
call stdlib${ii}$_ctrtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', p, 1_${ik}$,b( 1_${ik}$, n-p+1 ), ldb, d,&
p, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
! put the solution in x
call stdlib${ii}$_ccopy( p, d, 1_${ik}$, x( n-p+1 ), 1_${ik}$ )
! update c1
call stdlib${ii}$_cgemv( 'NO TRANSPOSE', n-p, p, -cone, a( 1_${ik}$, n-p+1 ), lda,d, 1_${ik}$, cone, c, &
1_${ik}$ )
end if
! solve r11*x1 = c1 for x1
if( n>p ) then
call stdlib${ii}$_ctrtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n-p, 1_${ik}$,a, lda, c, n-p, &
info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! put the solutions in x
call stdlib${ii}$_ccopy( n-p, c, 1_${ik}$, x, 1_${ik}$ )
end if
! compute the residual vector:
if( m<n ) then
nr = m + p - n
if( nr>0_${ik}$ )call stdlib${ii}$_cgemv( 'NO TRANSPOSE', nr, n-m, -cone, a( n-p+1, m+1 ),lda, d(&
nr+1 ), 1_${ik}$, cone, c( n-p+1 ), 1_${ik}$ )
else
nr = p
end if
if( nr>0_${ik}$ ) then
call stdlib${ii}$_ctrmv( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', nr,a( n-p+1, n-p+1 ), lda, &
d, 1_${ik}$ )
call stdlib${ii}$_caxpy( nr, -cone, d, 1_${ik}$, c( n-p+1 ), 1_${ik}$ )
end if
! backward transformation x = q**h*x
call stdlib${ii}$_cunmrq( 'LEFT', 'CONJUGATE TRANSPOSE', n, 1_${ik}$, p, b, ldb,work( 1_${ik}$ ), x, n, &
work( p+mn+1 ), lwork-p-mn, info )
work( 1_${ik}$ ) = p + mn + max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_cgglse
pure module subroutine stdlib${ii}$_zgglse( m, n, p, a, lda, b, ldb, c, d, x, work, lwork,info )
!! ZGGLSE solves the linear equality-constrained least squares (LSE)
!! problem:
!! minimize || c - A*x ||_2 subject to B*x = d
!! where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!! M-vector, and d is a given P-vector. It is assumed that
!! P <= N <= M+P, and
!! rank(B) = P and rank( (A) ) = N.
!! ( (B) )
!! These conditions ensure that the LSE problem has a unique solution,
!! which is obtained using a generalized RQ factorization of the
!! matrices (B, A) given by
!! B = (0 R)*Q, A = Z*T*Q.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
complex(dp), intent(inout) :: a(lda,*), b(ldb,*), c(*), d(*)
complex(dp), intent(out) :: work(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: lopt, lwkmin, lwkopt, mn, nb, nb1, nb2, nb3, nb4, nr
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
mn = min( m, n )
lquery = ( lwork==-1_${ik}$ )
if( m<0_${ik}$ ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p>n .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, m ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, p ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEQRF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGERQF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMQR', ' ', m, n, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMRQ', ' ', m, n, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = p + mn + max( m, n )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGGLSE', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! compute the grq factorization of matrices b and a:
! b*q**h = ( 0 t12 ) p z**h*a*q**h = ( r11 r12 ) n-p
! n-p p ( 0 r22 ) m+p-n
! n-p p
! where t12 and r11 are upper triangular, and q and z are
! unitary.
call stdlib${ii}$_zggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),work( p+mn+1 ), lwork-p-&
mn, info )
lopt = real( work( p+mn+1 ),KIND=dp)
! update c = z**h *c = ( c1 ) n-p
! ( c2 ) m+p-n
call stdlib${ii}$_zunmqr( 'LEFT', 'CONJUGATE TRANSPOSE', m, 1_${ik}$, mn, a, lda,work( p+1 ), c, &
max( 1_${ik}$, m ), work( p+mn+1 ),lwork-p-mn, info )
lopt = max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
! solve t12*x2 = d for x2
if( p>0_${ik}$ ) then
call stdlib${ii}$_ztrtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', p, 1_${ik}$,b( 1_${ik}$, n-p+1 ), ldb, d,&
p, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
! put the solution in x
call stdlib${ii}$_zcopy( p, d, 1_${ik}$, x( n-p+1 ), 1_${ik}$ )
! update c1
call stdlib${ii}$_zgemv( 'NO TRANSPOSE', n-p, p, -cone, a( 1_${ik}$, n-p+1 ), lda,d, 1_${ik}$, cone, c, &
1_${ik}$ )
end if
! solve r11*x1 = c1 for x1
if( n>p ) then
call stdlib${ii}$_ztrtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n-p, 1_${ik}$,a, lda, c, n-p, &
info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! put the solutions in x
call stdlib${ii}$_zcopy( n-p, c, 1_${ik}$, x, 1_${ik}$ )
end if
! compute the residual vector:
if( m<n ) then
nr = m + p - n
if( nr>0_${ik}$ )call stdlib${ii}$_zgemv( 'NO TRANSPOSE', nr, n-m, -cone, a( n-p+1, m+1 ),lda, d(&
nr+1 ), 1_${ik}$, cone, c( n-p+1 ), 1_${ik}$ )
else
nr = p
end if
if( nr>0_${ik}$ ) then
call stdlib${ii}$_ztrmv( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', nr,a( n-p+1, n-p+1 ), lda, &
d, 1_${ik}$ )
call stdlib${ii}$_zaxpy( nr, -cone, d, 1_${ik}$, c( n-p+1 ), 1_${ik}$ )
end if
! backward transformation x = q**h*x
call stdlib${ii}$_zunmrq( 'LEFT', 'CONJUGATE TRANSPOSE', n, 1_${ik}$, p, b, ldb,work( 1_${ik}$ ), x, n, &
work( p+mn+1 ), lwork-p-mn, info )
work( 1_${ik}$ ) = p + mn + max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_zgglse
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$gglse( m, n, p, a, lda, b, ldb, c, d, x, work, lwork,info )
!! ZGGLSE: solves the linear equality-constrained least squares (LSE)
!! problem:
!! minimize || c - A*x ||_2 subject to B*x = d
!! where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
!! M-vector, and d is a given P-vector. It is assumed that
!! P <= N <= M+P, and
!! rank(B) = P and rank( (A) ) = N.
!! ( (B) )
!! These conditions ensure that the LSE problem has a unique solution,
!! which is obtained using a generalized RQ factorization of the
!! matrices (B, A) given by
!! B = (0 R)*Q, A = Z*T*Q.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*), c(*), d(*)
complex(${ck}$), intent(out) :: work(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: lopt, lwkmin, lwkopt, mn, nb, nb1, nb2, nb3, nb4, nr
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
mn = min( m, n )
lquery = ( lwork==-1_${ik}$ )
if( m<0_${ik}$ ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p>n .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, m ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, p ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEQRF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGERQF', ' ', m, n, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMQR', ' ', m, n, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMRQ', ' ', m, n, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = p + mn + max( m, n )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGGLSE', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0 )return
! compute the grq factorization of matrices b and a:
! b*q**h = ( 0 t12 ) p z**h*a*q**h = ( r11 r12 ) n-p
! n-p p ( 0 r22 ) m+p-n
! n-p p
! where t12 and r11 are upper triangular, and q and z are
! unitary.
call stdlib${ii}$_${ci}$ggrqf( p, m, n, b, ldb, work, a, lda, work( p+1 ),work( p+mn+1 ), lwork-p-&
mn, info )
lopt = real( work( p+mn+1 ),KIND=${ck}$)
! update c = z**h *c = ( c1 ) n-p
! ( c2 ) m+p-n
call stdlib${ii}$_${ci}$unmqr( 'LEFT', 'CONJUGATE TRANSPOSE', m, 1_${ik}$, mn, a, lda,work( p+1 ), c, &
max( 1_${ik}$, m ), work( p+mn+1 ),lwork-p-mn, info )
lopt = max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
! solve t12*x2 = d for x2
if( p>0_${ik}$ ) then
call stdlib${ii}$_${ci}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', p, 1_${ik}$,b( 1_${ik}$, n-p+1 ), ldb, d,&
p, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
! put the solution in x
call stdlib${ii}$_${ci}$copy( p, d, 1_${ik}$, x( n-p+1 ), 1_${ik}$ )
! update c1
call stdlib${ii}$_${ci}$gemv( 'NO TRANSPOSE', n-p, p, -cone, a( 1_${ik}$, n-p+1 ), lda,d, 1_${ik}$, cone, c, &
1_${ik}$ )
end if
! solve r11*x1 = c1 for x1
if( n>p ) then
call stdlib${ii}$_${ci}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', n-p, 1_${ik}$,a, lda, c, n-p, &
info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! put the solutions in x
call stdlib${ii}$_${ci}$copy( n-p, c, 1_${ik}$, x, 1_${ik}$ )
end if
! compute the residual vector:
if( m<n ) then
nr = m + p - n
if( nr>0_${ik}$ )call stdlib${ii}$_${ci}$gemv( 'NO TRANSPOSE', nr, n-m, -cone, a( n-p+1, m+1 ),lda, d(&
nr+1 ), 1_${ik}$, cone, c( n-p+1 ), 1_${ik}$ )
else
nr = p
end if
if( nr>0_${ik}$ ) then
call stdlib${ii}$_${ci}$trmv( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', nr,a( n-p+1, n-p+1 ), lda, &
d, 1_${ik}$ )
call stdlib${ii}$_${ci}$axpy( nr, -cone, d, 1_${ik}$, c( n-p+1 ), 1_${ik}$ )
end if
! backward transformation x = q**h*x
call stdlib${ii}$_${ci}$unmrq( 'LEFT', 'CONJUGATE TRANSPOSE', n, 1_${ik}$, p, b, ldb,work( 1_${ik}$ ), x, n, &
work( p+mn+1 ), lwork-p-mn, info )
work( 1_${ik}$ ) = p + mn + max( lopt, int( work( p+mn+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_${ci}$gglse
#:endif
#:endfor
pure module subroutine stdlib${ii}$_sggglm( n, m, p, a, lda, b, ldb, d, x, y, work, lwork,info )
!! SGGGLM solves a general Gauss-Markov linear model (GLM) problem:
!! minimize || y ||_2 subject to d = A*x + B*y
!! x
!! where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!! given N-vector. It is assumed that M <= N <= M+P, and
!! rank(A) = M and rank( A B ) = N.
!! Under these assumptions, the constrained equation is always
!! consistent, and there is a unique solution x and a minimal 2-norm
!! solution y, which is obtained using a generalized QR factorization
!! of the matrices (A, B) given by
!! A = Q*(R), B = Q*T*Z.
!! (0)
!! In particular, if matrix B is square nonsingular, then the problem
!! GLM is equivalent to the following weighted linear least squares
!! problem
!! minimize || inv(B)*(d-A*x) ||_2
!! x
!! where inv(B) denotes the inverse of B.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
real(sp), intent(inout) :: a(lda,*), b(ldb,*), d(*)
real(sp), intent(out) :: work(*), x(*), y(*)
! ===================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: i, lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3, nb4, np
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
np = min( n, p )
lquery = ( lwork==-1_${ik}$ )
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( m<0_${ik}$ .or. m>n ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SGEQRF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SGERQF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SORMQR', ' ', n, m, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'SORMRQ', ' ', n, m, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = m + np + max( n, p )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGGGLM', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
do i = 1, m
x(i) = zero
end do
do i = 1, p
y(i) = zero
end do
return
end if
! compute the gqr factorization of matrices a and b:
! q**t*a = ( r11 ) m, q**t*b*z**t = ( t11 t12 ) m
! ( 0 ) n-m ( 0 t22 ) n-m
! m m+p-n n-m
! where r11 and t22 are upper triangular, and q and z are
! orthogonal.
call stdlib${ii}$_sggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),work( m+np+1 ), lwork-m-&
np, info )
lopt = work( m+np+1 )
! update left-hand-side vector d = q**t*d = ( d1 ) m
! ( d2 ) n-m
call stdlib${ii}$_sormqr( 'LEFT', 'TRANSPOSE', n, 1_${ik}$, m, a, lda, work, d,max( 1_${ik}$, n ), work( m+&
np+1 ), lwork-m-np, info )
lopt = max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
! solve t22*y2 = d2 for y2
if( n>m ) then
call stdlib${ii}$_strtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', n-m, 1_${ik}$,b( m+1, m+p-n+1 ), &
ldb, d( m+1 ), n-m, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
call stdlib${ii}$_scopy( n-m, d( m+1 ), 1_${ik}$, y( m+p-n+1 ), 1_${ik}$ )
end if
! set y1 = 0
do i = 1, m + p - n
y( i ) = zero
end do
! update d1 = d1 - t12*y2
call stdlib${ii}$_sgemv( 'NO TRANSPOSE', m, n-m, -one, b( 1_${ik}$, m+p-n+1 ), ldb,y( m+p-n+1 ), 1_${ik}$, &
one, d, 1_${ik}$ )
! solve triangular system: r11*x = d1
if( m>0_${ik}$ ) then
call stdlib${ii}$_strtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', m, 1_${ik}$, a, lda,d, m, info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! copy d to x
call stdlib${ii}$_scopy( m, d, 1_${ik}$, x, 1_${ik}$ )
end if
! backward transformation y = z**t *y
call stdlib${ii}$_sormrq( 'LEFT', 'TRANSPOSE', p, 1_${ik}$, np,b( max( 1_${ik}$, n-p+1 ), 1_${ik}$ ), ldb, work( &
m+1 ), y,max( 1_${ik}$, p ), work( m+np+1 ), lwork-m-np, info )
work( 1_${ik}$ ) = m + np + max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_sggglm
pure module subroutine stdlib${ii}$_dggglm( n, m, p, a, lda, b, ldb, d, x, y, work, lwork,info )
!! DGGGLM solves a general Gauss-Markov linear model (GLM) problem:
!! minimize || y ||_2 subject to d = A*x + B*y
!! x
!! where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!! given N-vector. It is assumed that M <= N <= M+P, and
!! rank(A) = M and rank( A B ) = N.
!! Under these assumptions, the constrained equation is always
!! consistent, and there is a unique solution x and a minimal 2-norm
!! solution y, which is obtained using a generalized QR factorization
!! of the matrices (A, B) given by
!! A = Q*(R), B = Q*T*Z.
!! (0)
!! In particular, if matrix B is square nonsingular, then the problem
!! GLM is equivalent to the following weighted linear least squares
!! problem
!! minimize || inv(B)*(d-A*x) ||_2
!! x
!! where inv(B) denotes the inverse of B.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
real(dp), intent(inout) :: a(lda,*), b(ldb,*), d(*)
real(dp), intent(out) :: work(*), x(*), y(*)
! ===================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: i, lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3, nb4, np
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
np = min( n, p )
lquery = ( lwork==-1_${ik}$ )
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( m<0_${ik}$ .or. m>n ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEQRF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGERQF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMQR', ' ', n, m, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMRQ', ' ', n, m, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = m + np + max( n, p )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGGGLM', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
do i = 1, m
x(i) = zero
end do
do i = 1, p
y(i) = zero
end do
return
end if
! compute the gqr factorization of matrices a and b:
! q**t*a = ( r11 ) m, q**t*b*z**t = ( t11 t12 ) m
! ( 0 ) n-m ( 0 t22 ) n-m
! m m+p-n n-m
! where r11 and t22 are upper triangular, and q and z are
! orthogonal.
call stdlib${ii}$_dggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),work( m+np+1 ), lwork-m-&
np, info )
lopt = work( m+np+1 )
! update left-hand-side vector d = q**t*d = ( d1 ) m
! ( d2 ) n-m
call stdlib${ii}$_dormqr( 'LEFT', 'TRANSPOSE', n, 1_${ik}$, m, a, lda, work, d,max( 1_${ik}$, n ), work( m+&
np+1 ), lwork-m-np, info )
lopt = max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
! solve t22*y2 = d2 for y2
if( n>m ) then
call stdlib${ii}$_dtrtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', n-m, 1_${ik}$,b( m+1, m+p-n+1 ), &
ldb, d( m+1 ), n-m, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
call stdlib${ii}$_dcopy( n-m, d( m+1 ), 1_${ik}$, y( m+p-n+1 ), 1_${ik}$ )
end if
! set y1 = 0
do i = 1, m + p - n
y( i ) = zero
end do
! update d1 = d1 - t12*y2
call stdlib${ii}$_dgemv( 'NO TRANSPOSE', m, n-m, -one, b( 1_${ik}$, m+p-n+1 ), ldb,y( m+p-n+1 ), 1_${ik}$, &
one, d, 1_${ik}$ )
! solve triangular system: r11*x = d1
if( m>0_${ik}$ ) then
call stdlib${ii}$_dtrtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', m, 1_${ik}$, a, lda,d, m, info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! copy d to x
call stdlib${ii}$_dcopy( m, d, 1_${ik}$, x, 1_${ik}$ )
end if
! backward transformation y = z**t *y
call stdlib${ii}$_dormrq( 'LEFT', 'TRANSPOSE', p, 1_${ik}$, np,b( max( 1_${ik}$, n-p+1 ), 1_${ik}$ ), ldb, work( &
m+1 ), y,max( 1_${ik}$, p ), work( m+np+1 ), lwork-m-np, info )
work( 1_${ik}$ ) = m + np + max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_dggglm
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$ggglm( n, m, p, a, lda, b, ldb, d, x, y, work, lwork,info )
!! DGGGLM: solves a general Gauss-Markov linear model (GLM) problem:
!! minimize || y ||_2 subject to d = A*x + B*y
!! x
!! where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!! given N-vector. It is assumed that M <= N <= M+P, and
!! rank(A) = M and rank( A B ) = N.
!! Under these assumptions, the constrained equation is always
!! consistent, and there is a unique solution x and a minimal 2-norm
!! solution y, which is obtained using a generalized QR factorization
!! of the matrices (A, B) given by
!! A = Q*(R), B = Q*T*Z.
!! (0)
!! In particular, if matrix B is square nonsingular, then the problem
!! GLM is equivalent to the following weighted linear least squares
!! problem
!! minimize || inv(B)*(d-A*x) ||_2
!! x
!! where inv(B) denotes the inverse of B.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*), d(*)
real(${rk}$), intent(out) :: work(*), x(*), y(*)
! ===================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: i, lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3, nb4, np
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
np = min( n, p )
lquery = ( lwork==-1_${ik}$ )
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( m<0_${ik}$ .or. m>n ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGEQRF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DGERQF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMQR', ' ', n, m, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'DORMRQ', ' ', n, m, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = m + np + max( n, p )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGGGLM', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
do i = 1, m
x(i) = zero
end do
do i = 1, p
y(i) = zero
end do
return
end if
! compute the gqr factorization of matrices a and b:
! q**t*a = ( r11 ) m, q**t*b*z**t = ( t11 t12 ) m
! ( 0 ) n-m ( 0 t22 ) n-m
! m m+p-n n-m
! where r11 and t22 are upper triangular, and q and z are
! orthogonal.
call stdlib${ii}$_${ri}$ggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),work( m+np+1 ), lwork-m-&
np, info )
lopt = work( m+np+1 )
! update left-hand-side vector d = q**t*d = ( d1 ) m
! ( d2 ) n-m
call stdlib${ii}$_${ri}$ormqr( 'LEFT', 'TRANSPOSE', n, 1_${ik}$, m, a, lda, work, d,max( 1_${ik}$, n ), work( m+&
np+1 ), lwork-m-np, info )
lopt = max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
! solve t22*y2 = d2 for y2
if( n>m ) then
call stdlib${ii}$_${ri}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', n-m, 1_${ik}$,b( m+1, m+p-n+1 ), &
ldb, d( m+1 ), n-m, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
call stdlib${ii}$_${ri}$copy( n-m, d( m+1 ), 1_${ik}$, y( m+p-n+1 ), 1_${ik}$ )
end if
! set y1 = 0
do i = 1, m + p - n
y( i ) = zero
end do
! update d1 = d1 - t12*y2
call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', m, n-m, -one, b( 1_${ik}$, m+p-n+1 ), ldb,y( m+p-n+1 ), 1_${ik}$, &
one, d, 1_${ik}$ )
! solve triangular system: r11*x = d1
if( m>0_${ik}$ ) then
call stdlib${ii}$_${ri}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', m, 1_${ik}$, a, lda,d, m, info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! copy d to x
call stdlib${ii}$_${ri}$copy( m, d, 1_${ik}$, x, 1_${ik}$ )
end if
! backward transformation y = z**t *y
call stdlib${ii}$_${ri}$ormrq( 'LEFT', 'TRANSPOSE', p, 1_${ik}$, np,b( max( 1_${ik}$, n-p+1 ), 1_${ik}$ ), ldb, work( &
m+1 ), y,max( 1_${ik}$, p ), work( m+np+1 ), lwork-m-np, info )
work( 1_${ik}$ ) = m + np + max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_${ri}$ggglm
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cggglm( n, m, p, a, lda, b, ldb, d, x, y, work, lwork,info )
!! CGGGLM solves a general Gauss-Markov linear model (GLM) problem:
!! minimize || y ||_2 subject to d = A*x + B*y
!! x
!! where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!! given N-vector. It is assumed that M <= N <= M+P, and
!! rank(A) = M and rank( A B ) = N.
!! Under these assumptions, the constrained equation is always
!! consistent, and there is a unique solution x and a minimal 2-norm
!! solution y, which is obtained using a generalized QR factorization
!! of the matrices (A, B) given by
!! A = Q*(R), B = Q*T*Z.
!! (0)
!! In particular, if matrix B is square nonsingular, then the problem
!! GLM is equivalent to the following weighted linear least squares
!! problem
!! minimize || inv(B)*(d-A*x) ||_2
!! x
!! where inv(B) denotes the inverse of B.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
complex(sp), intent(inout) :: a(lda,*), b(ldb,*), d(*)
complex(sp), intent(out) :: work(*), x(*), y(*)
! ===================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: i, lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3, nb4, np
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
np = min( n, p )
lquery = ( lwork==-1_${ik}$ )
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( m<0_${ik}$ .or. m>n ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CGEQRF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CGERQF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CUNMQR', ' ', n, m, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'CUNMRQ', ' ', n, m, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = m + np + max( n, p )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGGGLM', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
do i = 1, m
x(i) = czero
end do
do i = 1, p
y(i) = czero
end do
return
end if
! compute the gqr factorization of matrices a and b:
! q**h*a = ( r11 ) m, q**h*b*z**h = ( t11 t12 ) m
! ( 0 ) n-m ( 0 t22 ) n-m
! m m+p-n n-m
! where r11 and t22 are upper triangular, and q and z are
! unitary.
call stdlib${ii}$_cggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),work( m+np+1 ), lwork-m-&
np, info )
lopt = real( work( m+np+1 ),KIND=sp)
! update left-hand-side vector d = q**h*d = ( d1 ) m
! ( d2 ) n-m
call stdlib${ii}$_cunmqr( 'LEFT', 'CONJUGATE TRANSPOSE', n, 1_${ik}$, m, a, lda, work,d, max( 1_${ik}$, n )&
, work( m+np+1 ), lwork-m-np, info )
lopt = max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
! solve t22*y2 = d2 for y2
if( n>m ) then
call stdlib${ii}$_ctrtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', n-m, 1_${ik}$,b( m+1, m+p-n+1 ), &
ldb, d( m+1 ), n-m, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
call stdlib${ii}$_ccopy( n-m, d( m+1 ), 1_${ik}$, y( m+p-n+1 ), 1_${ik}$ )
end if
! set y1 = 0
do i = 1, m + p - n
y( i ) = czero
end do
! update d1 = d1 - t12*y2
call stdlib${ii}$_cgemv( 'NO TRANSPOSE', m, n-m, -cone, b( 1_${ik}$, m+p-n+1 ), ldb,y( m+p-n+1 ), 1_${ik}$,&
cone, d, 1_${ik}$ )
! solve triangular system: r11*x = d1
if( m>0_${ik}$ ) then
call stdlib${ii}$_ctrtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', m, 1_${ik}$, a, lda,d, m, info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! copy d to x
call stdlib${ii}$_ccopy( m, d, 1_${ik}$, x, 1_${ik}$ )
end if
! backward transformation y = z**h *y
call stdlib${ii}$_cunmrq( 'LEFT', 'CONJUGATE TRANSPOSE', p, 1_${ik}$, np,b( max( 1_${ik}$, n-p+1 ), 1_${ik}$ ), &
ldb, work( m+1 ), y,max( 1_${ik}$, p ), work( m+np+1 ), lwork-m-np, info )
work( 1_${ik}$ ) = m + np + max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_cggglm
pure module subroutine stdlib${ii}$_zggglm( n, m, p, a, lda, b, ldb, d, x, y, work, lwork,info )
!! ZGGGLM solves a general Gauss-Markov linear model (GLM) problem:
!! minimize || y ||_2 subject to d = A*x + B*y
!! x
!! where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!! given N-vector. It is assumed that M <= N <= M+P, and
!! rank(A) = M and rank( A B ) = N.
!! Under these assumptions, the constrained equation is always
!! consistent, and there is a unique solution x and a minimal 2-norm
!! solution y, which is obtained using a generalized QR factorization
!! of the matrices (A, B) given by
!! A = Q*(R), B = Q*T*Z.
!! (0)
!! In particular, if matrix B is square nonsingular, then the problem
!! GLM is equivalent to the following weighted linear least squares
!! problem
!! minimize || inv(B)*(d-A*x) ||_2
!! x
!! where inv(B) denotes the inverse of B.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
complex(dp), intent(inout) :: a(lda,*), b(ldb,*), d(*)
complex(dp), intent(out) :: work(*), x(*), y(*)
! ===================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: i, lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3, nb4, np
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
np = min( n, p )
lquery = ( lwork==-1_${ik}$ )
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( m<0_${ik}$ .or. m>n ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEQRF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGERQF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMQR', ' ', n, m, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMRQ', ' ', n, m, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = m + np + max( n, p )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGGGLM', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
do i = 1, m
x(i) = czero
end do
do i = 1, p
y(i) = czero
end do
return
end if
! compute the gqr factorization of matrices a and b:
! q**h*a = ( r11 ) m, q**h*b*z**h = ( t11 t12 ) m
! ( 0 ) n-m ( 0 t22 ) n-m
! m m+p-n n-m
! where r11 and t22 are upper triangular, and q and z are
! unitary.
call stdlib${ii}$_zggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),work( m+np+1 ), lwork-m-&
np, info )
lopt = real( work( m+np+1 ),KIND=dp)
! update left-hand-side vector d = q**h*d = ( d1 ) m
! ( d2 ) n-m
call stdlib${ii}$_zunmqr( 'LEFT', 'CONJUGATE TRANSPOSE', n, 1_${ik}$, m, a, lda, work,d, max( 1_${ik}$, n )&
, work( m+np+1 ), lwork-m-np, info )
lopt = max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
! solve t22*y2 = d2 for y2
if( n>m ) then
call stdlib${ii}$_ztrtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', n-m, 1_${ik}$,b( m+1, m+p-n+1 ), &
ldb, d( m+1 ), n-m, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
call stdlib${ii}$_zcopy( n-m, d( m+1 ), 1_${ik}$, y( m+p-n+1 ), 1_${ik}$ )
end if
! set y1 = 0
do i = 1, m + p - n
y( i ) = czero
end do
! update d1 = d1 - t12*y2
call stdlib${ii}$_zgemv( 'NO TRANSPOSE', m, n-m, -cone, b( 1_${ik}$, m+p-n+1 ), ldb,y( m+p-n+1 ), 1_${ik}$,&
cone, d, 1_${ik}$ )
! solve triangular system: r11*x = d1
if( m>0_${ik}$ ) then
call stdlib${ii}$_ztrtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', m, 1_${ik}$, a, lda,d, m, info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! copy d to x
call stdlib${ii}$_zcopy( m, d, 1_${ik}$, x, 1_${ik}$ )
end if
! backward transformation y = z**h *y
call stdlib${ii}$_zunmrq( 'LEFT', 'CONJUGATE TRANSPOSE', p, 1_${ik}$, np,b( max( 1_${ik}$, n-p+1 ), 1_${ik}$ ), &
ldb, work( m+1 ), y,max( 1_${ik}$, p ), work( m+np+1 ), lwork-m-np, info )
work( 1_${ik}$ ) = m + np + max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_zggglm
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$ggglm( n, m, p, a, lda, b, ldb, d, x, y, work, lwork,info )
!! ZGGGLM: solves a general Gauss-Markov linear model (GLM) problem:
!! minimize || y ||_2 subject to d = A*x + B*y
!! x
!! where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
!! given N-vector. It is assumed that M <= N <= M+P, and
!! rank(A) = M and rank( A B ) = N.
!! Under these assumptions, the constrained equation is always
!! consistent, and there is a unique solution x and a minimal 2-norm
!! solution y, which is obtained using a generalized QR factorization
!! of the matrices (A, B) given by
!! A = Q*(R), B = Q*T*Z.
!! (0)
!! In particular, if matrix B is square nonsingular, then the problem
!! GLM is equivalent to the following weighted linear least squares
!! problem
!! minimize || inv(B)*(d-A*x) ||_2
!! x
!! where inv(B) denotes the inverse of B.
! -- lapack driver 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, ldb, lwork, m, n, p
! Array Arguments
complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*), d(*)
complex(${ck}$), intent(out) :: work(*), x(*), y(*)
! ===================================================================
! Local Scalars
logical(lk) :: lquery
integer(${ik}$) :: i, lopt, lwkmin, lwkopt, nb, nb1, nb2, nb3, nb4, np
! Intrinsic Functions
! Executable Statements
! test the input parameters
info = 0_${ik}$
np = min( n, p )
lquery = ( lwork==-1_${ik}$ )
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( m<0_${ik}$ .or. m>n ) then
info = -2_${ik}$
else if( p<0_${ik}$ .or. p<n-m ) then
info = -3_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
! calculate workspace
if( info==0_${ik}$) then
if( n==0_${ik}$ ) then
lwkmin = 1_${ik}$
lwkopt = 1_${ik}$
else
nb1 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGEQRF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb2 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZGERQF', ' ', n, m, -1_${ik}$, -1_${ik}$ )
nb3 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMQR', ' ', n, m, p, -1_${ik}$ )
nb4 = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZUNMRQ', ' ', n, m, p, -1_${ik}$ )
nb = max( nb1, nb2, nb3, nb4 )
lwkmin = m + n + p
lwkopt = m + np + max( n, p )*nb
end if
work( 1_${ik}$ ) = lwkopt
if( lwork<lwkmin .and. .not.lquery ) then
info = -12_${ik}$
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGGGLM', -info )
return
else if( lquery ) then
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
do i = 1, m
x(i) = czero
end do
do i = 1, p
y(i) = czero
end do
return
end if
! compute the gqr factorization of matrices a and b:
! q**h*a = ( r11 ) m, q**h*b*z**h = ( t11 t12 ) m
! ( 0 ) n-m ( 0 t22 ) n-m
! m m+p-n n-m
! where r11 and t22 are upper triangular, and q and z are
! unitary.
call stdlib${ii}$_${ci}$ggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),work( m+np+1 ), lwork-m-&
np, info )
lopt = real( work( m+np+1 ),KIND=${ck}$)
! update left-hand-side vector d = q**h*d = ( d1 ) m
! ( d2 ) n-m
call stdlib${ii}$_${ci}$unmqr( 'LEFT', 'CONJUGATE TRANSPOSE', n, 1_${ik}$, m, a, lda, work,d, max( 1_${ik}$, n )&
, work( m+np+1 ), lwork-m-np, info )
lopt = max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
! solve t22*y2 = d2 for y2
if( n>m ) then
call stdlib${ii}$_${ci}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', n-m, 1_${ik}$,b( m+1, m+p-n+1 ), &
ldb, d( m+1 ), n-m, info )
if( info>0_${ik}$ ) then
info = 1_${ik}$
return
end if
call stdlib${ii}$_${ci}$copy( n-m, d( m+1 ), 1_${ik}$, y( m+p-n+1 ), 1_${ik}$ )
end if
! set y1 = 0
do i = 1, m + p - n
y( i ) = czero
end do
! update d1 = d1 - t12*y2
call stdlib${ii}$_${ci}$gemv( 'NO TRANSPOSE', m, n-m, -cone, b( 1_${ik}$, m+p-n+1 ), ldb,y( m+p-n+1 ), 1_${ik}$,&
cone, d, 1_${ik}$ )
! solve triangular system: r11*x = d1
if( m>0_${ik}$ ) then
call stdlib${ii}$_${ci}$trtrs( 'UPPER', 'NO TRANSPOSE', 'NON UNIT', m, 1_${ik}$, a, lda,d, m, info )
if( info>0_${ik}$ ) then
info = 2_${ik}$
return
end if
! copy d to x
call stdlib${ii}$_${ci}$copy( m, d, 1_${ik}$, x, 1_${ik}$ )
end if
! backward transformation y = z**h *y
call stdlib${ii}$_${ci}$unmrq( 'LEFT', 'CONJUGATE TRANSPOSE', p, 1_${ik}$, np,b( max( 1_${ik}$, n-p+1 ), 1_${ik}$ ), &
ldb, work( m+1 ), y,max( 1_${ik}$, p ), work( m+np+1 ), lwork-m-np, info )
work( 1_${ik}$ ) = m + np + max( lopt, int( work( m+np+1 ),KIND=${ik}$) )
return
end subroutine stdlib${ii}$_${ci}$ggglm
#:endif
#:endfor
#:endfor
end submodule stdlib_lapack_lsq_constrained
