LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dorbdb2()

subroutine dorbdb2 ( integer m,
integer p,
integer q,
double precision, dimension(ldx11,*) x11,
integer ldx11,
double precision, dimension(ldx21,*) x21,
integer ldx21,
double precision, dimension(*) theta,
double precision, dimension(*) phi,
double precision, dimension(*) taup1,
double precision, dimension(*) taup2,
double precision, dimension(*) tauq1,
double precision, dimension(*) work,
integer lwork,
integer info )

DORBDB2

Download DORBDB2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DORBDB2 simultaneously bidiagonalizes the blocks of a tall and skinny
!> matrix X with orthonormal columns:
!>
!>                            [ B11 ]
!>      [ X11 ]   [ P1 |    ] [  0  ]
!>      [-----] = [---------] [-----] Q1**T .
!>      [ X21 ]   [    | P2 ] [ B21 ]
!>                            [  0  ]
!>
!> X11 is P-by-Q, and X21 is (M-P)-by-Q. P must be no larger than M-P,
!> Q, or M-Q. Routines DORBDB1, DORBDB3, and DORBDB4 handle cases in
!> which P is not the minimum dimension.
!>
!> The orthogonal matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),
!> and (M-Q)-by-(M-Q), respectively. They are represented implicitly by
!> Householder vectors.
!>
!> B11 and B12 are P-by-P bidiagonal matrices represented implicitly by
!> angles THETA, PHI.
!>
!>
Parameters
[in]M
!>          M is INTEGER
!>           The number of rows X11 plus the number of rows in X21.
!>
[in]P
!>          P is INTEGER
!>           The number of rows in X11. 0 <= P <= min(M-P,Q,M-Q).
!>
[in]Q
!>          Q is INTEGER
!>           The number of columns in X11 and X21. 0 <= Q <= M.
!>
[in,out]X11
!>          X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
!>           On entry, the top block of the matrix X to be reduced. On
!>           exit, the columns of tril(X11) specify reflectors for P1 and
!>           the rows of triu(X11,1) specify reflectors for Q1.
!>
[in]LDX11
!>          LDX11 is INTEGER
!>           The leading dimension of X11. LDX11 >= P.
!>
[in,out]X21
!>          X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
!>           On entry, the bottom block of the matrix X to be reduced. On
!>           exit, the columns of tril(X21) specify reflectors for P2.
!>
[in]LDX21
!>          LDX21 is INTEGER
!>           The leading dimension of X21. LDX21 >= M-P.
!>
[out]THETA
!>          THETA is DOUBLE PRECISION array, dimension (Q)
!>           The entries of the bidiagonal blocks B11, B21 are defined by
!>           THETA and PHI. See Further Details.
!>
[out]PHI
!>          PHI is DOUBLE PRECISION array, dimension (Q-1)
!>           The entries of the bidiagonal blocks B11, B21 are defined by
!>           THETA and PHI. See Further Details.
!>
[out]TAUP1
!>          TAUP1 is DOUBLE PRECISION array, dimension (P-1)
!>           The scalar factors of the elementary reflectors that define
!>           P1.
!>
[out]TAUP2
!>          TAUP2 is DOUBLE PRECISION array, dimension (Q)
!>           The scalar factors of the elementary reflectors that define
!>           P2.
!>
[out]TAUQ1
!>          TAUQ1 is DOUBLE PRECISION array, dimension (Q)
!>           The scalar factors of the elementary reflectors that define
!>           Q1.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (LWORK)
!>
[in]LWORK
!>          LWORK is INTEGER
!>           The dimension of the array WORK. LWORK >= M-Q.
!>
!>           If LWORK = -1, then a workspace query is assumed; the routine
!>           only calculates the optimal size of the WORK array, returns
!>           this value as the first entry of the WORK array, and no error
!>           message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>           = 0:  successful exit.
!>           < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The upper-bidiagonal blocks B11, B21 are represented implicitly by
!>  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
!>  in each bidiagonal band is a product of a sine or cosine of a THETA
!>  with a sine or cosine of a PHI. See [1] or DORCSD for details.
!>
!>  P1, P2, and Q1 are represented as products of elementary reflectors.
!>  See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR
!>  and DORGLQ.
!>
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

Definition at line 198 of file dorbdb2.f.

