sorcsd.c File Reference

(r)

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

Function Documentation

sorcsd_()

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

{
    /* 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 *);
    real dummy[1];
    integer lbbcsdworkmin, itaup1, itaup2, itauq1, itauq2, lorbdbworkmin, lbbcsdworkopt;
    logical wantu1, wantu2;
    integer ibbcsd, lorbdbworkopt;
    extern /* Subroutine */
    int sbbcsd_(char *, char *, char *, char *, char * , integer *, integer *, integer *, real *, real *, real *, integer *, real *, integer *, real *, integer *, real *, integer * , real *, real *, real *, real *, real *, real *, real *, real *, real *, integer *, integer *);
    integer iorbdb, lorglqworkmin, lorgqrworkmin;
    extern /* Subroutine */
    int sorbdb_(char *, char *, integer *, integer *, integer *, real *, integer *, real *, integer *, real *, integer * , real *, integer *, real *, real *, real *, real *, real *, real *, real *, integer *, integer *), xerbla_(char *, integer *);
    integer lorglqworkopt, lorgqrworkopt;
    extern /* Subroutine */
    int slacpy_(char *, integer *, integer *, real *, integer *, real *, integer *);
    integer iorglq;
    extern /* Subroutine */
    int slapmr_(logical *, integer *, integer *, real *, integer *, integer *), slapmt_(logical *, integer *, integer *, real *, integer *, integer *);
    integer iorgqr;
    char signst[1];
    extern /* Subroutine */
    int sorglq_fla(integer *, integer *, integer *, real *, integer *, real *, real *, integer *, integer *);
    char transt[1];
    extern /* Subroutine */
    int sorgqr_fla(integer *, integer *, integer *, real *, integer *, real *, real *, integer *, integer *);
    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 Arrays .. */
    /* .. */
    /* .. 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';
        }
        sorcsd_(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;
        sorcsd_(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);
        sorgqr_fla(&i__1, &i__2, &i__3, dummy, &i__4, dummy, &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);
        sorglq_fla(&i__1, &i__2, &i__3, dummy, &i__4, dummy, &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);
        sorbdb_(trans, signs, m, p, q, &x11[x11_offset], ldx11, &x12[ x12_offset], ldx12, &x21[x21_offset], ldx21, &x22[x22_offset], ldx22, dummy, dummy, dummy, dummy, dummy, dummy, &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);
        sbbcsd_(jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, dummy, dummy, & u1[u1_offset], ldu1, &u2[u2_offset], ldu2, &v1t[v1t_offset], ldv1t, &v2t[v2t_offset], ldv2t, dummy, dummy, dummy, dummy, dummy, dummy, dummy, dummy, &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] = (real) 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_("SORCSD", &i__1);
        return 0;
    }
    else if (lquery)
    {
        return 0;
    }
    /* Transform to bidiagonal block form */
    sorbdb_(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)
        {
            slacpy_("L", p, q, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
            sorgqr_fla(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqrwork, info);
        }
        if (wantu2 && *m - *p > 0)
        {
            i__1 = *m - *p;
            slacpy_("L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
            i__1 = *m - *p;
            i__2 = *m - *p;
            sorgqr_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;
            slacpy_("U", &i__1, &i__2, &x11[(x11_dim1 << 1) + 1], ldx11, &v1t[ (v1t_dim1 << 1) + 2], ldv1t);
            v1t[v1t_dim1 + 1] = 1.f;
            i__1 = *q;
            for (j = 2;
                    j <= i__1;
                    ++j)
            {
                v1t[j * v1t_dim1 + 1] = 0.f;
                v1t[j + v1t_dim1] = 0.f;
            }
            i__1 = *q - 1;
            i__2 = *q - 1;
            i__3 = *q - 1;
            sorglq_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;
            slacpy_("U", p, &i__1, &x12[x12_offset], ldx12, &v2t[v2t_offset], ldv2t);
            i__1 = *m - *p - *q;
            i__2 = *m - *p - *q;
            slacpy_("U", &i__1, &i__2, &x22[*q + 1 + (*p + 1) * x22_dim1], ldx22, &v2t[*p + 1 + (*p + 1) * v2t_dim1], ldv2t);
            i__1 = *m - *q;
            i__2 = *m - *q;
            i__3 = *m - *q;
            sorglq_fla(&i__1, &i__2, &i__3, &v2t[v2t_offset], ldv2t, &work[ itauq2], &work[iorglq], &lorglqwork, info);
        }
    }
    else
    {
        if (wantu1 && *p > 0)
        {
            slacpy_("U", q, p, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
            sorglq_fla(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorglq], &lorglqwork, info);
        }
        if (wantu2 && *m - *p > 0)
        {
            i__1 = *m - *p;
            slacpy_("U", q, &i__1, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
            i__1 = *m - *p;
            i__2 = *m - *p;
            sorglq_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;
            slacpy_("L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &v1t[( v1t_dim1 << 1) + 2], ldv1t);
            v1t[v1t_dim1 + 1] = 1.f;
            i__1 = *q;
            for (j = 2;
                    j <= i__1;
                    ++j)
            {
                v1t[j * v1t_dim1 + 1] = 0.f;
                v1t[j + v1t_dim1] = 0.f;
            }
            i__1 = *q - 1;
            i__2 = *q - 1;
            i__3 = *q - 1;
            sorgqr_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;
            slacpy_("L", &i__1, p, &x12[x12_offset], ldx12, &v2t[v2t_offset], ldv2t);
            i__1 = *m - *p - *q;
            i__2 = *m - *p - *q;
            slacpy_("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;
            sorgqr_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 */
    sbbcsd_(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;
            slapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
        }
        else
        {
            i__1 = *m - *p;
            i__2 = *m - *p;
            slapmr_(&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;
            slapmt_(&c_false, &i__1, &i__2, &v2t[v2t_offset], ldv2t, &iwork[1] );
        }
        else
        {
            i__1 = *m - *q;
            i__2 = *m - *q;
            slapmr_(&c_false, &i__1, &i__2, &v2t[v2t_offset], ldv2t, &iwork[1] );
        }
    }
    return 0;
    /* End SORCSD */
}

References sorcsd_(), sorglq_fla(), and sorgqr_fla().

Referenced by sorcsd_().