dorcsd2by1.f

*> \brief \b DORCSD2BY1
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DORCSD2BY1 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd2by1.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd2by1.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd2by1.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
*                              X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
*                              LDV1T, WORK, LWORK, IWORK, INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          JOBU1, JOBU2, JOBV1T
*       INTEGER            INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
*      $                   M, P, Q
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   THETA(*)
*       DOUBLE PRECISION   U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
*      $                   X11(LDX11,*), X21(LDX21,*)
*       INTEGER            IWORK(*)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*>\verbatim
*>
*> 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).
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOBU1
*> \verbatim
*>          JOBU1 is CHARACTER
*>          = 'Y':      U1 is computed;
*>          otherwise:  U1 is not computed.
*> \endverbatim
*>
*> \param[in] JOBU2
*> \verbatim
*>          JOBU2 is CHARACTER
*>          = 'Y':      U2 is computed;
*>          otherwise:  U2 is not computed.
*> \endverbatim
*>
*> \param[in] JOBV1T
*> \verbatim
*>          JOBV1T is CHARACTER
*>          = 'Y':      V1T is computed;
*>          otherwise:  V1T is not computed.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows in X.
*> \endverbatim
*>
*> \param[in] P
*> \verbatim
*>          P is INTEGER
*>          The number of rows in X11. 0 <= P <= M.
*> \endverbatim
*>
*> \param[in] Q
*> \verbatim
*>          Q is INTEGER
*>          The number of columns in X11 and X21. 0 <= Q <= M.
*> \endverbatim
*>
*> \param[in,out] X11
*> \verbatim
*>          X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
*>          On entry, part of the orthogonal matrix whose CSD is desired.
*> \endverbatim
*>
*> \param[in] LDX11
*> \verbatim
*>          LDX11 is INTEGER
*>          The leading dimension of X11. LDX11 >= MAX(1,P).
*> \endverbatim
*>
*> \param[in,out] X21
*> \verbatim
*>          X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
*>          On entry, part of the orthogonal matrix whose CSD is desired.
*> \endverbatim
*>
*> \param[in] LDX21
*> \verbatim
*>          LDX21 is INTEGER
*>          The leading dimension of X21. LDX21 >= MAX(1,M-P).
*> \endverbatim
*>
*> \param[out] THETA
*> \verbatim
*>          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)) ).
*> \endverbatim
*>
*> \param[out] U1
*> \verbatim
*>          U1 is DOUBLE PRECISION array, dimension (P)
*>          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
*> \endverbatim
*>
*> \param[in] LDU1
*> \verbatim
*>          LDU1 is INTEGER
*>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
*>          MAX(1,P).
*> \endverbatim
*>
*> \param[out] U2
*> \verbatim
*>          U2 is DOUBLE PRECISION array, dimension (M-P)
*>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
*>          matrix U2.
*> \endverbatim
*>
*> \param[in] LDU2
*> \verbatim
*>          LDU2 is INTEGER
*>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
*>          MAX(1,M-P).
*> \endverbatim
*>
*> \param[out] V1T
*> \verbatim
*>          V1T is DOUBLE PRECISION array, dimension (Q)
*>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
*>          matrix V1**T.
*> \endverbatim
*>
*> \param[in] LDV1T
*> \verbatim
*>          LDV1T is INTEGER
*>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
*>          MAX(1,Q).
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          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.
*> \endverbatim
*
*> \par References:
*  ================
*>
*>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
*>      Algorithms, 50(1):33-65, 2009.
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup uncsd2by1
*
*  =====================================================================
      SUBROUTINE dorcsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
     $                       X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
     $                       LDV1T, WORK, LWORK, IWORK, INFO )
*
*  -- LAPACK computational routine (3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          JOBU1, JOBU2, JOBV1T
      INTEGER            INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
     $                   m, p, q
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   THETA(*)
      DOUBLE PRECISION   U1(ldu1,*), U2(ldu2,*), V1T(ldv1t,*), WORK(*),
     $                   x11(ldx11,*), x21(ldx21,*)
      INTEGER            IWORK(*)
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      parameter( one = 1.0d0, zero = 0.0d0 )
*     ..
*     .. Local Scalars ..
      INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
     $                   ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
     $                   iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
     $                   j, lbbcsd, lorbdb, lorglq, lorglqmin,
     $                   lorglqopt, lorgqr, lorgqrmin, lorgqropt,
     $                   lworkmin, lworkopt, r
      LOGICAL            LQUERY, WANTU1, WANTU2, WANTV1T
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   DUM1(1), DUM2(1,1)
*     ..
*     .. External Subroutines ..
      EXTERNAL           dbbcsd, dcopy, dlacpy, dlapmr, dlapmt, dorbdb1,
     $                   dorbdb2, dorbdb3, dorbdb4, dorglq, dorgqr,
     $                   xerbla
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. Intrinsic Function ..
      INTRINSIC          int, max, min
*     ..
*     .. Executable Statements ..
*
*     Test input arguments
*
      info = 0
      wantu1 = lsame( jobu1, 'Y' )
      wantu2 = lsame( jobu2, 'Y' )
      wantv1t = lsame( jobv1t, 'Y' )
      lquery = lwork .EQ. -1
*
      IF( m .LT. 0 ) THEN
         info = -4
      ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
         info = -5
      ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
         info = -6
      ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
         info = -8
      ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
         info = -10
      ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) ) THEN
         info = -13
      ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) ) THEN
         info = -15
      ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) ) THEN
         info = -17
      END IF
*
      r = min( p, m-p, q, m-q )
*
*     Compute workspace
*
*       WORK layout:
*     |-------------------------------------------------------|
*     | LWORKOPT (1)                                          |
*     |-------------------------------------------------------|
*     | PHI (MAX(1,R-1))                                      |
*     |-------------------------------------------------------|
*     | TAUP1 (MAX(1,P))                        | B11D (R)    |
*     | TAUP2 (MAX(1,M-P))                      | B11E (R-1)  |
*     | TAUQ1 (MAX(1,Q))                        | B12D (R)    |
*     |-----------------------------------------| B12E (R-1)  |
*     | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R)    |
*     |             |             |             | B21E (R-1)  |
*     |             |             |             | B22D (R)    |
*     |             |             |             | B22E (R-1)  |
*     |             |             |             | DBBCSD WORK |
*     |-------------------------------------------------------|
*
      IF( info .EQ. 0 ) THEN
         iphi = 2
         ib11d = iphi + max( 1, r-1 )
         ib11e = ib11d + max( 1, r )
         ib12d = ib11e + max( 1, r - 1 )
         ib12e = ib12d + max( 1, r )
         ib21d = ib12e + max( 1, r - 1 )
         ib21e = ib21d + max( 1, r )
         ib22d = ib21e + max( 1, r - 1 )
         ib22e = ib22d + max( 1, r )
         ibbcsd = ib22e + max( 1, r - 1 )
         itaup1 = iphi + max( 1, r-1 )
         itaup2 = itaup1 + max( 1, p )
         itauq1 = itaup2 + max( 1, m-p )
         iorbdb = itauq1 + max( 1, q )
         iorgqr = itauq1 + max( 1, q )
         iorglq = itauq1 + max( 1, q )
         lorgqrmin = 1
         lorgqropt = 1
         lorglqmin = 1
         lorglqopt = 1
         IF( r .EQ. q ) THEN
            CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                    dum1, dum1, dum1, dum1, work,
     $                    -1, childinfo )
            lorbdb = int( work(1) )
            IF( wantu1 .AND. p .GT. 0 ) THEN
               CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
     $                      childinfo )
               lorgqrmin = max( lorgqrmin, p )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            ENDIF
            IF( wantu2 .AND. m-p .GT. 0 ) THEN
               CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
     $                      -1, childinfo )
               lorgqrmin = max( lorgqrmin, m-p )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            END IF
            IF( wantv1t .AND. q .GT. 0 ) THEN
               CALL dorglq( q-1, q-1, q-1, v1t, ldv1t,
     $                      dum1, work(1), -1, childinfo )
               lorglqmin = max( lorglqmin, q-1 )
               lorglqopt = max( lorglqopt, int( work(1) ) )
            END IF
            CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
     $                   dum1, u1, ldu1, u2, ldu2, v1t, ldv1t,
     $                   dum2, 1, dum1, dum1, dum1,
     $                   dum1, dum1, dum1, dum1,
     $                   dum1, work(1), -1, childinfo )
            lbbcsd = int( work(1) )
         ELSE IF( r .EQ. p ) THEN
            CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                    dum1, dum1, dum1, dum1,
     $                    work(1), -1, childinfo )
            lorbdb = int( work(1) )
            IF( wantu1 .AND. p .GT. 0 ) THEN
               CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, dum1,
     $                      work(1), -1, childinfo )
               lorgqrmin = max( lorgqrmin, p-1 )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            END IF
            IF( wantu2 .AND. m-p .GT. 0 ) THEN
               CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
     $                      -1, childinfo )
               lorgqrmin = max( lorgqrmin, m-p )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            END IF
            IF( wantv1t .AND. q .GT. 0 ) THEN
               CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
     $                      childinfo )
               lorglqmin = max( lorglqmin, q )
               lorglqopt = max( lorglqopt, int( work(1) ) )
            END IF
            CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
     $                   dum1, v1t, ldv1t, dum2, 1, u1, ldu1,
     $                   u2, ldu2, dum1, dum1, dum1,
     $                   dum1, dum1, dum1, dum1,
     $                   dum1, work(1), -1, childinfo )
            lbbcsd = int( work(1) )
         ELSE IF( r .EQ. m-p ) THEN
            CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                    dum1, dum1, dum1, dum1,
     $                    work(1), -1, childinfo )
            lorbdb = int( work(1) )
            IF( wantu1 .AND. p .GT. 0 ) THEN
               CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
     $                      childinfo )
               lorgqrmin = max( lorgqrmin, p )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            END IF
            IF( wantu2 .AND. m-p .GT. 0 ) THEN
               CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
     $                      dum1, work(1), -1, childinfo )
               lorgqrmin = max( lorgqrmin, m-p-1 )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            END IF
            IF( wantv1t .AND. q .GT. 0 ) THEN
               CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
     $                      childinfo )
               lorglqmin = max( lorglqmin, q )
               lorglqopt = max( lorglqopt, int( work(1) ) )
            END IF
            CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
     $                   theta, dum1, dum2, 1, v1t, ldv1t, u2,
     $                   ldu2, u1, ldu1, dum1, dum1, dum1,
     $                   dum1, dum1, dum1, dum1,
     $                   dum1, work(1), -1, childinfo )
            lbbcsd = int( work(1) )
         ELSE
            CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                    dum1, dum1, dum1, dum1,
     $                    dum1, work(1), -1, childinfo )
            lorbdb = m + int( work(1) )
            IF( wantu1 .AND. p .GT. 0 ) THEN
               CALL dorgqr( p, p, m-q, u1, ldu1, dum1, work(1), -1,
     $                      childinfo )
               lorgqrmin = max( lorgqrmin, p )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            END IF
            IF( wantu2 .AND. m-p .GT. 0 ) THEN
               CALL dorgqr( m-p, m-p, m-q, u2, ldu2, dum1, work(1),
     $                      -1, childinfo )
               lorgqrmin = max( lorgqrmin, m-p )
               lorgqropt = max( lorgqropt, int( work(1) ) )
            END IF
            IF( wantv1t .AND. q .GT. 0 ) THEN
               CALL dorglq( q, q, q, v1t, ldv1t, dum1, work(1), -1,
     $                      childinfo )
               lorglqmin = max( lorglqmin, q )
               lorglqopt = max( lorglqopt, int( work(1) ) )
            END IF
            CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
     $                   theta, dum1, u2, ldu2, u1, ldu1, dum2,
     $                   1, v1t, ldv1t, dum1, dum1, dum1,
     $                   dum1, dum1, dum1, dum1,
     $                   dum1, work(1), -1, childinfo )
            lbbcsd = int( work(1) )
         END IF
         lworkmin = max( iorbdb+lorbdb-1,
     $                   iorgqr+lorgqrmin-1,
     $                   iorglq+lorglqmin-1,
     $                   ibbcsd+lbbcsd-1 )
         lworkopt = max( iorbdb+lorbdb-1,
     $                   iorgqr+lorgqropt-1,
     $                   iorglq+lorglqopt-1,
     $                   ibbcsd+lbbcsd-1 )
         work(1) = lworkopt
         IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
            info = -19
         END IF
      END IF
      IF( info .NE. 0 ) THEN
         CALL xerbla( 'DORCSD2BY1', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
      lorgqr = lwork-iorgqr+1
      lorglq = lwork-iorglq+1
*
*     Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
*     in which R = MIN(P,M-P,Q,M-Q)
*
      IF( r .EQ. q ) THEN
*
*        Case 1: R = Q
*
*        Simultaneously bidiagonalize X11 and X21
*
         CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                 work(iphi), work(itaup1), work(itaup2),
     $                 work(itauq1), work(iorbdb), lorbdb, childinfo )
*
*        Accumulate Householder reflectors
*
         IF( wantu1 .AND. p .GT. 0 ) THEN
            CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
            CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
     $                   lorgqr, childinfo )
         END IF
         IF( wantu2 .AND. m-p .GT. 0 ) THEN
            CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
            CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
     $                   work(iorgqr), lorgqr, childinfo )
         END IF
         IF( wantv1t .AND. q .GT. 0 ) THEN
            v1t(1,1) = one
            DO j = 2, q
               v1t(1,j) = zero
               v1t(j,1) = zero
            END DO
            CALL dlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
     $                   ldv1t )
            CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
     $                   work(iorglq), lorglq, childinfo )
         END IF
*
*        Simultaneously diagonalize X11 and X21.
*
         CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
     $                work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t,
     $                dum2, 1, work(ib11d), work(ib11e),
     $                work(ib12d), work(ib12e), work(ib21d),
     $                work(ib21e), work(ib22d), work(ib22e),
     $                work(ibbcsd), lbbcsd, childinfo )
*
*        Permute rows and columns to place zero submatrices in
*        preferred positions
*
         IF( q .GT. 0 .AND. wantu2 ) THEN
            DO i = 1, q
               iwork(i) = m - p - q + i
            END DO
            DO i = q + 1, m - p
               iwork(i) = i - q
            END DO
            CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
         END IF
      ELSE IF( r .EQ. p ) THEN
*
*        Case 2: R = P
*
*        Simultaneously bidiagonalize X11 and X21
*
         CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                 work(iphi), work(itaup1), work(itaup2),
     $                 work(itauq1), work(iorbdb), lorbdb, childinfo )
*
*        Accumulate Householder reflectors
*
         IF( wantu1 .AND. p .GT. 0 ) THEN
            u1(1,1) = one
            DO j = 2, p
               u1(1,j) = zero
               u1(j,1) = zero
            END DO
            CALL dlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
            CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
     $                   work(iorgqr), lorgqr, childinfo )
         END IF
         IF( wantu2 .AND. m-p .GT. 0 ) THEN
            CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
            CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
     $                   work(iorgqr), lorgqr, childinfo )
         END IF
         IF( wantv1t .AND. q .GT. 0 ) THEN
            CALL dlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
            CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
     $                   work(iorglq), lorglq, childinfo )
         END IF
*
*        Simultaneously diagonalize X11 and X21.
*
         CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
     $                work(iphi), v1t, ldv1t, dum1, 1, u1, ldu1, u2,
     $                ldu2, work(ib11d), work(ib11e), work(ib12d),
     $                work(ib12e), work(ib21d), work(ib21e),
     $                work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
     $                childinfo )
*
*        Permute rows and columns to place identity submatrices in
*        preferred positions
*
         IF( q .GT. 0 .AND. wantu2 ) THEN
            DO i = 1, q
               iwork(i) = m - p - q + i
            END DO
            DO i = q + 1, m - p
               iwork(i) = i - q
            END DO
            CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
         END IF
      ELSE IF( r .EQ. m-p ) THEN
*
*        Case 3: R = M-P
*
*        Simultaneously bidiagonalize X11 and X21
*
         CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                 work(iphi), work(itaup1), work(itaup2),
     $                 work(itauq1), work(iorbdb), lorbdb, childinfo )
*
*        Accumulate Householder reflectors
*
         IF( wantu1 .AND. p .GT. 0 ) THEN
            CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
            CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
     $                   lorgqr, childinfo )
         END IF
         IF( wantu2 .AND. m-p .GT. 0 ) THEN
            u2(1,1) = one
            DO j = 2, m-p
               u2(1,j) = zero
               u2(j,1) = zero
            END DO
            CALL dlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
     $                   ldu2 )
            CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
     $                   work(itaup2), work(iorgqr), lorgqr, childinfo )
         END IF
         IF( wantv1t .AND. q .GT. 0 ) THEN
            CALL dlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
            CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
     $                   work(iorglq), lorglq, childinfo )
         END IF
*
*        Simultaneously diagonalize X11 and X21.
*
         CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
     $                theta, work(iphi), dum1, 1, v1t, ldv1t, u2,
     $                ldu2, u1, ldu1, work(ib11d), work(ib11e),
     $                work(ib12d), work(ib12e), work(ib21d),
     $                work(ib21e), work(ib22d), work(ib22e),
     $                work(ibbcsd), lbbcsd, childinfo )
*
*        Permute rows and columns to place identity submatrices in
*        preferred positions
*
         IF( q .GT. r ) THEN
            DO i = 1, r
               iwork(i) = q - r + i
            END DO
            DO i = r + 1, q
               iwork(i) = i - r
            END DO
            IF( wantu1 ) THEN
               CALL dlapmt( .false., p, q, u1, ldu1, iwork )
            END IF
            IF( wantv1t ) THEN
               CALL dlapmr( .false., q, q, v1t, ldv1t, iwork )
            END IF
         END IF
      ELSE
*
*        Case 4: R = M-Q
*
*        Simultaneously bidiagonalize X11 and X21
*
         CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
     $                 work(iphi), work(itaup1), work(itaup2),
     $                 work(itauq1), work(iorbdb), work(iorbdb+m),
     $                 lorbdb-m, childinfo )
*
*        Accumulate Householder reflectors
*
         IF( wantu2 .AND. m-p .GT. 0 ) THEN
            CALL dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
         END IF
         IF( wantu1 .AND. p .GT. 0 ) THEN
            CALL dcopy( p, work(iorbdb), 1, u1, 1 )
            DO j = 2, p
               u1(1,j) = zero
            END DO
            CALL dlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
     $                   ldu1 )
            CALL dorgqr( p, p, m-q, u1, ldu1, work(itaup1),
     $                   work(iorgqr), lorgqr, childinfo )
         END IF
         IF( wantu2 .AND. m-p .GT. 0 ) THEN
            DO j = 2, m-p
               u2(1,j) = zero
            END DO
            CALL dlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
     $                   ldu2 )
            CALL dorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
     $                   work(iorgqr), lorgqr, childinfo )
         END IF
         IF( wantv1t .AND. q .GT. 0 ) THEN
            CALL dlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
            CALL dlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
     $                   v1t(m-q+1,m-q+1), ldv1t )
            CALL dlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
     $                   v1t(p+1,p+1), ldv1t )
            CALL dorglq( q, q, q, v1t, ldv1t, work(itauq1),
     $                   work(iorglq), lorglq, childinfo )
         END IF
*
*        Simultaneously diagonalize X11 and X21.
*
         CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
     $                theta, work(iphi), u2, ldu2, u1, ldu1, dum1,
     $                1, v1t, ldv1t, work(ib11d), work(ib11e),
     $                work(ib12d), work(ib12e), work(ib21d),
     $                work(ib21e), work(ib22d), work(ib22e),
     $                work(ibbcsd), lbbcsd, childinfo )
*
*        Permute rows and columns to place identity submatrices in
*        preferred positions
*
         IF( p .GT. r ) THEN
            DO i = 1, r
               iwork(i) = p - r + i
            END DO
            DO i = r + 1, p
               iwork(i) = i - r
            END DO
            IF( wantu1 ) THEN
               CALL dlapmt( .false., p, p, u1, ldu1, iwork )
            END IF
            IF( wantv1t ) THEN
               CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
            END IF
         END IF
      END IF
*
      RETURN
*
*     End of DORCSD2BY1
*
      END