{
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;
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
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
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
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
int dorglq_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *);
integer lorgqrworkopt, iorglq;
extern
int dorgqr_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *);
integer iorgqr;
char signst[1], transt[1];
logical lquery, wantv1t, wantv2t;
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;
*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
{
i__1 = 1;
i__2 = *m - *q;
if (! colmajor && *ldx12 < max(i__1,i__2))
{
*info = -13;
}
else
{
i__1 = 1;
i__2 = *m - *p;
if (colmajor && *ldx21 < max(i__1,i__2))
{
*info = -15;
}
else if (! colmajor && *ldx21 < max(1,*q))
{
*info = -15;
}
else
{
i__1 = 1;
i__2 = *m - *p;
if (colmajor && *ldx22 < max(i__1,i__2))
{
*info = -17;
}
else
{
i__1 = 1;
i__2 = *m - *q;
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;
}
}
}
}
}
i__1 = *p;
i__2 = *m - *p;
i__3 = *q;
i__4 = *m - *q;
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;
}
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;
}
if (*info == 0)
{
iphi = 2;
i__1 = 1;
i__2 = *q - 1;
itaup1 = iphi + max(i__1,i__2);
itaup2 = itaup1 + max(1,*p);
i__1 = 1;
i__2 = *m - *p;
itauq1 = itaup2 + max(i__1,i__2);
itauq2 = itauq1 + max(1,*q);
i__1 = 1;
i__2 = *m - *q;
iorgqr = itauq2 + max(i__1,i__2);
i__1 = *m - *q;
i__2 = *m - *q;
i__3 = *m - *q;
i__5 = 1;
i__6 = *m - *q;
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];
i__1 = 1;
i__2 = *m - *q;
lorgqrworkmin = max(i__1,i__2);
i__1 = 1;
i__2 = *m - *q;
iorglq = itauq2 + max(i__1,i__2);
i__1 = *m - *q;
i__2 = *m - *q;
i__3 = *m - *q;
i__5 = 1;
i__6 = *m - *q;
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];
i__1 = 1;
i__2 = *m - *q;
lorglqworkmin = max(i__1,i__2);
i__1 = 1;
i__2 = *m - *q;
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;
i__1 = 1;
i__2 = *m - *q;
ib11d = itauq2 + max(i__1,i__2);
ib11e = ib11d + max(1,*q);
i__1 = 1;
i__2 = *q - 1;
ib12d = ib11e + max(i__1,i__2);
ib12e = ib12d + max(1,*q);
i__1 = 1;
i__2 = *q - 1;
ib21d = ib12e + max(i__1,i__2);
ib21e = ib21d + max(1,*q);
i__1 = 1;
i__2 = *q - 1;
ib22d = ib21e + max(i__1,i__2);
ib22e = ib22d + max(1,*q);
i__1 = 1;
i__2 = *q - 1;
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;
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;
lworkopt = max(i__1,i__2) - 1;
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;
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;
}
}
if (*info != 0)
{
i__1 = -(*info);
xerbla_("DORCSD", &i__1);
return 0;
}
else if (lquery)
{
return 0;
}
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);
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);
}
}
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);
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;
}