libflame  revision_anchor
Functions
dorcsd.c File Reference

(r)

Functions

int dorcsd_ (char *jobu1, char *jobu2, char *jobv1t, char *jobv2t, char *trans, char *signs, integer *m, integer *p, integer *q, doublereal *x11, integer *ldx11, doublereal *x12, integer *ldx12, doublereal *x21, integer *ldx21, doublereal *x22, integer *ldx22, doublereal *theta, doublereal *u1, integer *ldu1, doublereal *u2, integer *ldu2, doublereal *v1t, integer *ldv1t, doublereal *v2t, integer *ldv2t, doublereal *work, integer *lwork, integer *iwork, integer *info)

Function Documentation

int dorcsd_ ( char *  jobu1,
char *  jobu2,
char *  jobv1t,
char *  jobv2t,
char *  trans,
char *  signs,
integer m,
integer p,
integer q,
doublereal x11,
integer ldx11,
doublereal x12,
integer ldx12,
doublereal x21,
integer ldx21,
doublereal x22,
integer ldx22,
doublereal theta,
doublereal u1,
integer ldu1,
doublereal u2,
integer ldu2,
doublereal v1t,
integer ldv1t,
doublereal v2t,
integer ldv2t,
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, v2t_dim1, v2t_offset, x11_dim1, x11_offset, x12_dim1, x12_offset, x21_dim1, x21_offset, x22_dim1, x22_offset, i__1, i__2, i__3, i__4, i__5, i__6;
    /* Local variables */
    logical colmajor;
    integer lworkmin, lworkopt, i__, j, childinfo, lbbcsdwork, lorbdbwork, lorglqwork, lorgqrwork, ib11d, ib11e, ib12d, ib12e, ib21d, ib21e, ib22d, ib22e, iphi;
    logical defaultsigns;
    extern logical lsame_(char *, char *);
    integer lbbcsdworkmin, itaup1, itaup2, itauq1, itauq2, lorbdbworkmin, lbbcsdworkopt;
    logical wantu1, wantu2;
    extern /* Subroutine */
    int dbbcsd_(char *, char *, char *, char *, char * , integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *);
    integer ibbcsd, lorbdbworkopt;
    extern /* Subroutine */
    int dorbdb_(char *, char *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *);
    integer iorbdb, lorglqworkmin, lorgqrworkmin;
    extern /* Subroutine */
    int dlacpy_(char *, integer *, integer *, doublereal *, integer *, doublereal *, integer *), dlapmr_(logical *, integer *, integer *, doublereal *, integer *, integer *), xerbla_(char *, integer *), dlapmt_(logical *, integer *, integer *, doublereal *, integer *, integer *);
    integer lorglqworkopt;
    extern /* Subroutine */
    int dorglq_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *);
    integer lorgqrworkopt, iorglq;
    extern /* Subroutine */
    int dorgqr_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *);
    integer iorgqr;
    char signst[1], transt[1];
    logical lquery, wantv1t, wantv2t;
    /* -- LAPACK computational routine (version 3.5.0) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* November 2013 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* =================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. Intrinsic Functions */
    /* .. */
    /* .. Executable Statements .. */
    /* Test input arguments */
    /* Parameter adjustments */
    x11_dim1 = *ldx11;
    x11_offset = 1 + x11_dim1;
    x11 -= x11_offset;
    x12_dim1 = *ldx12;
    x12_offset = 1 + x12_dim1;
    x12 -= x12_offset;
    x21_dim1 = *ldx21;
    x21_offset = 1 + x21_dim1;
    x21 -= x21_offset;
    x22_dim1 = *ldx22;
    x22_offset = 1 + x22_dim1;
    x22 -= x22_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;
    v2t_dim1 = *ldv2t;
    v2t_offset = 1 + v2t_dim1;
    v2t -= v2t_offset;
    --work;
    --iwork;
    /* Function Body */
    *info = 0;
    wantu1 = lsame_(jobu1, "Y");
    wantu2 = lsame_(jobu2, "Y");
    wantv1t = lsame_(jobv1t, "Y");
    wantv2t = lsame_(jobv2t, "Y");
    colmajor = ! lsame_(trans, "T");
    defaultsigns = ! lsame_(signs, "O");
    lquery = *lwork == -1;
    if (*m < 0)
    {
        *info = -7;
    }
    else if (*p < 0 || *p > *m)
    {
        *info = -8;
    }
    else if (*q < 0 || *q > *m)
    {
        *info = -9;
    }
    else if (colmajor && *ldx11 < max(1,*p))
    {
        *info = -11;
    }
    else if (! colmajor && *ldx11 < max(1,*q))
    {
        *info = -11;
    }
    else if (colmajor && *ldx12 < max(1,*p))
    {
        *info = -13;
    }
    else /* if(complicated condition) */
    {
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *q; // , expr subst
        if (! colmajor && *ldx12 < max(i__1,i__2))
        {
            *info = -13;
        }
        else /* if(complicated condition) */
        {
            /* Computing MAX */
            i__1 = 1;
            i__2 = *m - *p; // , expr subst
            if (colmajor && *ldx21 < max(i__1,i__2))
            {
                *info = -15;
            }
            else if (! colmajor && *ldx21 < max(1,*q))
            {
                *info = -15;
            }
            else /* if(complicated condition) */
            {
                /* Computing MAX */
                i__1 = 1;
                i__2 = *m - *p; // , expr subst
                if (colmajor && *ldx22 < max(i__1,i__2))
                {
                    *info = -17;
                }
                else /* if(complicated condition) */
                {
                    /* Computing MAX */
                    i__1 = 1;
                    i__2 = *m - *q; // , expr subst
                    if (! colmajor && *ldx22 < max(i__1,i__2))
                    {
                        *info = -17;
                    }
                    else if (wantu1 && *ldu1 < *p)
                    {
                        *info = -20;
                    }
                    else if (wantu2 && *ldu2 < *m - *p)
                    {
                        *info = -22;
                    }
                    else if (wantv1t && *ldv1t < *q)
                    {
                        *info = -24;
                    }
                    else if (wantv2t && *ldv2t < *m - *q)
                    {
                        *info = -26;
                    }
                }
            }
        }
    }
    /* Work with transpose if convenient */
    /* Computing MIN */
    i__1 = *p;
    i__2 = *m - *p; // , expr subst
    /* Computing MIN */
    i__3 = *q;
    i__4 = *m - *q; // , expr subst
    if (*info == 0 && min(i__1,i__2) < min(i__3,i__4))
    {
        if (colmajor)
        {
            *(unsigned char *)transt = 'T';
        }
        else
        {
            *(unsigned char *)transt = 'N';
        }
        if (defaultsigns)
        {
            *(unsigned char *)signst = 'O';
        }
        else
        {
            *(unsigned char *)signst = 'D';
        }
        dorcsd_(jobv1t, jobv2t, jobu1, jobu2, transt, signst, m, q, p, &x11[ x11_offset], ldx11, &x21[x21_offset], ldx21, &x12[x12_offset], ldx12, &x22[x22_offset], ldx22, &theta[1], &v1t[v1t_offset], ldv1t, &v2t[v2t_offset], ldv2t, &u1[u1_offset], ldu1, &u2[ u2_offset], ldu2, &work[1], lwork, &iwork[1], info);
        return 0;
    }
    /* Work with permutation [ 0 I;
    I 0 ] * X * [ 0 I;
    I 0 ] if */
    /* convenient */
    if (*info == 0 && *m - *q < *q)
    {
        if (defaultsigns)
        {
            *(unsigned char *)signst = 'O';
        }
        else
        {
            *(unsigned char *)signst = 'D';
        }
        i__1 = *m - *p;
        i__2 = *m - *q;
        dorcsd_(jobu2, jobu1, jobv2t, jobv1t, trans, signst, m, &i__1, &i__2, &x22[x22_offset], ldx22, &x21[x21_offset], ldx21, &x12[ x12_offset], ldx12, &x11[x11_offset], ldx11, &theta[1], &u2[ u2_offset], ldu2, &u1[u1_offset], ldu1, &v2t[v2t_offset], ldv2t, &v1t[v1t_offset], ldv1t, &work[1], lwork, &iwork[1], info);
        return 0;
    }
    /* Compute workspace */
    if (*info == 0)
    {
        iphi = 2;
        /* Computing MAX */
        i__1 = 1;
        i__2 = *q - 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);
        itauq2 = itauq1 + max(1,*q);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *q; // , expr subst
        iorgqr = itauq2 + max(i__1,i__2);
        i__1 = *m - *q;
        i__2 = *m - *q;
        i__3 = *m - *q;
        /* Computing MAX */
        i__5 = 1;
        i__6 = *m - *q; // , expr subst
        i__4 = max(i__5,i__6);
        dorgqr_fla(&i__1, &i__2, &i__3, &u1[u1_offset], &i__4, &u1[u1_offset], & work[1], &c_n1, &childinfo);
        lorgqrworkopt = (integer) work[1];
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *q; // , expr subst
        lorgqrworkmin = max(i__1,i__2);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *q; // , expr subst
        iorglq = itauq2 + max(i__1,i__2);
        i__1 = *m - *q;
        i__2 = *m - *q;
        i__3 = *m - *q;
        /* Computing MAX */
        i__5 = 1;
        i__6 = *m - *q; // , expr subst
        i__4 = max(i__5,i__6);
        dorglq_fla(&i__1, &i__2, &i__3, &u1[u1_offset], &i__4, &u1[u1_offset], & work[1], &c_n1, &childinfo);
        lorglqworkopt = (integer) work[1];
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *q; // , expr subst
        lorglqworkmin = max(i__1,i__2);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *q; // , expr subst
        iorbdb = itauq2 + max(i__1,i__2);
        dorbdb_(trans, signs, m, p, q, &x11[x11_offset], ldx11, &x12[ x12_offset], ldx12, &x21[x21_offset], ldx21, &x22[x22_offset], ldx22, &theta[1], &v1t[v1t_offset], &u1[u1_offset], &u2[ u2_offset], &v1t[v1t_offset], &v2t[v2t_offset], &work[1], & c_n1, &childinfo);
        lorbdbworkopt = (integer) work[1];
        lorbdbworkmin = lorbdbworkopt;
        /* Computing MAX */
        i__1 = 1;
        i__2 = *m - *q; // , expr subst
        ib11d = itauq2 + max(i__1,i__2);
        ib11e = ib11d + max(1,*q);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *q - 1; // , expr subst
        ib12d = ib11e + max(i__1,i__2);
        ib12e = ib12d + max(1,*q);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *q - 1; // , expr subst
        ib21d = ib12e + max(i__1,i__2);
        ib21e = ib21d + max(1,*q);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *q - 1; // , expr subst
        ib22d = ib21e + max(i__1,i__2);
        ib22e = ib22d + max(1,*q);
        /* Computing MAX */
        i__1 = 1;
        i__2 = *q - 1; // , expr subst
        ibbcsd = ib22e + max(i__1,i__2);
        dbbcsd_(jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, &theta[1], & theta[1], &u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &v1t[ v1t_offset], ldv1t, &v2t[v2t_offset], ldv2t, &u1[u1_offset], & u1[u1_offset], &u1[u1_offset], &u1[u1_offset], &u1[u1_offset], &u1[u1_offset], &u1[u1_offset], &u1[u1_offset], &work[1], & c_n1, &childinfo);
        lbbcsdworkopt = (integer) work[1];
        lbbcsdworkmin = lbbcsdworkopt;
        /* Computing MAX */
        i__1 = iorgqr + lorgqrworkopt, i__2 = iorglq + lorglqworkopt, i__1 = max(i__1,i__2), i__2 = iorbdb + lorbdbworkopt;
        i__1 = max( i__1,i__2);
        i__2 = ibbcsd + lbbcsdworkopt; // ; expr subst
        lworkopt = max(i__1,i__2) - 1;
        /* Computing MAX */
        i__1 = iorgqr + lorgqrworkmin, i__2 = iorglq + lorglqworkmin, i__1 = max(i__1,i__2), i__2 = iorbdb + lorbdbworkopt;
        i__1 = max( i__1,i__2);
        i__2 = ibbcsd + lbbcsdworkmin; // ; expr subst
        lworkmin = max(i__1,i__2) - 1;
        work[1] = (doublereal) max(lworkopt,lworkmin);
        if (*lwork < lworkmin && ! lquery)
        {
            *info = -22;
        }
        else
        {
            lorgqrwork = *lwork - iorgqr + 1;
            lorglqwork = *lwork - iorglq + 1;
            lorbdbwork = *lwork - iorbdb + 1;
            lbbcsdwork = *lwork - ibbcsd + 1;
        }
    }
    /* Abort if any illegal arguments */
    if (*info != 0)
    {
        i__1 = -(*info);
        xerbla_("DORCSD", &i__1);
        return 0;
    }
    else if (lquery)
    {
        return 0;
    }
    /* Transform to bidiagonal block form */
    dorbdb_(trans, signs, m, p, q, &x11[x11_offset], ldx11, &x12[x12_offset], ldx12, &x21[x21_offset], ldx21, &x22[x22_offset], ldx22, &theta[1] , &work[iphi], &work[itaup1], &work[itaup2], &work[itauq1], &work[ itauq2], &work[iorbdb], &lorbdbwork, &childinfo);
    /* Accumulate Householder reflectors */
    if (colmajor)
    {
        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], &lorgqrwork, info);
        }
        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], &lorgqrwork, info);
        }
        if (wantv1t && *q > 0)
        {
            i__1 = *q - 1;
            i__2 = *q - 1;
            dlacpy_("U", &i__1, &i__2, &x11[(x11_dim1 << 1) + 1], ldx11, &v1t[ (v1t_dim1 << 1) + 2], ldv1t);
            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;
            i__3 = *q - 1;
            dorglq_fla(&i__1, &i__2, &i__3, &v1t[(v1t_dim1 << 1) + 2], ldv1t, & work[itauq1], &work[iorglq], &lorglqwork, info);
        }
        if (wantv2t && *m - *q > 0)
        {
            i__1 = *m - *q;
            dlacpy_("U", p, &i__1, &x12[x12_offset], ldx12, &v2t[v2t_offset], ldv2t);
            if (*m - *p > *q)
            {
                i__1 = *m - *p - *q;
                i__2 = *m - *p - *q;
                dlacpy_("U", &i__1, &i__2, &x22[*q + 1 + (*p + 1) * x22_dim1], ldx22, &v2t[*p + 1 + (*p + 1) * v2t_dim1], ldv2t);
            }
            if (*m > *q)
            {
                i__1 = *m - *q;
                i__2 = *m - *q;
                i__3 = *m - *q;
                dorglq_fla(&i__1, &i__2, &i__3, &v2t[v2t_offset], ldv2t, &work[ itauq2], &work[iorglq], &lorglqwork, info);
            }
        }
    }
    else
    {
        if (wantu1 && *p > 0)
        {
            dlacpy_("U", q, p, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
            dorglq_fla(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorglq], &lorglqwork, info);
        }
        if (wantu2 && *m - *p > 0)
        {
            i__1 = *m - *p;
            dlacpy_("U", q, &i__1, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
            i__1 = *m - *p;
            i__2 = *m - *p;
            dorglq_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], & work[iorglq], &lorglqwork, info);
        }
        if (wantv1t && *q > 0)
        {
            i__1 = *q - 1;
            i__2 = *q - 1;
            dlacpy_("L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &v1t[( v1t_dim1 << 1) + 2], ldv1t);
            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;
            i__3 = *q - 1;
            dorgqr_fla(&i__1, &i__2, &i__3, &v1t[(v1t_dim1 << 1) + 2], ldv1t, & work[itauq1], &work[iorgqr], &lorgqrwork, info);
        }
        if (wantv2t && *m - *q > 0)
        {
            i__1 = *m - *q;
            dlacpy_("L", &i__1, p, &x12[x12_offset], ldx12, &v2t[v2t_offset], ldv2t);
            i__1 = *m - *p - *q;
            i__2 = *m - *p - *q;
            dlacpy_("L", &i__1, &i__2, &x22[*p + 1 + (*q + 1) * x22_dim1], ldx22, &v2t[*p + 1 + (*p + 1) * v2t_dim1], ldv2t);
            i__1 = *m - *q;
            i__2 = *m - *q;
            i__3 = *m - *q;
            dorgqr_fla(&i__1, &i__2, &i__3, &v2t[v2t_offset], ldv2t, &work[ itauq2], &work[iorgqr], &lorgqrwork, info);
        }
    }
    /* Compute the CSD of the matrix in bidiagonal-block form */
    dbbcsd_(jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, &theta[1], &work[ iphi], &u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &v1t[ v1t_offset], ldv1t, &v2t[v2t_offset], ldv2t, &work[ib11d], &work[ ib11e], &work[ib12d], &work[ib12e], &work[ib21d], &work[ib21e], & work[ib22d], &work[ib22e], &work[ibbcsd], &lbbcsdwork, info);
    /* Permute rows and columns to place identity submatrices in top- */
    /* left corner of (1,1)-block and/or bottom-right corner of (1,2)- */
    /* block and/or bottom-right corner of (2,1)-block and/or top-left */
    /* corner of (2,2)-block */
    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;
        }
        if (colmajor)
        {
            i__1 = *m - *p;
            i__2 = *m - *p;
            dlapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
        }
        else
        {
            i__1 = *m - *p;
            i__2 = *m - *p;
            dlapmr_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
        }
    }
    if (*m > 0 && wantv2t)
    {
        i__1 = *p;
        for (i__ = 1;
                i__ <= i__1;
                ++i__)
        {
            iwork[i__] = *m - *p - *q + i__;
        }
        i__1 = *m - *q;
        for (i__ = *p + 1;
                i__ <= i__1;
                ++i__)
        {
            iwork[i__] = i__ - *p;
        }
        if (! colmajor)
        {
            i__1 = *m - *q;
            i__2 = *m - *q;
            dlapmt_(&c_false, &i__1, &i__2, &v2t[v2t_offset], ldv2t, &iwork[1] );
        }
        else
        {
            i__1 = *m - *q;
            i__2 = *m - *q;
            dlapmr_(&c_false, &i__1, &i__2, &v2t[v2t_offset], ldv2t, &iwork[1] );
        }
    }
    return 0;
    /* End DORCSD */
}