LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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, 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 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 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.
[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.
Date
June 2016

Definition at line 256 of file cuncsd2by1.f.

256 *
257 * -- LAPACK computational routine (version 3.6.1) --
258 * -- LAPACK is a software package provided by Univ. of Tennessee, --
259 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
260 * June 2016
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,
297  $ cunbdb2, cunbdb3, cunbdb4, cunglq, cungqr,
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
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  END IF
516  IF( info .NE. 0 ) THEN
517  CALL xerbla( 'CUNCSD2BY1', -info )
518  RETURN
519  ELSE IF( lquery ) THEN
520  RETURN
521  END IF
522  lorgqr = lwork-iorgqr+1
523  lorglq = lwork-iorglq+1
524 *
525 * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
526 * in which R = MIN(P,M-P,Q,M-Q)
527 *
528  IF( r .EQ. q ) THEN
529 *
530 * Case 1: R = Q
531 *
532 * Simultaneously bidiagonalize X11 and X21
533 *
534  CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
535  $ rwork(iphi), work(itaup1), work(itaup2),
536  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
537 *
538 * Accumulate Householder reflectors
539 *
540  IF( wantu1 .AND. p .GT. 0 ) THEN
541  CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
542  CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
543  $ lorgqr, childinfo )
544  END IF
545  IF( wantu2 .AND. m-p .GT. 0 ) THEN
546  CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
547  CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
548  $ work(iorgqr), lorgqr, childinfo )
549  END IF
550  IF( wantv1t .AND. q .GT. 0 ) THEN
551  v1t(1,1) = one
552  DO j = 2, q
553  v1t(1,j) = zero
554  v1t(j,1) = zero
555  END DO
556  CALL clacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
557  $ ldv1t )
558  CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
559  $ work(iorglq), lorglq, childinfo )
560  END IF
561 *
562 * Simultaneously diagonalize X11 and X21.
563 *
564  CALL cbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
565  $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
566  $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
567  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
568  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
569  $ childinfo )
570 *
571 * Permute rows and columns to place zero submatrices in
572 * preferred positions
573 *
574  IF( q .GT. 0 .AND. wantu2 ) THEN
575  DO i = 1, q
576  iwork(i) = m - p - q + i
577  END DO
578  DO i = q + 1, m - p
579  iwork(i) = i - q
580  END DO
581  CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
582  END IF
583  ELSE IF( r .EQ. p ) THEN
584 *
585 * Case 2: R = P
586 *
587 * Simultaneously bidiagonalize X11 and X21
588 *
589  CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
590  $ rwork(iphi), work(itaup1), work(itaup2),
591  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
592 *
593 * Accumulate Householder reflectors
594 *
595  IF( wantu1 .AND. p .GT. 0 ) THEN
596  u1(1,1) = one
597  DO j = 2, p
598  u1(1,j) = zero
599  u1(j,1) = zero
600  END DO
601  CALL clacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
602  CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
603  $ work(iorgqr), lorgqr, childinfo )
604  END IF
605  IF( wantu2 .AND. m-p .GT. 0 ) THEN
606  CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
607  CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
608  $ work(iorgqr), lorgqr, childinfo )
609  END IF
610  IF( wantv1t .AND. q .GT. 0 ) THEN
611  CALL clacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
612  CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
613  $ work(iorglq), lorglq, childinfo )
614  END IF
615 *
616 * Simultaneously diagonalize X11 and X21.
617 *
618  CALL cbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
619  $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
620  $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
621  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
622  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
623  $ childinfo )
624 *
625 * Permute rows and columns to place identity submatrices in
626 * preferred positions
627 *
628  IF( q .GT. 0 .AND. wantu2 ) THEN
629  DO i = 1, q
630  iwork(i) = m - p - q + i
631  END DO
632  DO i = q + 1, m - p
633  iwork(i) = i - q
634  END DO
635  CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
636  END IF
637  ELSE IF( r .EQ. m-p ) THEN
638 *
639 * Case 3: R = M-P
640 *
641 * Simultaneously bidiagonalize X11 and X21
642 *
643  CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
644  $ rwork(iphi), work(itaup1), work(itaup2),
645  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
646 *
647 * Accumulate Householder reflectors
648 *
649  IF( wantu1 .AND. p .GT. 0 ) THEN
650  CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
651  CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
652  $ lorgqr, childinfo )
653  END IF
654  IF( wantu2 .AND. m-p .GT. 0 ) THEN
655  u2(1,1) = one
656  DO j = 2, m-p
657  u2(1,j) = zero
658  u2(j,1) = zero
659  END DO
660  CALL clacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
661  $ ldu2 )
662  CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
663  $ work(itaup2), work(iorgqr), lorgqr, childinfo )
664  END IF
665  IF( wantv1t .AND. q .GT. 0 ) THEN
666  CALL clacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
667  CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
668  $ work(iorglq), lorglq, childinfo )
669  END IF
670 *
671 * Simultaneously diagonalize X11 and X21.
672 *
673  CALL cbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
674  $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
675  $ u1, ldu1, rwork(ib11d), rwork(ib11e),
676  $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
677  $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
678  $ rwork(ibbcsd), lbbcsd, childinfo )
679 *
680 * Permute rows and columns to place identity submatrices in
681 * preferred positions
682 *
683  IF( q .GT. r ) THEN
684  DO i = 1, r
685  iwork(i) = q - r + i
686  END DO
687  DO i = r + 1, q
688  iwork(i) = i - r
689  END DO
690  IF( wantu1 ) THEN
691  CALL clapmt( .false., p, q, u1, ldu1, iwork )
692  END IF
693  IF( wantv1t ) THEN
694  CALL clapmr( .false., q, q, v1t, ldv1t, iwork )
695  END IF
696  END IF
697  ELSE
698 *
699 * Case 4: R = M-Q
700 *
701 * Simultaneously bidiagonalize X11 and X21
702 *
703  CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
704  $ rwork(iphi), work(itaup1), work(itaup2),
705  $ work(itauq1), work(iorbdb), work(iorbdb+m),
706  $ lorbdb-m, childinfo )
707 *
708 * Accumulate Householder reflectors
709 *
710  IF( wantu1 .AND. p .GT. 0 ) THEN
711  CALL ccopy( p, work(iorbdb), 1, u1, 1 )
712  DO j = 2, p
713  u1(1,j) = zero
714  END DO
715  CALL clacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
716  $ ldu1 )
717  CALL cungqr( p, p, m-q, u1, ldu1, work(itaup1),
718  $ work(iorgqr), lorgqr, childinfo )
719  END IF
720  IF( wantu2 .AND. m-p .GT. 0 ) THEN
721  CALL ccopy( m-p, work(iorbdb+p), 1, u2, 1 )
722  DO j = 2, m-p
723  u2(1,j) = zero
724  END DO
725  CALL clacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
726  $ ldu2 )
727  CALL cungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
728  $ work(iorgqr), lorgqr, childinfo )
729  END IF
730  IF( wantv1t .AND. q .GT. 0 ) THEN
731  CALL clacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
732  CALL clacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
733  $ v1t(m-q+1,m-q+1), ldv1t )
734  CALL clacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
735  $ v1t(p+1,p+1), ldv1t )
736  CALL cunglq( q, q, q, v1t, ldv1t, work(itauq1),
737  $ work(iorglq), lorglq, childinfo )
738  END IF
739 *
740 * Simultaneously diagonalize X11 and X21.
741 *
742  CALL cbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
743  $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
744  $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
745  $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
746  $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
747  $ rwork(ibbcsd), lbbcsd, childinfo )
748 *
749 * Permute rows and columns to place identity submatrices in
750 * preferred positions
751 *
752  IF( p .GT. r ) THEN
753  DO i = 1, r
754  iwork(i) = p - r + i
755  END DO
756  DO i = r + 1, p
757  iwork(i) = i - r
758  END DO
759  IF( wantu1 ) THEN
760  CALL clapmt( .false., p, p, u1, ldu1, iwork )
761  END IF
762  IF( wantv1t ) THEN
763  CALL clapmr( .false., p, q, v1t, ldv1t, iwork )
764  END IF
765  END IF
766  END IF
767 *
768  RETURN
769 *
770 * End of CUNCSD2BY1
771 *
cunbdb4
subroutine cunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
CUNBDB4
Definition: cunbdb4.f:215
cbbcsd
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cunbdb1
subroutine cunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB1
Definition: cunbdb1.f:204
cunglq
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
Definition: cunglq.f:129
cunbdb2
subroutine cunbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB2
Definition: cunbdb2.f:204
clapmr
subroutine clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: clapmr.f:106
clacpy
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
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
cunbdb3
subroutine cunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB3
Definition: cunbdb3.f:204
cungqr
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:130
clapmt
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
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: