LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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, and V1 are P-by-P,
 (M-P)-by-(M-P), and Q-by-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 in X.
[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.
[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 DOUBLE PRECISION 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.
[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.
[out]IWORK
          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 255 of file zuncsd2by1.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: