LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dorcsd2by1()

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:
 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:

                                [  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 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). 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 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 235 of file dorcsd2by1.f.

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