LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cuncsd2by1()

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:

                                [  I1 0  0 ]
                                [  0  C  0 ]
          [ X11 ]   [ U1 |    ] [  0  0  0 ]
      X = [-----] = [---------] [----------] V1**T .
          [ X21 ]   [    | U2 ] [  0  0  0 ]
                                [  0  S  0 ]
                                [  0  0  I2]

 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). I1 is a K1-by-K1 identity matrix and I2 is a
 K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
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 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 REAL 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.
[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 and RWORK
          arrays, returns this value as the first entry of the WORK
          and RWORK array, respectively, and no error message related
          to LWORK or LRWORK 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 WORK and RWORK
          arrays, returns this value as the first entry of the WORK
          and RWORK array, respectively, and no error message related
          to LWORK or 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:  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.

Definition at line 253 of file cuncsd2by1.f.

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