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

◆ zunbdb4()

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

ZUNBDB4

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

Purpose:
!>
!> ZUNBDB4 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. M-Q must be no larger than P,
!> M-P, or Q. Routines ZUNBDB1, ZUNBDB2, and ZUNBDB3 handle cases in
!> which M-Q is not the minimum dimension.
!>
!> The unitary 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 (M-Q)-by-(M-Q) 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 <= M.
!>
[in]Q
!>          Q is INTEGER
!>           The number of columns in X11 and X21. 0 <= Q <= M and
!>           M-Q <= min(P,M-P,Q).
!>
[in,out]X11
!>          X11 is COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (M-Q)
!>           The scalar factors of the elementary reflectors that define
!>           P1.
!>
[out]TAUP2
!>          TAUP2 is COMPLEX*16 array, dimension (M-Q)
!>           The scalar factors of the elementary reflectors that define
!>           P2.
!>
[out]TAUQ1
!>          TAUQ1 is COMPLEX*16 array, dimension (Q)
!>           The scalar factors of the elementary reflectors that define
!>           Q1.
!>
[out]PHANTOM
!>          PHANTOM is COMPLEX*16 array, dimension (M)
!>           The routine computes an M-by-1 column vector Y that is
!>           orthogonal to the columns of [ X11; X21 ]. PHANTOM(1:P) and
!>           PHANTOM(P+1:M) contain Householder vectors for Y(1:P) and
!>           Y(P+1:M), respectively.
!>
[out]WORK
!>          WORK is COMPLEX*16 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 ZUNCSD for details.
!>
!>  P1, P2, and Q1 are represented as products of elementary reflectors.
!>  See ZUNCSD2BY1 for details on generating P1, P2, and Q1 using ZUNGQR
!>  and ZUNGLQ.
!>
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

Definition at line 208 of file zunbdb4.f.

212*
213* -- LAPACK computational routine --
214* -- LAPACK is a software package provided by Univ. of Tennessee, --
215* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216*
217* .. Scalar Arguments ..
218 INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
219* ..
220* .. Array Arguments ..
221 DOUBLE PRECISION PHI(*), THETA(*)
222 COMPLEX*16 PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
223 $ WORK(*), X11(LDX11,*), X21(LDX21,*)
224* ..
225*
226* ====================================================================
227*
228* .. Parameters ..
229 COMPLEX*16 NEGONE, ONE, ZERO
230 parameter( negone = (-1.0d0,0.0d0), one = (1.0d0,0.0d0),
231 $ zero = (0.0d0,0.0d0) )
232* ..
233* .. Local Scalars ..
234 DOUBLE PRECISION C, S
235 INTEGER CHILDINFO, I, ILARF, IORBDB5, J, LLARF,
236 $ LORBDB5, LWORKMIN, LWORKOPT
237 LOGICAL LQUERY
238* ..
239* .. External Subroutines ..
240 EXTERNAL zlarf1f, zlarfgp, zunbdb5, zdrot, zscal,
241 $ zlacgv,
242 $ xerbla
243* ..
244* .. External Functions ..
245 DOUBLE PRECISION DZNRM2
246 EXTERNAL dznrm2
247* ..
248* .. Intrinsic Function ..
249 INTRINSIC atan2, cos, max, sin, sqrt
250* ..
251* .. Executable Statements ..
252*
253* Test input arguments
254*
255 info = 0
256 lquery = lwork .EQ. -1
257*
258 IF( m .LT. 0 ) THEN
259 info = -1
260 ELSE IF( p .LT. m-q .OR. m-p .LT. m-q ) THEN
261 info = -2
262 ELSE IF( q .LT. m-q .OR. q .GT. m ) THEN
263 info = -3
264 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
265 info = -5
266 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
267 info = -7
268 END IF
269*
270* Compute workspace
271*
272 IF( info .EQ. 0 ) THEN
273 ilarf = 2
274 llarf = max( q-1, p-1, m-p-1 )
275 iorbdb5 = 2
276 lorbdb5 = q
277 lworkopt = ilarf + llarf - 1
278 lworkopt = max( lworkopt, iorbdb5 + lorbdb5 - 1 )
279 lworkmin = lworkopt
280 work(1) = lworkopt
281 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
282 info = -14
283 END IF
284 END IF
285 IF( info .NE. 0 ) THEN
286 CALL xerbla( 'ZUNBDB4', -info )
287 RETURN
288 ELSE IF( lquery ) THEN
289 RETURN
290 END IF
291*
292* Reduce columns 1, ..., M-Q of X11 and X21
293*
294 DO i = 1, m-q
295*
296 IF( i .EQ. 1 ) THEN
297 DO j = 1, m
298 phantom(j) = zero
299 END DO
300 CALL zunbdb5( p, m-p, q, phantom(1), 1, phantom(p+1), 1,
301 $ x11, ldx11, x21, ldx21, work(iorbdb5),
302 $ lorbdb5, childinfo )
303 CALL zscal( p, negone, phantom(1), 1 )
304 CALL zlarfgp( p, phantom(1), phantom(2), 1, taup1(1) )
305 CALL zlarfgp( m-p, phantom(p+1), phantom(p+2), 1,
306 $ taup2(1) )
307 theta(i) = atan2( dble( phantom(1) ), dble( phantom(p+1) ) )
308 c = cos( theta(i) )
309 s = sin( theta(i) )
310 CALL zlarf1f( 'L', p, q, phantom(1), 1, conjg(taup1(1)),
311 $ x11,
312 $ ldx11, work(ilarf) )
313 CALL zlarf1f( 'L', m-p, q, phantom(p+1), 1,
314 $ conjg(taup2(1)),
315 $ x21, ldx21, work(ilarf) )
316 ELSE
317 CALL zunbdb5( p-i+1, m-p-i+1, q-i+1, x11(i,i-1), 1,
318 $ x21(i,i-1), 1, x11(i,i), ldx11, x21(i,i),
319 $ ldx21, work(iorbdb5), lorbdb5, childinfo )
320 CALL zscal( p-i+1, negone, x11(i,i-1), 1 )
321 CALL zlarfgp( p-i+1, x11(i,i-1), x11(i+1,i-1), 1,
322 $ taup1(i) )
323 CALL zlarfgp( m-p-i+1, x21(i,i-1), x21(i+1,i-1), 1,
324 $ taup2(i) )
325 theta(i) = atan2( dble( x11(i,i-1) ), dble( x21(i,i-1) ) )
326 c = cos( theta(i) )
327 s = sin( theta(i) )
328 CALL zlarf1f( 'L', p-i+1, q-i+1, x11(i,i-1), 1,
329 $ conjg(taup1(i)), x11(i,i), ldx11,
330 $ work(ilarf) )
331 CALL zlarf1f( 'L', m-p-i+1, q-i+1, x21(i,i-1), 1,
332 $ conjg(taup2(i)), x21(i,i), ldx21,
333 $ work(ilarf) )
334 END IF
335*
336 CALL zdrot( q-i+1, x11(i,i), ldx11, x21(i,i), ldx21, s, -c )
337 CALL zlacgv( q-i+1, x21(i,i), ldx21 )
338 CALL zlarfgp( q-i+1, x21(i,i), x21(i,i+1), ldx21, tauq1(i) )
339 c = dble( x21(i,i) )
340 CALL zlarf1f( 'R', p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
341 $ x11(i+1,i), ldx11, work(ilarf) )
342 CALL zlarf1f( 'R', m-p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
343 $ x21(i+1,i), ldx21, work(ilarf) )
344 CALL zlacgv( q-i+1, x21(i,i), ldx21 )
345 IF( i .LT. m-q ) THEN
346 s = sqrt( dznrm2( p-i, x11(i+1,i), 1 )**2
347 $ + dznrm2( m-p-i, x21(i+1,i), 1 )**2 )
348 phi(i) = atan2( s, c )
349 END IF
350*
351 END DO
352*
353* Reduce the bottom-right portion of X11 to [ I 0 ]
354*
355 DO i = m - q + 1, p
356 CALL zlacgv( q-i+1, x11(i,i), ldx11 )
357 CALL zlarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
358 CALL zlarf1f( 'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
359 $ x11(i+1,i), ldx11, work(ilarf) )
360 CALL zlarf1f( 'R', q-p, q-i+1, x11(i,i), ldx11, tauq1(i),
361 $ x21(m-q+1,i), ldx21, work(ilarf) )
362 CALL zlacgv( q-i+1, x11(i,i), ldx11 )
363 END DO
364*
365* Reduce the bottom-right portion of X21 to [ 0 I ]
366*
367 DO i = p + 1, q
368 CALL zlacgv( q-i+1, x21(m-q+i-p,i), ldx21 )
369 CALL zlarfgp( q-i+1, x21(m-q+i-p,i), x21(m-q+i-p,i+1),
370 $ ldx21,
371 $ tauq1(i) )
372 CALL zlarf1f( 'R', q-i, q-i+1, x21(m-q+i-p,i), ldx21,
373 $ tauq1(i),
374 $ x21(m-q+i-p+1,i), ldx21, work(ilarf) )
375 CALL zlacgv( q-i+1, x21(m-q+i-p,i), ldx21 )
376 END DO
377*
378 RETURN
379*
380* End of ZUNBDB4
381*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zlacgv
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:72
zlarf1f
subroutine zlarf1f(side, m, n, v, incv, tau, c, ldc, work)
ZLARF1F applies an elementary reflector to a general rectangular
Definition zlarf1f.f:157
zlarfgp
subroutine zlarfgp(n, alpha, x, incx, tau)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition zlarfgp.f:102
dznrm2
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
zdrot
subroutine zdrot(n, zx, incx, zy, incy, c, s)
ZDROT
Definition zdrot.f:98
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
zunbdb5
subroutine zunbdb5(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
ZUNBDB5
Definition zunbdb5.f:155
Here is the call graph for this function:
Here is the caller graph for this function: