LAPACK  3.10.1
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.

Definition at line 230 of file dorcsd2by1.f.

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