201*
202* -- LAPACK computational routine --
203* -- LAPACK is a software package provided by Univ. of Tennessee, --
204* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
205*
206* .. Scalar Arguments ..
207 INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
208* ..
209* .. Array Arguments ..
210 DOUBLE PRECISION PHI(*), THETA(*)
211 DOUBLE PRECISION TAUP1(*), TAUP2(*), TAUQ1(*), WORK(*),
212 $ X11(LDX11,*), X21(LDX21,*)
213* ..
214*
215* ====================================================================
216*
217* .. Parameters ..
218 DOUBLE PRECISION NEGONE, ONE
219 parameter( negone = -1.0d0, one = 1.0d0 )
220* ..
221* .. Local Scalars ..
222 DOUBLE PRECISION C, S
223 INTEGER CHILDINFO, I, ILARF, IORBDB5, LLARF, LORBDB5,
224 $ LWORKMIN, LWORKOPT
225 LOGICAL LQUERY
226* ..
227* .. External Subroutines ..
228 EXTERNAL dlarf1f, dlarfgp, dorbdb5, drot, dscal,
229 $ xerbla
230* ..
231* .. External Functions ..
232 DOUBLE PRECISION DNRM2
233 EXTERNAL dnrm2
234* ..
235* .. Intrinsic Function ..
236 INTRINSIC atan2, cos, max, sin, sqrt
237* ..
238* .. Executable Statements ..
239*
240* Test input arguments
241*
242 info = 0
243 lquery = lwork .EQ. -1
244*
245 IF( m .LT. 0 ) THEN
246 info = -1
247 ELSE IF( p .LT. 0 .OR. p .GT. m-p ) THEN
248 info = -2
249 ELSE IF( q .LT. 0 .OR. q .LT. p .OR. m-q .LT. p ) THEN
250 info = -3
251 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
252 info = -5
253 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
254 info = -7
255 END IF
256*
257* Compute workspace
258*
259 IF( info .EQ. 0 ) THEN
260 ilarf = 2
261 llarf = max( p-1, m-p, q-1 )
262 iorbdb5 = 2
263 lorbdb5 = q-1
264 lworkopt = max( ilarf+llarf-1, iorbdb5+lorbdb5-1 )
265 lworkmin = lworkopt
266 work(1) = lworkopt
267 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
268 info = -14
269 END IF
270 END IF
271 IF( info .NE. 0 ) THEN
272 CALL xerbla( 'DORBDB2', -info )
273 RETURN
274 ELSE IF( lquery ) THEN
275 RETURN
276 END IF
277*
278* Reduce rows 1, ..., P of X11 and X21
279*
280 DO i = 1, p
281*
282 IF( i .GT. 1 ) THEN
283 CALL drot( q-i+1, x11(i,i), ldx11, x21(i-1,i), ldx21, c,
284 $ s )
285 END IF
286 CALL dlarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
287 c = x11(i,i)
288 CALL dlarf1f( 'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
289 $ x11(i+1,i), ldx11, work(ilarf) )
290 CALL dlarf1f( 'R', m-p-i+1, q-i+1, x11(i,i), ldx11,
291 $ tauq1(i), x21(i,i), ldx21, work(ilarf) )
292 s = sqrt( dnrm2( p-i, x11(i+1,i), 1 )**2
293 $ + dnrm2( m-p-i+1, x21(i,i), 1 )**2 )
294 theta(i) = atan2( s, c )
295*
296 CALL dorbdb5( p-i, m-p-i+1, q-i, x11(i+1,i), 1, x21(i,i), 1,
297 $ x11(i+1,i+1), ldx11, x21(i,i+1), ldx21,
298 $ work(iorbdb5), lorbdb5, childinfo )
299 CALL dscal( p-i, negone, x11(i+1,i), 1 )
300 CALL dlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1, taup2(i) )
301 IF( i .LT. p ) THEN
302 CALL dlarfgp( p-i, x11(i+1,i), x11(i+2,i), 1, taup1(i) )
303 phi(i) = atan2( x11(i+1,i), x21(i,i) )
304 c = cos( phi(i) )
305 s = sin( phi(i) )
306 CALL dlarf1f( 'L', p-i, q-i, x11(i+1,i), 1, taup1(i),
307 $ x11(i+1,i+1), ldx11, work(ilarf) )
308 END IF
309 CALL dlarf1f( 'L', m-p-i+1, q-i, x21(i,i), 1, taup2(i),
310 $ x21(i,i+1), ldx21, work(ilarf) )
311*
312 END DO
313*
314* Reduce the bottom-right portion of X21 to the identity matrix
315*
316 DO i = p + 1, q
317 CALL dlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1, taup2(i) )
318 CALL dlarf1f( 'L', m-p-i+1, q-i, x21(i,i), 1, taup2(i),
319 $ x21(i,i+1), ldx21, work(ilarf) )
320 END DO
321*
322 RETURN
323*
324* End of DORBDB2
325*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlarf1f
subroutine dlarf1f(side, m, n, v, incv, tau, c, ldc, work)
DLARF1F applies an elementary reflector to a general rectangular
Definition dlarf1f.f:157
dlarfgp
subroutine dlarfgp(n, alpha, x, incx, tau)
DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition dlarfgp.f:102
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dorbdb5
subroutine dorbdb5(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
DORBDB5
Definition dorbdb5.f:155
Here is the call graph for this function:
Here is the caller graph for this function: