236 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;
238 integer lworkmin, lworkopt, i__, j, r__, childinfo, lorglqmin, lorgqrmin, lorglqopt, lorgqropt, ib11d, ib11e, ib12d, ib12e, ib21d, ib21e, ib22d, ib22e, iphi;
239 extern logical lsame_(
char *,
char *);
242 integer itaup1, itaup2, itauq1;
246 integer ibbcsd, lbbcsd, iorbdb, lorbdb;
248 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 *);
253 integer lorglq, iorgqr, lorgqr;
255 int dorbdb1_(), dorbdb2_(), dorbdb3_(), dorbdb4_() ;
280 x11_offset = 1 + x11_dim1;
283 x21_offset = 1 + x21_dim1;
287 u1_offset = 1 + u1_dim1;
290 u2_offset = 1 + u2_dim1;
293 v1t_offset = 1 + v1t_dim1;
299 wantu1 = lsame_(jobu1,
"Y");
300 wantu2 = lsame_(jobu2,
"Y");
301 wantv1t = lsame_(jobv1t,
"Y");
302 lquery = *lwork == -1;
307 else if (*p < 0 || *p > *m)
311 else if (*q < 0 || *q > *m)
315 else if (*ldx11 < max(1,*p))
324 if (*ldx21 < max(i__1,i__2))
328 else if (wantu1 && *ldu1 < *p)
332 else if (wantu2 && *ldu2 < *m - *p)
336 else if (wantv1t && *ldv1t < *q)
342 i__1 = *p, i__2 = *m - *p, i__1 = min(i__1,i__2);
345 r__ = min(i__1,i__2);
369 ib11d = iphi + max(i__1,i__2);
370 ib11e = ib11d + max(1,r__);
374 ib12d = ib11e + max(i__1,i__2);
375 ib12e = ib12d + max(1,r__);
379 ib21d = ib12e + max(i__1,i__2);
380 ib21e = ib21d + max(1,r__);
384 ib22d = ib21e + max(i__1,i__2);
385 ib22e = ib22d + max(1,r__);
389 ibbcsd = ib22e + max(i__1,i__2);
393 itaup1 = iphi + max(i__1,i__2);
394 itaup2 = itaup1 + max(1,*p);
398 itauq1 = itaup2 + max(i__1,i__2);
399 iorbdb = itauq1 + max(1,*q);
400 iorgqr = itauq1 + max(1,*q);
401 iorglq = itauq1 + max(1,*q);
404 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);
409 lorgqrmin = max(1,*p);
416 dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, (
doublereal*)&c__0, &work[1] , &c_n1, &childinfo);
420 lorgqrmin = max(i__1,i__2);
426 i__1 = max(i__2,i__3);
430 i__4 = max(i__5,i__6);
434 i__7 = max(i__8,i__9);
435 dorglq_fla(&i__1, &i__4, &i__7, &v1t[v1t_offset], ldv1t, (
doublereal*)&c__0, & work[1], &c_n1, &childinfo);
439 lorglqmin = max(i__1,i__2);
441 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);
446 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);
448 if (*p - 1 >= *m - *p)
453 dorgqr_fla(&i__1, &i__2, &i__3, &u1[(u1_dim1 << 1) + 2], ldu1, (
doublereal*)&c__0, &work[1], &c_n1, &childinfo);
457 lorgqrmin = max(i__1,i__2);
464 dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, (
doublereal*)&c__0, &work[1] , &c_n1, &childinfo);
468 lorgqrmin = max(i__1,i__2);
471 dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, (
doublereal*)&c__0, &work[1], & c_n1, &childinfo);
472 lorglqmin = max(1,*q);
474 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);
477 else if (r__ == *m - *p)
479 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);
481 if (*p >= *m - *p - 1)
484 lorgqrmin = max(1,*p);
492 dorgqr_fla(&i__1, &i__2, &i__3, &u2[(u2_dim1 << 1) + 2], ldu2, (
doublereal*)&c__0, &work[1], &c_n1, &childinfo);
496 lorgqrmin = max(i__1,i__2);
499 dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, (
doublereal*)&c__0, &work[1], & c_n1, &childinfo);
500 lorglqmin = max(1,*q);
504 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);
509 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);
510 lorbdb = *m + (
integer) work[1];
514 dorgqr_fla(p, p, &i__1, &u1[u1_offset], ldu1, (
doublereal*)&c__0, &work[1], & c_n1, &childinfo);
515 lorgqrmin = max(1,*p);
523 dorgqr_fla(&i__1, &i__2, &i__3, &u2[u2_offset], ldu2, (
doublereal*)&c__0, & work[1], &c_n1, &childinfo);
527 lorgqrmin = max(i__1,i__2);
531 lorglqmin = max(1,*q);
535 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);
539 i__1 = iorbdb + lorbdb - 1, i__2 = iorgqr + lorgqrmin - 1, i__1 = max( i__1,i__2), i__2 = iorglq + lorglqmin - 1;
540 i__1 = max(i__1, i__2);
541 i__2 = ibbcsd + lbbcsd - 1;
542 lworkmin = max(i__1,i__2);
544 i__1 = iorbdb + lorbdb - 1, i__2 = iorgqr + lorgqropt - 1, i__1 = max( i__1,i__2), i__2 = iorglq + lorglqopt - 1;
545 i__1 = max(i__1, i__2);
546 i__2 = ibbcsd + lbbcsd - 1;
547 lworkopt = max(i__1,i__2);
549 if (*lwork < lworkmin && ! lquery)
557 xerbla_(
"DORCSD2BY1", &i__1);
564 lorgqr = *lwork - iorgqr + 1;
565 lorglq = *lwork - iorglq + 1;
572 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);
574 if (wantu1 && *p > 0)
576 dlacpy_(
"L", p, q, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
577 dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqr, &childinfo);
579 if (wantu2 && *m - *p > 0)
582 dlacpy_(
"L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
585 dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], & work[iorgqr], &lorgqr, &childinfo);
587 if (wantv1t && *q > 0)
589 v1t[v1t_dim1 + 1] = 1.;
595 v1t[j * v1t_dim1 + 1] = 0.;
596 v1t[j + v1t_dim1] = 0.;
600 dlacpy_(
"U", &i__1, &i__2, &x21[(x21_dim1 << 1) + 1], ldx21, &v1t[ (v1t_dim1 << 1) + 2], ldv1t);
604 dorglq_fla(&i__1, &i__2, &i__3, &v1t[(v1t_dim1 << 1) + 2], ldv1t, & work[itauq1], &work[iorglq], &lorglq, &childinfo);
607 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);
610 if (*q > 0 && wantu2)
617 iwork[i__] = *m - *p - *q + i__;
624 iwork[i__] = i__ - *q;
628 dlapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
635 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);
637 if (wantu1 && *p > 0)
639 u1[u1_dim1 + 1] = 1.;
645 u1[j * u1_dim1 + 1] = 0.;
646 u1[j + u1_dim1] = 0.;
650 dlacpy_(
"L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &u1[( u1_dim1 << 1) + 2], ldu1);
654 dorgqr_fla(&i__1, &i__2, &i__3, &u1[(u1_dim1 << 1) + 2], ldu1, &work[ itaup1], &work[iorgqr], &lorgqr, &childinfo);
656 if (wantu2 && *m - *p > 0)
659 dlacpy_(
"L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
662 dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], & work[iorgqr], &lorgqr, &childinfo);
664 if (wantv1t && *q > 0)
666 dlacpy_(
"U", p, q, &x11[x11_offset], ldx11, &v1t[v1t_offset], ldv1t);
667 dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, &work[itauq1], &work[ iorglq], &lorglq, &childinfo);
670 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);
673 if (*q > 0 && wantu2)
680 iwork[i__] = *m - *p - *q + i__;
687 iwork[i__] = i__ - *q;
691 dlapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
694 else if (r__ == *m - *p)
698 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);
700 if (wantu1 && *p > 0)
702 dlacpy_(
"L", p, q, &x11[x11_offset], ldx11, &u1[u1_offset], ldu1);
703 dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqr, &childinfo);
705 if (wantu2 && *m - *p > 0)
707 u2[u2_dim1 + 1] = 1.;
713 u2[j * u2_dim1 + 1] = 0.;
714 u2[j + u2_dim1] = 0.;
718 dlacpy_(
"L", &i__1, &i__2, &x21[x21_dim1 + 2], ldx21, &u2[( u2_dim1 << 1) + 2], ldu2);
722 dorgqr_fla(&i__1, &i__2, &i__3, &u2[(u2_dim1 << 1) + 2], ldu2, &work[ itaup2], &work[iorgqr], &lorgqr, &childinfo);
724 if (wantv1t && *q > 0)
727 dlacpy_(
"U", &i__1, q, &x21[x21_offset], ldx21, &v1t[v1t_offset], ldv1t);
728 dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, &work[itauq1], &work[ iorglq], &lorglq, &childinfo);
733 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);
743 iwork[i__] = *q - r__ + i__;
750 iwork[i__] = i__ - r__;
754 dlapmt_(&c_false, p, q, &u1[u1_offset], ldu1, &iwork[1]);
758 dlapmr_(&c_false, q, q, &v1t[v1t_offset], ldv1t, &iwork[1]);
767 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) ;
769 if (wantu1 && *p > 0)
771 dcopy_(p, &work[iorbdb], &c__1, &u1[u1_offset], &c__1);
777 u1[j * u1_dim1 + 1] = 0.;
781 dlacpy_(
"L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &u1[( u1_dim1 << 1) + 2], ldu1);
783 dorgqr_fla(p, p, &i__1, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqr, &childinfo);
785 if (wantu2 && *m - *p > 0)
788 dcopy_(&i__1, &work[iorbdb + *p], &c__1, &u2[u2_offset], &c__1);
794 u2[j * u2_dim1 + 1] = 0.;
798 dlacpy_(
"L", &i__1, &i__2, &x21[x21_dim1 + 2], ldx21, &u2[( u2_dim1 << 1) + 2], ldu2);
802 dorgqr_fla(&i__1, &i__2, &i__3, &u2[u2_offset], ldu2, &work[itaup2], &work[iorgqr], &lorgqr, &childinfo);
804 if (wantv1t && *q > 0)
807 dlacpy_(
"U", &i__1, q, &x21[x21_offset], ldx21, &v1t[v1t_offset], ldv1t);
808 i__1 = *p - (*m - *q);
809 i__2 = *q - (*m - *q);
810 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);
813 dlacpy_(
"U", &i__1, &i__2, &x21[*m - *q + 1 + (*p + 1) * x21_dim1] , ldx21, &v1t[*p + 1 + (*p + 1) * v1t_dim1], ldv1t);
814 dorglq_fla(q, q, q, &v1t[v1t_offset], ldv1t, &work[itauq1], &work[ iorglq], &lorglq, &childinfo);
819 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);
829 iwork[i__] = *p - r__ + i__;
836 iwork[i__] = i__ - r__;
840 dlapmt_(&c_false, p, p, &u1[u1_offset], ldu1, &iwork[1]);
844 dlapmr_(&c_false, p, q, &v1t[v1t_offset], ldv1t, &iwork[1]);
double doublereal
Definition: FLA_f2c.h:31
int integer
Definition: FLA_f2c.h:25
int logical
Definition: FLA_f2c.h:36
int dorglq_fla(integer *m, integer *n, integer *k, doublereal *a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, integer *info)
Definition: dorglq.c:122
int dorgqr_fla(integer *m, integer *n, integer *k, doublereal *a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, integer *info)
Definition: dorgqr.c:123