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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

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

CUNCSD2BY1

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

Purpose:

 CUNCSD2BY1 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 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 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 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 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 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 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 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 REAL 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:  CBBCSD 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 264 of file cuncsd2by1.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: