LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
zuncsd2by1.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine zuncsd2by1 (JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
 ZUNCSD2BY1 More...
 

Function/Subroutine Documentation

subroutine zuncsd2by1 ( character  JOBU1,
character  JOBU2,
character  JOBV1T,
integer  M,
integer  P,
integer  Q,
complex*16, dimension(ldx11,*)  X11,
integer  LDX11,
complex*16, dimension(ldx21,*)  X21,
integer  LDX21,
double precision, dimension(*)  THETA,
complex*16, dimension(ldu1,*)  U1,
integer  LDU1,
complex*16, dimension(ldu2,*)  U2,
integer  LDU2,
complex*16, dimension(ldv1t,*)  V1T,
integer  LDV1T,
complex*16, dimension(*)  WORK,
integer  LWORK,
double precision, dimension(*)  RWORK,
integer  LRWORK,
integer, dimension(*)  IWORK,
integer  INFO 
)

ZUNCSD2BY1

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

Purpose:

 ZUNCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
 orthonormal columns that has been partitioned into a 2-by-1 block
 structure:

                                [  I  0  0 ]
                                [  0  C  0 ]
          [ X11 ]   [ U1 |    ] [  0  0  0 ]
      X = [-----] = [---------] [----------] V1**T .
          [ X21 ]   [    | U2 ] [  0  0  0 ]
                                [  0  S  0 ]
                                [  0  0  I ]
 
 X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P,
 (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
 R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
 which R = MIN(P,M-P,Q,M-Q).
Parameters
[in]JOBU1
          JOBU1 is CHARACTER
           = 'Y':      U1 is computed;
           otherwise:  U1 is not computed.
[in]JOBU2
          JOBU2 is CHARACTER
           = 'Y':      U2 is computed;
           otherwise:  U2 is not computed.
[in]JOBV1T
          JOBV1T is CHARACTER
           = 'Y':      V1T is computed;
           otherwise:  V1T is not computed.
[in]M
          M is INTEGER
           The number of rows and columns in X.
[in]P
          P is INTEGER
           The number of rows in X11 and X12. 0 <= P <= M.
[in]Q
          Q is INTEGER
           The number of columns in X11 and X21. 0 <= Q <= M.
[in,out]X11
          X11 is COMPLEX*16 array, dimension (LDX11,Q)
           On entry, part of the unitary matrix whose CSD is
           desired.
[in]LDX11
          LDX11 is INTEGER
           The leading dimension of X11. LDX11 >= MAX(1,P).
[in,out]X21
          X21 is COMPLEX*16 array, dimension (LDX21,Q)
           On entry, part of the unitary matrix whose CSD is
           desired.
[in]LDX21
          LDX21 is INTEGER
           The leading dimension of X21. LDX21 >= MAX(1,M-P).
[out]THETA
          THETA is COMPLEX*16 array, dimension (R), in which R =
           MIN(P,M-P,Q,M-Q).
           C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
           S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
[out]U1
          U1 is COMPLEX*16 array, dimension (P)
           If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
[in]LDU1
          LDU1 is INTEGER
           The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
           MAX(1,P).
[out]U2
          U2 is COMPLEX*16 array, dimension (M-P)
           If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
           matrix U2.
[in]LDU2
          LDU2 is INTEGER
           The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
           MAX(1,M-P).
[out]V1T
          V1T is COMPLEX*16 array, dimension (Q)
           If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
           matrix V1**T.
[in]LDV1T
          LDV1T is INTEGER
           The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
           MAX(1,Q).
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
           On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
           If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
           ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
           define the matrix in intermediate bidiagonal-block form
           remaining after nonconvergence. INFO specifies the number
           of nonzero PHI's.
[in]LWORK
          LWORK is INTEGER
           The dimension of the array WORK.
           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]RWORK
          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
           On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
           If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
           ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
           define the matrix in intermediate bidiagonal-block form
           remaining after nonconvergence. INFO specifies the number
           of nonzero PHI's.
[in]LRWORK
          LRWORK is INTEGER
           The dimension of the array RWORK.
 
           If LRWORK = -1, then a workspace query is assumed; the routine
           only calculates the optimal size of the RWORK array, returns
           this value as the first entry of the work array, and no error
           message related to LRWORK is issued by XERBLA.
 \param[out] IWORK
 \verbatim
          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
[out]INFO
          INFO is INTEGER
           = 0:  successful exit.
           < 0:  if INFO = -i, the i-th argument had an illegal value.
           > 0:  ZBBCSD did not converge. See the description of WORK
                above for details.
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
July 2012

Definition at line 263 of file zuncsd2by1.f.

263 *
264 * -- LAPACK computational routine (version 3.5.0) --
265 * -- LAPACK is a software package provided by Univ. of Tennessee, --
266 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
267 * July 2012
268 *
269 * .. Scalar Arguments ..
270  CHARACTER jobu1, jobu2, jobv1t
271  INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
272  $ m, p, q
273  INTEGER lrwork, lrworkmin, lrworkopt
274 * ..
275 * .. Array Arguments ..
276  DOUBLE PRECISION rwork(*)
277  DOUBLE PRECISION theta(*)
278  COMPLEX*16 u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
279  $ x11(ldx11,*), x21(ldx21,*)
280  INTEGER iwork(*)
281 * ..
282 *
283 * =====================================================================
284 *
285 * .. Parameters ..
286  COMPLEX*16 one, zero
287  parameter( one = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
288 * ..
289 * .. Local Scalars ..
290  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
291  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
292  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
293  $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
294  $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
295  $ lworkmin, lworkopt, r
296  LOGICAL lquery, wantu1, wantu2, wantv1t
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL zbbcsd, zcopy, zlacpy, zlapmr, zlapmt, zunbdb1,
301  $ xerbla
302 * ..
303 * .. External Functions ..
304  LOGICAL lsame
305  EXTERNAL lsame
306 * ..
307 * .. Intrinsic Function ..
308  INTRINSIC int, max, min
309 * ..
310 * .. Executable Statements ..
311 *
312 * Test input arguments
313 *
314  info = 0
315  wantu1 = lsame( jobu1, 'Y' )
316  wantu2 = lsame( jobu2, 'Y' )
317  wantv1t = lsame( jobv1t, 'Y' )
318  lquery = lwork .EQ. -1
319 *
320  IF( m .LT. 0 ) THEN
321  info = -4
322  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
323  info = -5
324  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
325  info = -6
326  ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
327  info = -8
328  ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
329  info = -10
330  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
331  info = -13
332  ELSE IF( wantu2 .AND. ldu2 .LT. m - p ) THEN
333  info = -15
334  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
335  info = -17
336  END IF
337 *
338  r = min( p, m-p, q, m-q )
339 *
340 * Compute workspace
341 *
342 * WORK layout:
343 * |-----------------------------------------|
344 * | LWORKOPT (1) |
345 * |-----------------------------------------|
346 * | TAUP1 (MAX(1,P)) |
347 * | TAUP2 (MAX(1,M-P)) |
348 * | TAUQ1 (MAX(1,Q)) |
349 * |-----------------------------------------|
350 * | ZUNBDB WORK | ZUNGQR WORK | ZUNGLQ WORK |
351 * | | | |
352 * | | | |
353 * | | | |
354 * | | | |
355 * |-----------------------------------------|
356 * RWORK layout:
357 * |------------------|
358 * | LRWORKOPT (1) |
359 * |------------------|
360 * | PHI (MAX(1,R-1)) |
361 * |------------------|
362 * | B11D (R) |
363 * | B11E (R-1) |
364 * | B12D (R) |
365 * | B12E (R-1) |
366 * | B21D (R) |
367 * | B21E (R-1) |
368 * | B22D (R) |
369 * | B22E (R-1) |
370 * | ZBBCSD RWORK |
371 * |------------------|
372 *
373  IF( info .EQ. 0 ) THEN
374  iphi = 2
375  ib11d = iphi + max( 1, r-1 )
376  ib11e = ib11d + max( 1, r )
377  ib12d = ib11e + max( 1, r - 1 )
378  ib12e = ib12d + max( 1, r )
379  ib21d = ib12e + max( 1, r - 1 )
380  ib21e = ib21d + max( 1, r )
381  ib22d = ib21e + max( 1, r - 1 )
382  ib22e = ib22d + max( 1, r )
383  ibbcsd = ib22e + max( 1, r - 1 )
384  itaup1 = 2
385  itaup2 = itaup1 + max( 1, p )
386  itauq1 = itaup2 + max( 1, m-p )
387  iorbdb = itauq1 + max( 1, q )
388  iorgqr = itauq1 + max( 1, q )
389  iorglq = itauq1 + max( 1, q )
390  IF( r .EQ. q ) THEN
391  CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
392  $ 0, 0, work, -1, childinfo )
393  lorbdb = int( work(1) )
394  IF( p .GE. m-p ) THEN
395  CALL zungqr( p, p, q, u1, ldu1, 0, work(1), -1,
396  $ childinfo )
397  lorgqrmin = max( 1, p )
398  lorgqropt = int( work(1) )
399  ELSE
400  CALL zungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
401  $ childinfo )
402  lorgqrmin = max( 1, m-p )
403  lorgqropt = int( work(1) )
404  END IF
405  CALL zunglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
406  $ 0, work(1), -1, childinfo )
407  lorglqmin = max( 1, q-1 )
408  lorglqopt = int( work(1) )
409  CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
410  $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
411  $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
412  lbbcsd = int( rwork(1) )
413  ELSE IF( r .EQ. p ) THEN
414  CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
415  $ 0, 0, work(1), -1, childinfo )
416  lorbdb = int( work(1) )
417  IF( p-1 .GE. m-p ) THEN
418  CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
419  $ -1, childinfo )
420  lorgqrmin = max( 1, p-1 )
421  lorgqropt = int( work(1) )
422  ELSE
423  CALL zungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
424  $ childinfo )
425  lorgqrmin = max( 1, m-p )
426  lorgqropt = int( work(1) )
427  END IF
428  CALL zunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
429  $ childinfo )
430  lorglqmin = max( 1, q )
431  lorglqopt = int( work(1) )
432  CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
433  $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
434  $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
435  lbbcsd = int( rwork(1) )
436  ELSE IF( r .EQ. m-p ) THEN
437  CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
438  $ 0, 0, work(1), -1, childinfo )
439  lorbdb = int( work(1) )
440  IF( p .GE. m-p-1 ) THEN
441  CALL zungqr( p, p, q, u1, ldu1, 0, work(1), -1,
442  $ childinfo )
443  lorgqrmin = max( 1, p )
444  lorgqropt = int( work(1) )
445  ELSE
446  CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
447  $ work(1), -1, childinfo )
448  lorgqrmin = max( 1, m-p-1 )
449  lorgqropt = int( work(1) )
450  END IF
451  CALL zunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
452  $ childinfo )
453  lorglqmin = max( 1, q )
454  lorglqopt = int( work(1) )
455  CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
456  $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
457  $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
458  $ childinfo )
459  lbbcsd = int( rwork(1) )
460  ELSE
461  CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
462  $ 0, 0, 0, work(1), -1, childinfo )
463  lorbdb = m + int( work(1) )
464  IF( p .GE. m-p ) THEN
465  CALL zungqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
466  $ childinfo )
467  lorgqrmin = max( 1, p )
468  lorgqropt = int( work(1) )
469  ELSE
470  CALL zungqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
471  $ childinfo )
472  lorgqrmin = max( 1, m-p )
473  lorgqropt = int( work(1) )
474  END IF
475  CALL zunglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
476  $ childinfo )
477  lorglqmin = max( 1, q )
478  lorglqopt = int( work(1) )
479  CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
480  $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
481  $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
482  $ childinfo )
483  lbbcsd = int( rwork(1) )
484  END IF
485  lrworkmin = ibbcsd+lbbcsd-1
486  lrworkopt = lrworkmin
487  rwork(1) = lrworkopt
488  lworkmin = max( iorbdb+lorbdb-1,
489  $ iorgqr+lorgqrmin-1,
490  $ iorglq+lorglqmin-1 )
491  lworkopt = max( iorbdb+lorbdb-1,
492  $ iorgqr+lorgqropt-1,
493  $ iorglq+lorglqopt-1 )
494  work(1) = lworkopt
495  IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
496  info = -19
497  END IF
498  END IF
499  IF( info .NE. 0 ) THEN
500  CALL xerbla( 'ZUNCSD2BY1', -info )
501  RETURN
502  ELSE IF( lquery ) THEN
503  RETURN
504  END IF
505  lorgqr = lwork-iorgqr+1
506  lorglq = lwork-iorglq+1
507 *
508 * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
509 * in which R = MIN(P,M-P,Q,M-Q)
510 *
511  IF( r .EQ. q ) THEN
512 *
513 * Case 1: R = Q
514 *
515 * Simultaneously bidiagonalize X11 and X21
516 *
517  CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
518  $ rwork(iphi), work(itaup1), work(itaup2),
519  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
520 *
521 * Accumulate Householder reflectors
522 *
523  IF( wantu1 .AND. p .GT. 0 ) THEN
524  CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
525  CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
526  $ lorgqr, childinfo )
527  END IF
528  IF( wantu2 .AND. m-p .GT. 0 ) THEN
529  CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
530  CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
531  $ work(iorgqr), lorgqr, childinfo )
532  END IF
533  IF( wantv1t .AND. q .GT. 0 ) THEN
534  v1t(1,1) = one
535  DO j = 2, q
536  v1t(1,j) = zero
537  v1t(j,1) = zero
538  END DO
539  CALL zlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
540  $ ldv1t )
541  CALL zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
542  $ work(iorglq), lorglq, childinfo )
543  END IF
544 *
545 * Simultaneously diagonalize X11 and X21.
546 *
547  CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
548  $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
549  $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
550  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
551  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
552  $ childinfo )
553 *
554 * Permute rows and columns to place zero submatrices in
555 * preferred positions
556 *
557  IF( q .GT. 0 .AND. wantu2 ) THEN
558  DO i = 1, q
559  iwork(i) = m - p - q + i
560  END DO
561  DO i = q + 1, m - p
562  iwork(i) = i - q
563  END DO
564  CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
565  END IF
566  ELSE IF( r .EQ. p ) THEN
567 *
568 * Case 2: R = P
569 *
570 * Simultaneously bidiagonalize X11 and X21
571 *
572  CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
573  $ rwork(iphi), work(itaup1), work(itaup2),
574  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
575 *
576 * Accumulate Householder reflectors
577 *
578  IF( wantu1 .AND. p .GT. 0 ) THEN
579  u1(1,1) = one
580  DO j = 2, p
581  u1(1,j) = zero
582  u1(j,1) = zero
583  END DO
584  CALL zlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
585  CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
586  $ work(iorgqr), lorgqr, childinfo )
587  END IF
588  IF( wantu2 .AND. m-p .GT. 0 ) THEN
589  CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
590  CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
591  $ work(iorgqr), lorgqr, childinfo )
592  END IF
593  IF( wantv1t .AND. q .GT. 0 ) THEN
594  CALL zlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
595  CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
596  $ work(iorglq), lorglq, childinfo )
597  END IF
598 *
599 * Simultaneously diagonalize X11 and X21.
600 *
601  CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
602  $ rwork(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
603  $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
604  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
605  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
606  $ childinfo )
607 *
608 * Permute rows and columns to place identity submatrices in
609 * preferred positions
610 *
611  IF( q .GT. 0 .AND. wantu2 ) THEN
612  DO i = 1, q
613  iwork(i) = m - p - q + i
614  END DO
615  DO i = q + 1, m - p
616  iwork(i) = i - q
617  END DO
618  CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
619  END IF
620  ELSE IF( r .EQ. m-p ) THEN
621 *
622 * Case 3: R = M-P
623 *
624 * Simultaneously bidiagonalize X11 and X21
625 *
626  CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
627  $ rwork(iphi), work(itaup1), work(itaup2),
628  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
629 *
630 * Accumulate Householder reflectors
631 *
632  IF( wantu1 .AND. p .GT. 0 ) THEN
633  CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
634  CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
635  $ lorgqr, childinfo )
636  END IF
637  IF( wantu2 .AND. m-p .GT. 0 ) THEN
638  u2(1,1) = one
639  DO j = 2, m-p
640  u2(1,j) = zero
641  u2(j,1) = zero
642  END DO
643  CALL zlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
644  $ ldu2 )
645  CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
646  $ work(itaup2), work(iorgqr), lorgqr, childinfo )
647  END IF
648  IF( wantv1t .AND. q .GT. 0 ) THEN
649  CALL zlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
650  CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
651  $ work(iorglq), lorglq, childinfo )
652  END IF
653 *
654 * Simultaneously diagonalize X11 and X21.
655 *
656  CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
657  $ theta, rwork(iphi), 0, 1, v1t, ldv1t, u2, ldu2,
658  $ u1, ldu1, rwork(ib11d), rwork(ib11e),
659  $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
660  $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
661  $ rwork(ibbcsd), lbbcsd, childinfo )
662 *
663 * Permute rows and columns to place identity submatrices in
664 * preferred positions
665 *
666  IF( q .GT. r ) THEN
667  DO i = 1, r
668  iwork(i) = q - r + i
669  END DO
670  DO i = r + 1, q
671  iwork(i) = i - r
672  END DO
673  IF( wantu1 ) THEN
674  CALL zlapmt( .false., p, q, u1, ldu1, iwork )
675  END IF
676  IF( wantv1t ) THEN
677  CALL zlapmr( .false., q, q, v1t, ldv1t, iwork )
678  END IF
679  END IF
680  ELSE
681 *
682 * Case 4: R = M-Q
683 *
684 * Simultaneously bidiagonalize X11 and X21
685 *
686  CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
687  $ rwork(iphi), work(itaup1), work(itaup2),
688  $ work(itauq1), work(iorbdb), work(iorbdb+m),
689  $ lorbdb-m, childinfo )
690 *
691 * Accumulate Householder reflectors
692 *
693  IF( wantu1 .AND. p .GT. 0 ) THEN
694  CALL zcopy( p, work(iorbdb), 1, u1, 1 )
695  DO j = 2, p
696  u1(1,j) = zero
697  END DO
698  CALL zlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
699  $ ldu1 )
700  CALL zungqr( p, p, m-q, u1, ldu1, work(itaup1),
701  $ work(iorgqr), lorgqr, childinfo )
702  END IF
703  IF( wantu2 .AND. m-p .GT. 0 ) THEN
704  CALL zcopy( m-p, work(iorbdb+p), 1, u2, 1 )
705  DO j = 2, m-p
706  u2(1,j) = zero
707  END DO
708  CALL zlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
709  $ ldu2 )
710  CALL zungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
711  $ work(iorgqr), lorgqr, childinfo )
712  END IF
713  IF( wantv1t .AND. q .GT. 0 ) THEN
714  CALL zlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
715  CALL zlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
716  $ v1t(m-q+1,m-q+1), ldv1t )
717  CALL zlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
718  $ v1t(p+1,p+1), ldv1t )
719  CALL zunglq( q, q, q, v1t, ldv1t, work(itauq1),
720  $ work(iorglq), lorglq, childinfo )
721  END IF
722 *
723 * Simultaneously diagonalize X11 and X21.
724 *
725  CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
726  $ theta, rwork(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
727  $ ldv1t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
728  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
729  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
730  $ childinfo )
731 *
732 * Permute rows and columns to place identity submatrices in
733 * preferred positions
734 *
735  IF( p .GT. r ) THEN
736  DO i = 1, r
737  iwork(i) = p - r + i
738  END DO
739  DO i = r + 1, p
740  iwork(i) = i - r
741  END DO
742  IF( wantu1 ) THEN
743  CALL zlapmt( .false., p, p, u1, ldu1, iwork )
744  END IF
745  IF( wantv1t ) THEN
746  CALL zlapmr( .false., p, q, v1t, ldv1t, iwork )
747  END IF
748  END IF
749  END IF
750 *
751  RETURN
752 *
753 * End of ZUNCSD2BY1
754 *
subroutine zunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
ZUNBDB4
Definition: zunbdb4.f:215
subroutine zlapmt(FORWRD, M, N, X, LDX, K)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: zlapmt.f:106
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB1
Definition: zunbdb1.f:205
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zunbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB2
Definition: zunbdb2.f:203
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB3
Definition: zunbdb3.f:203
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
Definition: zunglq.f:129
subroutine zbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, RWORK, LRWORK, INFO)
ZBBCSD
Definition: zbbcsd.f:334
subroutine zlapmr(FORWRD, M, N, X, LDX, K)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: zlapmr.f:106

Here is the call graph for this function:

Here is the caller graph for this function: