libflame  revision_anchor
Functions
dorcsd2by1.c File Reference

(r)

Functions

int dorcsd2by1_ (char *jobu1, char *jobu2, char *jobv1t, integer *m, integer *p, integer *q, doublereal *x11, integer *ldx11, doublereal *x21, integer *ldx21, doublereal *theta, doublereal *u1, integer *ldu1, doublereal *u2, integer *ldu2, doublereal *v1t, integer *ldv1t, doublereal *work, integer *lwork, integer *iwork, integer *info)

Function Documentation

int dorcsd2by1_ ( char *  jobu1,
char *  jobu2,
char *  jobv1t,
integer m,
integer p,
integer q,
doublereal x11,
integer ldx11,
doublereal x21,
integer ldx21,
doublereal theta,
doublereal u1,
integer ldu1,
doublereal u2,
integer ldu2,
doublereal v1t,
integer ldv1t,
doublereal work,
integer lwork,
integer iwork,
integer info 
)

References dorglq_fla(), and dorgqr_fla().

{
    /* System generated locals */
    integer u1_dim1, u1_offset, u2_dim1, u2_offset, v1t_dim1, v1t_offset, x11_dim1, x11_offset, x21_dim1, x21_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, i__9;
    /* Local variables */
    integer lworkmin, lworkopt, i__, j, r__, childinfo, lorglqmin, lorgqrmin, lorglqopt, lorgqropt, ib11d, ib11e, ib12d, ib12e, ib21d, ib21e, ib22d, ib22e, iphi;
    extern logical lsame_(char *, char *);
    extern /* Subroutine */
    int dcopy_(integer *, doublereal *, integer *, doublereal *, integer *);
    integer itaup1, itaup2, itauq1;
    logical wantu1, wantu2;
    extern /* Subroutine */
    int dbbcsd_();
    integer ibbcsd, lbbcsd, iorbdb, lorbdb;
    extern /* Subroutine */
    int dlacpy_(char *, integer *, integer *, doublereal *, integer *, doublereal *, integer *), xerbla_(char *, integer *), dlapmr_(logical *, integer *, integer *, doublereal *, integer *, integer *), dlapmt_(logical *, integer *, integer *, doublereal *, integer *, integer *);
    integer iorglq;
    extern int /* Subroutine */
      dorglq_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *),
      dorgqr_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *);
    integer lorglq, iorgqr, lorgqr;
    extern /* Subroutine */
    int dorbdb1_(), dorbdb2_(), dorbdb3_(), dorbdb4_() ;
    logical lquery, wantv1t;
    /* -- 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..-- */
    /* July 2012 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. Intrinsic Function .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Test input arguments */
    /* Parameter adjustments */
    x11_dim1 = *ldx11;
    x11_offset = 1 + x11_dim1;
    x11 -= x11_offset;
    x21_dim1 = *ldx21;
    x21_offset = 1 + x21_dim1;
    x21 -= x21_offset;
    --theta;
    u1_dim1 = *ldu1;
    u1_offset = 1 + u1_dim1;
    u1 -= u1_offset;
    u2_dim1 = *ldu2;
    u2_offset = 1 + u2_dim1;
    u2 -= u2_offset;
    v1t_dim1 = *ldv1t;
    v1t_offset = 1 + v1t_dim1;
    v1t -= v1t_offset;
    --work;
    --iwork;
    /* Function Body */
    *info = 0;
    wantu1 = lsame_(jobu1, "Y");
    wantu2 = lsame_(jobu2, "Y");
    wantv1t = lsame_(jobv1t, "Y");
    lquery = *lwork == -1;
    if (*m < 0)
    {
        *info = -4;
    }
    else if (*p < 0 || *p > *m)
    {
        *info = -5;
    }
    else if (*q < 0 || *q > *m)
    {
        *info = -6;
    }
    else if (*ldx11 < max(1,*p))
    {
        *info = -8;
    }
    else /* if(complicated condition) */
    {
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *p; // , expr subst
        if (*ldx21 < max(i__1,i__2))
        {
            *info = -10;
        }
        else if (wantu1 && *ldu1 < *p)
        {
            *info = -13;
        }
        else if (wantu2 && *ldu2 < *m - *p)
        {
            *info = -15;
        }
        else if (wantv1t && *ldv1t < *q)
        {
            *info = -17;
        }
    }
    /* Computing MIN */
    i__1 = *p, i__2 = *m - *p, i__1 = min(i__1,i__2);
    i__1 = min(i__1,*q);
    i__2 = *m - *q; // ; expr subst
    r__ = min(i__1,i__2);
    /* 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 == 0)
    {
        iphi = 2;
        /* Computing MAX */
        i__1 = 1;
        i__2 = r__ - 1; // , expr subst
        ib11d = iphi + max(i__1,i__2);
        ib11e = ib11d + max(1,r__);
        /* Computing MAX */
        i__1 = 1;
        i__2 = r__ - 1; // , expr subst
        ib12d = ib11e + max(i__1,i__2);
        ib12e = ib12d + max(1,r__);
        /* Computing MAX */
        i__1 = 1;
        i__2 = r__ - 1; // , expr subst
        ib21d = ib12e + max(i__1,i__2);
        ib21e = ib21d + max(1,r__);
        /* Computing MAX */
        i__1 = 1;
        i__2 = r__ - 1; // , expr subst
        ib22d = ib21e + max(i__1,i__2);
        ib22e = ib22d + max(1,r__);
        /* Computing MAX */
        i__1 = 1;
        i__2 = r__ - 1; // , expr subst
        ibbcsd = ib22e + max(i__1,i__2);
        /* Computing MAX */
        i__1 = 1;
        i__2 = r__ - 1; // , expr subst
        itaup1 = iphi + max(i__1,i__2);
        itaup2 = itaup1 + max(1,*p);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *p; // , expr subst
        itauq1 = itaup2 + max(i__1,i__2);
        iorbdb = itauq1 + max(1,*q);
        iorgqr = itauq1 + max(1,*q);
        iorglq = itauq1 + max(1,*q);
        if (r__ == *q)
        {
            dorbdb1_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, &theta[1], &c__0, &c__0, &c__0, &c__0, &work[1], & c_n1, &childinfo);
            lorbdb = (integer) work[1];
            if (*p >= *m - *p)
            {
                dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
                lorgqrmin = max(1,*p);
                lorgqropt = (integer) work[1];
            }
            else
            {
                i__1 = *m - *p;
                i__2 = *m - *p;
                dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, (doublereal*)&c__0, &work[1] , &c_n1, &childinfo);
                /* Computing MAX */
                i__1 = 1;
                i__2 = *m - *p; // , expr subst
                lorgqrmin = max(i__1,i__2);
                lorgqropt = (integer) work[1];
            }
            /* Computing MAX */
            i__2 = 0;
            i__3 = *q - 1; // , expr subst
            i__1 = max(i__2,i__3);
            /* Computing MAX */
            i__5 = 0;
            i__6 = *q - 1; // , expr subst
            i__4 = max(i__5,i__6);
            /* Computing MAX */
            i__8 = 0;
            i__9 = *q - 1; // , expr subst
            i__7 = max(i__8,i__9);
            dorglq_fla(&i__1, &i__4, &i__7, &v1t[v1t_offset], ldv1t, (doublereal*)&c__0, & work[1], &c_n1, &childinfo);
            /* Computing MAX */
            i__1 = 1;
            i__2 = *q - 1; // , expr subst
            lorglqmin = max(i__1,i__2);
            lorglqopt = (integer) work[1];
            dbbcsd_(jobu1, jobu2, jobv1t, "N", "N", m, p, q, &theta[1], &c__0, &u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &v1t[ v1t_offset], ldv1t, &c__0, &c__1, &c__0, &c__0, &c__0, & c__0, &c__0, &c__0, &c__0, &c__0, &work[1], &c_n1, & childinfo);
            lbbcsd = (integer) work[1];
        }
        else if (r__ == *p)
        {
            dorbdb2_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, &theta[1], &c__0, &c__0, &c__0, &c__0, &work[1], & c_n1, &childinfo);
            lorbdb = (integer) work[1];
            if (*p - 1 >= *m - *p)
            {
                i__1 = *p - 1;
                i__2 = *p - 1;
                i__3 = *p - 1;
                dorgqr_fla(&i__1, &i__2, &i__3, &u1[(u1_dim1 << 1) + 2], ldu1, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
                /* Computing MAX */
                i__1 = 1;
                i__2 = *p - 1; // , expr subst
                lorgqrmin = max(i__1,i__2);
                lorgqropt = (integer) work[1];
            }
            else
            {
                i__1 = *m - *p;
                i__2 = *m - *p;
                dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, (doublereal*)&c__0, &work[1] , &c_n1, &childinfo);
                /* Computing MAX */
                i__1 = 1;
                i__2 = *m - *p; // , expr subst
                lorgqrmin = max(i__1,i__2);
                lorgqropt = (integer) work[1];
            }
            dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, (doublereal*)&c__0, &work[1], & c_n1, &childinfo);
            lorglqmin = max(1,*q);
            lorglqopt = (integer) work[1];
            dbbcsd_(jobv1t, "N", jobu1, jobu2, "T", m, q, p, &theta[1], &c__0, &v1t[v1t_offset], ldv1t, &c__0, &c__1, &u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &c__0, &c__0, &c__0, &c__0, & c__0, &c__0, &c__0, &c__0, &work[1], &c_n1, &childinfo);
            lbbcsd = (integer) work[1];
        }
        else if (r__ == *m - *p)
        {
            dorbdb3_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, &theta[1], &c__0, &c__0, &c__0, &c__0, &work[1], & c_n1, &childinfo);
            lorbdb = (integer) work[1];
            if (*p >= *m - *p - 1)
            {
                dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
                lorgqrmin = max(1,*p);
                lorgqropt = (integer) work[1];
            }
            else
            {
                i__1 = *m - *p - 1;
                i__2 = *m - *p - 1;
                i__3 = *m - *p - 1;
                dorgqr_fla(&i__1, &i__2, &i__3, &u2[(u2_dim1 << 1) + 2], ldu2, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
                /* Computing MAX */
                i__1 = 1;
                i__2 = *m - *p - 1; // , expr subst
                lorgqrmin = max(i__1,i__2);
                lorgqropt = (integer) work[1];
            }
            dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, (doublereal*)&c__0, &work[1], & c_n1, &childinfo);
            lorglqmin = max(1,*q);
            lorglqopt = (integer) work[1];
            i__1 = *m - *q;
            i__2 = *m - *p;
            dbbcsd_("N", jobv1t, jobu2, jobu1, "T", m, &i__1, &i__2, &theta[1] , &c__0, &c__0, &c__1, &v1t[v1t_offset], ldv1t, &u2[ u2_offset], ldu2, &u1[u1_offset], ldu1, &c__0, &c__0, & c__0, &c__0, &c__0, &c__0, &c__0, &c__0, &work[1], &c_n1, &childinfo);
            lbbcsd = (integer) work[1];
        }
        else
        {
            dorbdb4_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, &theta[1], &c__0, &c__0, &c__0, &c__0, &c__0, & work[1], &c_n1, &childinfo);
            lorbdb = *m + (integer) work[1];
            if (*p >= *m - *p)
            {
                i__1 = *m - *q;
                dorgqr_fla(p, p, &i__1, &u1[u1_offset], ldu1, (doublereal*)&c__0, &work[1], & c_n1, &childinfo);
                lorgqrmin = max(1,*p);
                lorgqropt = (integer) work[1];
            }
            else
            {
                i__1 = *m - *p;
                i__2 = *m - *p;
                i__3 = *m - *q;
                dorgqr_fla(&i__1, &i__2, &i__3, &u2[u2_offset], ldu2, (doublereal*)&c__0, & work[1], &c_n1, &childinfo);
                /* Computing MAX */
                i__1 = 1;
                i__2 = *m - *p; // , expr subst
                lorgqrmin = max(i__1,i__2);
                lorgqropt = (integer) work[1];
            }
            dorglq_fla(q, q, q, &v1t[v1t_offset], ldv1t, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
            lorglqmin = max(1,*q);
            lorglqopt = (integer) work[1];
            i__1 = *m - *p;
            i__2 = *m - *q;
            dbbcsd_(jobu2, jobu1, "N", jobv1t, "N", m, &i__1, &i__2, &theta[1] , &c__0, &u2[u2_offset], ldu2, &u1[u1_offset], ldu1, & c__0, &c__1, &v1t[v1t_offset], ldv1t, &c__0, &c__0, &c__0, &c__0, &c__0, &c__0, &c__0, &c__0, &work[1], &c_n1, & childinfo);
            lbbcsd = (integer) work[1];
        }
        /* Computing MAX */
        i__1 = iorbdb + lorbdb - 1, i__2 = iorgqr + lorgqrmin - 1, i__1 = max( i__1,i__2), i__2 = iorglq + lorglqmin - 1;
        i__1 = max(i__1, i__2);
        i__2 = ibbcsd + lbbcsd - 1; // ; expr subst
        lworkmin = max(i__1,i__2);
        /* Computing MAX */
        i__1 = iorbdb + lorbdb - 1, i__2 = iorgqr + lorgqropt - 1, i__1 = max( i__1,i__2), i__2 = iorglq + lorglqopt - 1;
        i__1 = max(i__1, i__2);
        i__2 = ibbcsd + lbbcsd - 1; // ; expr subst
        lworkopt = max(i__1,i__2);
        work[1] = (doublereal) lworkopt;
        if (*lwork < lworkmin && ! lquery)
        {
            *info = -19;
        }
    }
    if (*info != 0)
    {
        i__1 = -(*info);
        xerbla_("DORCSD2BY1", &i__1);
        return 0;
    }
    else if (lquery)
    {
        return 0;
    }
    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__ == *q)
    {
        /* Case 1: R = Q */
        /* Simultaneously bidiagonalize X11 and X21 */
        dorbdb1_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, & theta[1], &work[iphi], &work[itaup1], &work[itaup2], &work[ itauq1], &work[iorbdb], &lorbdb, &childinfo);
        /* Accumulate Householder reflectors */
        if (wantu1 && *p > 0)
        {
            dlacpy_("L", p, q, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
            dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqr, &childinfo);
        }
        if (wantu2 && *m - *p > 0)
        {
            i__1 = *m - *p;
            dlacpy_("L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
            i__1 = *m - *p;
            i__2 = *m - *p;
            dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], & work[iorgqr], &lorgqr, &childinfo);
        }
        if (wantv1t && *q > 0)
        {
            v1t[v1t_dim1 + 1] = 1.;
            i__1 = *q;
            for (j = 2;
                    j <= i__1;
                    ++j)
            {
                v1t[j * v1t_dim1 + 1] = 0.;
                v1t[j + v1t_dim1] = 0.;
            }
            i__1 = *q - 1;
            i__2 = *q - 1;
            dlacpy_("U", &i__1, &i__2, &x21[(x21_dim1 << 1) + 1], ldx21, &v1t[ (v1t_dim1 << 1) + 2], ldv1t);
            i__1 = *q - 1;
            i__2 = *q - 1;
            i__3 = *q - 1;
            dorglq_fla(&i__1, &i__2, &i__3, &v1t[(v1t_dim1 << 1) + 2], ldv1t, & work[itauq1], &work[iorglq], &lorglq, &childinfo);
        }
        /* Simultaneously diagonalize X11 and X21. */
        dbbcsd_(jobu1, jobu2, jobv1t, "N", "N", m, p, q, &theta[1], &work[ iphi], &u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &v1t[ v1t_offset], ldv1t, &c__0, &c__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 > 0 && wantu2)
        {
            i__1 = *q;
            for (i__ = 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = *m - *p - *q + i__;
            }
            i__1 = *m - *p;
            for (i__ = *q + 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = i__ - *q;
            }
            i__1 = *m - *p;
            i__2 = *m - *p;
            dlapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
        }
    }
    else if (r__ == *p)
    {
        /* Case 2: R = P */
        /* Simultaneously bidiagonalize X11 and X21 */
        dorbdb2_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, & theta[1], &work[iphi], &work[itaup1], &work[itaup2], &work[ itauq1], &work[iorbdb], &lorbdb, &childinfo);
        /* Accumulate Householder reflectors */
        if (wantu1 && *p > 0)
        {
            u1[u1_dim1 + 1] = 1.;
            i__1 = *p;
            for (j = 2;
                    j <= i__1;
                    ++j)
            {
                u1[j * u1_dim1 + 1] = 0.;
                u1[j + u1_dim1] = 0.;
            }
            i__1 = *p - 1;
            i__2 = *p - 1;
            dlacpy_("L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &u1[( u1_dim1 << 1) + 2], ldu1);
            i__1 = *p - 1;
            i__2 = *p - 1;
            i__3 = *p - 1;
            dorgqr_fla(&i__1, &i__2, &i__3, &u1[(u1_dim1 << 1) + 2], ldu1, &work[ itaup1], &work[iorgqr], &lorgqr, &childinfo);
        }
        if (wantu2 && *m - *p > 0)
        {
            i__1 = *m - *p;
            dlacpy_("L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
            i__1 = *m - *p;
            i__2 = *m - *p;
            dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], & work[iorgqr], &lorgqr, &childinfo);
        }
        if (wantv1t && *q > 0)
        {
            dlacpy_("U", p, q, &x11[x11_offset], ldx11, &v1t[v1t_offset], ldv1t);
            dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, &work[itauq1], &work[ iorglq], &lorglq, &childinfo);
        }
        /* Simultaneously diagonalize X11 and X21. */
        dbbcsd_(jobv1t, "N", jobu1, jobu2, "T", m, q, p, &theta[1], &work[ iphi], &v1t[v1t_offset], ldv1t, &c__0, &c__1, &u1[u1_offset], ldu1, &u2[u2_offset], 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 > 0 && wantu2)
        {
            i__1 = *q;
            for (i__ = 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = *m - *p - *q + i__;
            }
            i__1 = *m - *p;
            for (i__ = *q + 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = i__ - *q;
            }
            i__1 = *m - *p;
            i__2 = *m - *p;
            dlapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
        }
    }
    else if (r__ == *m - *p)
    {
        /* Case 3: R = M-P */
        /* Simultaneously bidiagonalize X11 and X21 */
        dorbdb3_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, & theta[1], &work[iphi], &work[itaup1], &work[itaup2], &work[ itauq1], &work[iorbdb], &lorbdb, &childinfo);
        /* Accumulate Householder reflectors */
        if (wantu1 && *p > 0)
        {
            dlacpy_("L", p, q, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
            dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqr, &childinfo);
        }
        if (wantu2 && *m - *p > 0)
        {
            u2[u2_dim1 + 1] = 1.;
            i__1 = *m - *p;
            for (j = 2;
                    j <= i__1;
                    ++j)
            {
                u2[j * u2_dim1 + 1] = 0.;
                u2[j + u2_dim1] = 0.;
            }
            i__1 = *m - *p - 1;
            i__2 = *m - *p - 1;
            dlacpy_("L", &i__1, &i__2, &x21[x21_dim1 + 2], ldx21, &u2[( u2_dim1 << 1) + 2], ldu2);
            i__1 = *m - *p - 1;
            i__2 = *m - *p - 1;
            i__3 = *m - *p - 1;
            dorgqr_fla(&i__1, &i__2, &i__3, &u2[(u2_dim1 << 1) + 2], ldu2, &work[ itaup2], &work[iorgqr], &lorgqr, &childinfo);
        }
        if (wantv1t && *q > 0)
        {
            i__1 = *m - *p;
            dlacpy_("U", &i__1, q, &x21[x21_offset], ldx21, &v1t[v1t_offset], ldv1t);
            dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, &work[itauq1], &work[ iorglq], &lorglq, &childinfo);
        }
        /* Simultaneously diagonalize X11 and X21. */
        i__1 = *m - *q;
        i__2 = *m - *p;
        dbbcsd_("N", jobv1t, jobu2, jobu1, "T", m, &i__1, &i__2, &theta[1], & work[iphi], &c__0, &c__1, &v1t[v1t_offset], ldv1t, &u2[ u2_offset], ldu2, &u1[u1_offset], 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 > r__)
        {
            i__1 = r__;
            for (i__ = 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = *q - r__ + i__;
            }
            i__1 = *q;
            for (i__ = r__ + 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = i__ - r__;
            }
            if (wantu1)
            {
                dlapmt_(&c_false, p, q, &u1[u1_offset], ldu1, &iwork[1]);
            }
            if (wantv1t)
            {
                dlapmr_(&c_false, q, q, &v1t[v1t_offset], ldv1t, &iwork[1]);
            }
        }
    }
    else
    {
        /* Case 4: R = M-Q */
        /* Simultaneously bidiagonalize X11 and X21 */
        i__1 = lorbdb - *m;
        dorbdb4_(m, p, q, &x11[x11_offset], ldx11, &x21[x21_offset], ldx21, & theta[1], &work[iphi], &work[itaup1], &work[itaup2], &work[ itauq1], &work[iorbdb], &work[iorbdb + *m], &i__1, &childinfo) ;
        /* Accumulate Householder reflectors */
        if (wantu1 && *p > 0)
        {
            dcopy_(p, &work[iorbdb], &c__1, &u1[u1_offset], &c__1);
            i__1 = *p;
            for (j = 2;
                    j <= i__1;
                    ++j)
            {
                u1[j * u1_dim1 + 1] = 0.;
            }
            i__1 = *p - 1;
            i__2 = *m - *q - 1;
            dlacpy_("L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &u1[( u1_dim1 << 1) + 2], ldu1);
            i__1 = *m - *q;
            dorgqr_fla(p, p, &i__1, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqr, &childinfo);
        }
        if (wantu2 && *m - *p > 0)
        {
            i__1 = *m - *p;
            dcopy_(&i__1, &work[iorbdb + *p], &c__1, &u2[u2_offset], &c__1);
            i__1 = *m - *p;
            for (j = 2;
                    j <= i__1;
                    ++j)
            {
                u2[j * u2_dim1 + 1] = 0.;
            }
            i__1 = *m - *p - 1;
            i__2 = *m - *q - 1;
            dlacpy_("L", &i__1, &i__2, &x21[x21_dim1 + 2], ldx21, &u2[( u2_dim1 << 1) + 2], ldu2);
            i__1 = *m - *p;
            i__2 = *m - *p;
            i__3 = *m - *q;
            dorgqr_fla(&i__1, &i__2, &i__3, &u2[u2_offset], ldu2, &work[itaup2], &work[iorgqr], &lorgqr, &childinfo);
        }
        if (wantv1t && *q > 0)
        {
            i__1 = *m - *q;
            dlacpy_("U", &i__1, q, &x21[x21_offset], ldx21, &v1t[v1t_offset], ldv1t);
            i__1 = *p - (*m - *q);
            i__2 = *q - (*m - *q);
            dlacpy_("U", &i__1, &i__2, &x11[*m - *q + 1 + (*m - *q + 1) * x11_dim1], ldx11, &v1t[*m - *q + 1 + (*m - *q + 1) * v1t_dim1], ldv1t);
            i__1 = -(*p) + *q;
            i__2 = *q - *p;
            dlacpy_("U", &i__1, &i__2, &x21[*m - *q + 1 + (*p + 1) * x21_dim1] , ldx21, &v1t[*p + 1 + (*p + 1) * v1t_dim1], ldv1t);
            dorglq_fla(q, q, q, &v1t[v1t_offset], ldv1t, &work[itauq1], &work[ iorglq], &lorglq, &childinfo);
        }
        /* Simultaneously diagonalize X11 and X21. */
        i__1 = *m - *p;
        i__2 = *m - *q;
        dbbcsd_(jobu2, jobu1, "N", jobv1t, "N", m, &i__1, &i__2, &theta[1], & work[iphi], &u2[u2_offset], ldu2, &u1[u1_offset], ldu1, &c__0, &c__1, &v1t[v1t_offset], 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 > r__)
        {
            i__1 = r__;
            for (i__ = 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = *p - r__ + i__;
            }
            i__1 = *p;
            for (i__ = r__ + 1;
                    i__ <= i__1;
                    ++i__)
            {
                iwork[i__] = i__ - r__;
            }
            if (wantu1)
            {
                dlapmt_(&c_false, p, p, &u1[u1_offset], ldu1, &iwork[1]);
            }
            if (wantv1t)
            {
                dlapmr_(&c_false, p, q, &v1t[v1t_offset], ldv1t, &iwork[1]);
            }
        }
    }
    return 0;
    /* End of DORCSD2BY1 */
}