LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dorcsd2by1 ( character  JOBU1,
character  JOBU2,
character  JOBV1T,
integer  M,
integer  P,
integer  Q,
double precision, dimension(ldx11,*)  X11,
integer  LDX11,
double precision, dimension(ldx21,*)  X21,
integer  LDX21,
double precision, dimension(*)  THETA,
double precision, dimension(ldu1,*)  U1,
integer  LDU1,
double precision, dimension(ldu2,*)  U2,
integer  LDU2,
double precision, dimension(ldv1t,*)  V1T,
integer  LDV1T,
double precision, dimension(*)  WORK,
integer  LWORK,
integer, dimension(*)  IWORK,
integer  INFO 
)

DORCSD2BY1

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

Purpose: 

 Purpose:
 ========

 DORCSD2BY1 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 orthogonal 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 DOUBLE PRECISION array, dimension (LDX11,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.
[in]LDX11
          LDX11 is INTEGER
          The leading dimension of X11. LDX11 >= MAX(1,P).
[in,out]X21
          X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
          On entry, part of the orthogonal 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 DOUBLE PRECISION array, dimension (P)
          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
[in]LDU1
          LDU1 is INTEGER
          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
          MAX(1,P).
[out]U2
          U2 is DOUBLE PRECISION array, dimension (M-P)
          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
          matrix U2.
[in]LDU2
          LDU2 is INTEGER
          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
          MAX(1,M-P).
[out]V1T
          V1T is DOUBLE PRECISION array, dimension (Q)
          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
          matrix V1**T.
[in]LDV1T
          LDV1T is INTEGER
          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
          MAX(1,Q).
[out]WORK
          WORK is DOUBLE PRECISION 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]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:  DBBCSD 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 236 of file dorcsd2by1.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: