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

◆ dorcsd2by1_()

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().

234 {
235  /* System generated locals */
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;
237  /* Local variables */
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 *);
240  extern /* Subroutine */
241  int dcopy_(integer *, doublereal *, integer *, doublereal *, integer *);
242  integer itaup1, itaup2, itauq1;
243  logical wantu1, wantu2;
244  extern /* Subroutine */
245  int dbbcsd_();
246  integer ibbcsd, lbbcsd, iorbdb, lorbdb;
247  extern /* Subroutine */
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 *);
249  integer iorglq;
250  extern int /* Subroutine */
253  integer lorglq, iorgqr, lorgqr;
254  extern /* Subroutine */
255  int dorbdb1_(), dorbdb2_(), dorbdb3_(), dorbdb4_() ;
256  logical lquery, wantv1t;
257  /* -- LAPACK computational routine (3.5.0) -- */
258  /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
259  /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
260  /* July 2012 */
261  /* .. Scalar Arguments .. */
262  /* .. */
263  /* .. Array Arguments .. */
264  /* .. */
265  /* ===================================================================== */
266  /* .. Parameters .. */
267  /* .. */
268  /* .. Local Scalars .. */
269  /* .. */
270  /* .. External Subroutines .. */
271  /* .. */
272  /* .. External Functions .. */
273  /* .. */
274  /* .. Intrinsic Function .. */
275  /* .. */
276  /* .. Executable Statements .. */
277  /* Test input arguments */
278  /* Parameter adjustments */
279  x11_dim1 = *ldx11;
280  x11_offset = 1 + x11_dim1;
281  x11 -= x11_offset;
282  x21_dim1 = *ldx21;
283  x21_offset = 1 + x21_dim1;
284  x21 -= x21_offset;
285  --theta;
286  u1_dim1 = *ldu1;
287  u1_offset = 1 + u1_dim1;
288  u1 -= u1_offset;
289  u2_dim1 = *ldu2;
290  u2_offset = 1 + u2_dim1;
291  u2 -= u2_offset;
292  v1t_dim1 = *ldv1t;
293  v1t_offset = 1 + v1t_dim1;
294  v1t -= v1t_offset;
295  --work;
296  --iwork;
297  /* Function Body */
298  *info = 0;
299  wantu1 = lsame_(jobu1, "Y");
300  wantu2 = lsame_(jobu2, "Y");
301  wantv1t = lsame_(jobv1t, "Y");
302  lquery = *lwork == -1;
303  if (*m < 0)
304  {
305  *info = -4;
306  }
307  else if (*p < 0 || *p > *m)
308  {
309  *info = -5;
310  }
311  else if (*q < 0 || *q > *m)
312  {
313  *info = -6;
314  }
315  else if (*ldx11 < max(1,*p))
316  {
317  *info = -8;
318  }
319  else /* if(complicated condition) */
320  {
321  /* Computing MAX */
322  i__1 = 1;
323  i__2 = *m - *p; // , expr subst
324  if (*ldx21 < max(i__1,i__2))
325  {
326  *info = -10;
327  }
328  else if (wantu1 && *ldu1 < *p)
329  {
330  *info = -13;
331  }
332  else if (wantu2 && *ldu2 < *m - *p)
333  {
334  *info = -15;
335  }
336  else if (wantv1t && *ldv1t < *q)
337  {
338  *info = -17;
339  }
340  }
341  /* Computing MIN */
342  i__1 = *p, i__2 = *m - *p, i__1 = min(i__1,i__2);
343  i__1 = min(i__1,*q);
344  i__2 = *m - *q; // ; expr subst
345  r__ = min(i__1,i__2);
346  /* Compute workspace */
347  /* WORK layout: */
348  /* |-------------------------------------------------------| */
349  /* | LWORKOPT (1) | */
350  /* |-------------------------------------------------------| */
351  /* | PHI (MAX(1,R-1)) | */
352  /* |-------------------------------------------------------| */
353  /* | TAUP1 (MAX(1,P)) | B11D (R) | */
354  /* | TAUP2 (MAX(1,M-P)) | B11E (R-1) | */
355  /* | TAUQ1 (MAX(1,Q)) | B12D (R) | */
356  /* |-----------------------------------------| B12E (R-1) | */
357  /* | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R) | */
358  /* | | | | B21E (R-1) | */
359  /* | | | | B22D (R) | */
360  /* | | | | B22E (R-1) | */
361  /* | | | | DBBCSD WORK | */
362  /* |-------------------------------------------------------| */
363  if (*info == 0)
364  {
365  iphi = 2;
366  /* Computing MAX */
367  i__1 = 1;
368  i__2 = r__ - 1; // , expr subst
369  ib11d = iphi + max(i__1,i__2);
370  ib11e = ib11d + max(1,r__);
371  /* Computing MAX */
372  i__1 = 1;
373  i__2 = r__ - 1; // , expr subst
374  ib12d = ib11e + max(i__1,i__2);
375  ib12e = ib12d + max(1,r__);
376  /* Computing MAX */
377  i__1 = 1;
378  i__2 = r__ - 1; // , expr subst
379  ib21d = ib12e + max(i__1,i__2);
380  ib21e = ib21d + max(1,r__);
381  /* Computing MAX */
382  i__1 = 1;
383  i__2 = r__ - 1; // , expr subst
384  ib22d = ib21e + max(i__1,i__2);
385  ib22e = ib22d + max(1,r__);
386  /* Computing MAX */
387  i__1 = 1;
388  i__2 = r__ - 1; // , expr subst
389  ibbcsd = ib22e + max(i__1,i__2);
390  /* Computing MAX */
391  i__1 = 1;
392  i__2 = r__ - 1; // , expr subst
393  itaup1 = iphi + max(i__1,i__2);
394  itaup2 = itaup1 + max(1,*p);
395  /* Computing MAX */
396  i__1 = 1;
397  i__2 = *m - *p; // , expr subst
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);
402  if (r__ == *q)
403  {
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);
405  lorbdb = (integer) work[1];
406  if (*p >= *m - *p)
407  {
408  dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
409  lorgqrmin = max(1,*p);
410  lorgqropt = (integer) work[1];
411  }
412  else
413  {
414  i__1 = *m - *p;
415  i__2 = *m - *p;
416  dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, (doublereal*)&c__0, &work[1] , &c_n1, &childinfo);
417  /* Computing MAX */
418  i__1 = 1;
419  i__2 = *m - *p; // , expr subst
420  lorgqrmin = max(i__1,i__2);
421  lorgqropt = (integer) work[1];
422  }
423  /* Computing MAX */
424  i__2 = 0;
425  i__3 = *q - 1; // , expr subst
426  i__1 = max(i__2,i__3);
427  /* Computing MAX */
428  i__5 = 0;
429  i__6 = *q - 1; // , expr subst
430  i__4 = max(i__5,i__6);
431  /* Computing MAX */
432  i__8 = 0;
433  i__9 = *q - 1; // , expr subst
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);
436  /* Computing MAX */
437  i__1 = 1;
438  i__2 = *q - 1; // , expr subst
439  lorglqmin = max(i__1,i__2);
440  lorglqopt = (integer) work[1];
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);
442  lbbcsd = (integer) work[1];
443  }
444  else if (r__ == *p)
445  {
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);
447  lorbdb = (integer) work[1];
448  if (*p - 1 >= *m - *p)
449  {
450  i__1 = *p - 1;
451  i__2 = *p - 1;
452  i__3 = *p - 1;
453  dorgqr_fla(&i__1, &i__2, &i__3, &u1[(u1_dim1 << 1) + 2], ldu1, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
454  /* Computing MAX */
455  i__1 = 1;
456  i__2 = *p - 1; // , expr subst
457  lorgqrmin = max(i__1,i__2);
458  lorgqropt = (integer) work[1];
459  }
460  else
461  {
462  i__1 = *m - *p;
463  i__2 = *m - *p;
464  dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, (doublereal*)&c__0, &work[1] , &c_n1, &childinfo);
465  /* Computing MAX */
466  i__1 = 1;
467  i__2 = *m - *p; // , expr subst
468  lorgqrmin = max(i__1,i__2);
469  lorgqropt = (integer) work[1];
470  }
471  dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, (doublereal*)&c__0, &work[1], & c_n1, &childinfo);
472  lorglqmin = max(1,*q);
473  lorglqopt = (integer) work[1];
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);
475  lbbcsd = (integer) work[1];
476  }
477  else if (r__ == *m - *p)
478  {
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);
480  lorbdb = (integer) work[1];
481  if (*p >= *m - *p - 1)
482  {
483  dorgqr_fla(p, p, q, &u1[u1_offset], ldu1, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
484  lorgqrmin = max(1,*p);
485  lorgqropt = (integer) work[1];
486  }
487  else
488  {
489  i__1 = *m - *p - 1;
490  i__2 = *m - *p - 1;
491  i__3 = *m - *p - 1;
492  dorgqr_fla(&i__1, &i__2, &i__3, &u2[(u2_dim1 << 1) + 2], ldu2, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
493  /* Computing MAX */
494  i__1 = 1;
495  i__2 = *m - *p - 1; // , expr subst
496  lorgqrmin = max(i__1,i__2);
497  lorgqropt = (integer) work[1];
498  }
499  dorglq_fla(q, q, &r__, &v1t[v1t_offset], ldv1t, (doublereal*)&c__0, &work[1], & c_n1, &childinfo);
500  lorglqmin = max(1,*q);
501  lorglqopt = (integer) work[1];
502  i__1 = *m - *q;
503  i__2 = *m - *p;
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);
505  lbbcsd = (integer) work[1];
506  }
507  else
508  {
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];
511  if (*p >= *m - *p)
512  {
513  i__1 = *m - *q;
514  dorgqr_fla(p, p, &i__1, &u1[u1_offset], ldu1, (doublereal*)&c__0, &work[1], & c_n1, &childinfo);
515  lorgqrmin = max(1,*p);
516  lorgqropt = (integer) work[1];
517  }
518  else
519  {
520  i__1 = *m - *p;
521  i__2 = *m - *p;
522  i__3 = *m - *q;
523  dorgqr_fla(&i__1, &i__2, &i__3, &u2[u2_offset], ldu2, (doublereal*)&c__0, & work[1], &c_n1, &childinfo);
524  /* Computing MAX */
525  i__1 = 1;
526  i__2 = *m - *p; // , expr subst
527  lorgqrmin = max(i__1,i__2);
528  lorgqropt = (integer) work[1];
529  }
530  dorglq_fla(q, q, q, &v1t[v1t_offset], ldv1t, (doublereal*)&c__0, &work[1], &c_n1, &childinfo);
531  lorglqmin = max(1,*q);
532  lorglqopt = (integer) work[1];
533  i__1 = *m - *p;
534  i__2 = *m - *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);
536  lbbcsd = (integer) work[1];
537  }
538  /* Computing MAX */
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; // ; expr subst
542  lworkmin = max(i__1,i__2);
543  /* Computing MAX */
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; // ; expr subst
547  lworkopt = max(i__1,i__2);
548  work[1] = (doublereal) lworkopt;
549  if (*lwork < lworkmin && ! lquery)
550  {
551  *info = -19;
552  }
553  }
554  if (*info != 0)
555  {
556  i__1 = -(*info);
557  xerbla_("DORCSD2BY1", &i__1);
558  return 0;
559  }
560  else if (lquery)
561  {
562  return 0;
563  }
564  lorgqr = *lwork - iorgqr + 1;
565  lorglq = *lwork - iorglq + 1;
566  /* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q, */
567  /* in which R = MIN(P,M-P,Q,M-Q) */
568  if (r__ == *q)
569  {
570  /* Case 1: R = Q */
571  /* Simultaneously bidiagonalize X11 and X21 */
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);
573  /* Accumulate Householder reflectors */
574  if (wantu1 && *p > 0)
575  {
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);
578  }
579  if (wantu2 && *m - *p > 0)
580  {
581  i__1 = *m - *p;
582  dlacpy_("L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
583  i__1 = *m - *p;
584  i__2 = *m - *p;
585  dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], & work[iorgqr], &lorgqr, &childinfo);
586  }
587  if (wantv1t && *q > 0)
588  {
589  v1t[v1t_dim1 + 1] = 1.;
590  i__1 = *q;
591  for (j = 2;
592  j <= i__1;
593  ++j)
594  {
595  v1t[j * v1t_dim1 + 1] = 0.;
596  v1t[j + v1t_dim1] = 0.;
597  }
598  i__1 = *q - 1;
599  i__2 = *q - 1;
600  dlacpy_("U", &i__1, &i__2, &x21[(x21_dim1 << 1) + 1], ldx21, &v1t[ (v1t_dim1 << 1) + 2], ldv1t);
601  i__1 = *q - 1;
602  i__2 = *q - 1;
603  i__3 = *q - 1;
604  dorglq_fla(&i__1, &i__2, &i__3, &v1t[(v1t_dim1 << 1) + 2], ldv1t, & work[itauq1], &work[iorglq], &lorglq, &childinfo);
605  }
606  /* Simultaneously diagonalize X11 and X21. */
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);
608  /* Permute rows and columns to place zero submatrices in */
609  /* preferred positions */
610  if (*q > 0 && wantu2)
611  {
612  i__1 = *q;
613  for (i__ = 1;
614  i__ <= i__1;
615  ++i__)
616  {
617  iwork[i__] = *m - *p - *q + i__;
618  }
619  i__1 = *m - *p;
620  for (i__ = *q + 1;
621  i__ <= i__1;
622  ++i__)
623  {
624  iwork[i__] = i__ - *q;
625  }
626  i__1 = *m - *p;
627  i__2 = *m - *p;
628  dlapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
629  }
630  }
631  else if (r__ == *p)
632  {
633  /* Case 2: R = P */
634  /* Simultaneously bidiagonalize X11 and X21 */
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);
636  /* Accumulate Householder reflectors */
637  if (wantu1 && *p > 0)
638  {
639  u1[u1_dim1 + 1] = 1.;
640  i__1 = *p;
641  for (j = 2;
642  j <= i__1;
643  ++j)
644  {
645  u1[j * u1_dim1 + 1] = 0.;
646  u1[j + u1_dim1] = 0.;
647  }
648  i__1 = *p - 1;
649  i__2 = *p - 1;
650  dlacpy_("L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &u1[( u1_dim1 << 1) + 2], ldu1);
651  i__1 = *p - 1;
652  i__2 = *p - 1;
653  i__3 = *p - 1;
654  dorgqr_fla(&i__1, &i__2, &i__3, &u1[(u1_dim1 << 1) + 2], ldu1, &work[ itaup1], &work[iorgqr], &lorgqr, &childinfo);
655  }
656  if (wantu2 && *m - *p > 0)
657  {
658  i__1 = *m - *p;
659  dlacpy_("L", &i__1, q, &x21[x21_offset], ldx21, &u2[u2_offset], ldu2);
660  i__1 = *m - *p;
661  i__2 = *m - *p;
662  dorgqr_fla(&i__1, &i__2, q, &u2[u2_offset], ldu2, &work[itaup2], & work[iorgqr], &lorgqr, &childinfo);
663  }
664  if (wantv1t && *q > 0)
665  {
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);
668  }
669  /* Simultaneously diagonalize X11 and X21. */
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);
671  /* Permute rows and columns to place identity submatrices in */
672  /* preferred positions */
673  if (*q > 0 && wantu2)
674  {
675  i__1 = *q;
676  for (i__ = 1;
677  i__ <= i__1;
678  ++i__)
679  {
680  iwork[i__] = *m - *p - *q + i__;
681  }
682  i__1 = *m - *p;
683  for (i__ = *q + 1;
684  i__ <= i__1;
685  ++i__)
686  {
687  iwork[i__] = i__ - *q;
688  }
689  i__1 = *m - *p;
690  i__2 = *m - *p;
691  dlapmt_(&c_false, &i__1, &i__2, &u2[u2_offset], ldu2, &iwork[1]);
692  }
693  }
694  else if (r__ == *m - *p)
695  {
696  /* Case 3: R = M-P */
697  /* Simultaneously bidiagonalize X11 and X21 */
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);
699  /* Accumulate Householder reflectors */
700  if (wantu1 && *p > 0)
701  {
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);
704  }
705  if (wantu2 && *m - *p > 0)
706  {
707  u2[u2_dim1 + 1] = 1.;
708  i__1 = *m - *p;
709  for (j = 2;
710  j <= i__1;
711  ++j)
712  {
713  u2[j * u2_dim1 + 1] = 0.;
714  u2[j + u2_dim1] = 0.;
715  }
716  i__1 = *m - *p - 1;
717  i__2 = *m - *p - 1;
718  dlacpy_("L", &i__1, &i__2, &x21[x21_dim1 + 2], ldx21, &u2[( u2_dim1 << 1) + 2], ldu2);
719  i__1 = *m - *p - 1;
720  i__2 = *m - *p - 1;
721  i__3 = *m - *p - 1;
722  dorgqr_fla(&i__1, &i__2, &i__3, &u2[(u2_dim1 << 1) + 2], ldu2, &work[ itaup2], &work[iorgqr], &lorgqr, &childinfo);
723  }
724  if (wantv1t && *q > 0)
725  {
726  i__1 = *m - *p;
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);
729  }
730  /* Simultaneously diagonalize X11 and X21. */
731  i__1 = *m - *q;
732  i__2 = *m - *p;
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);
734  /* Permute rows and columns to place identity submatrices in */
735  /* preferred positions */
736  if (*q > r__)
737  {
738  i__1 = r__;
739  for (i__ = 1;
740  i__ <= i__1;
741  ++i__)
742  {
743  iwork[i__] = *q - r__ + i__;
744  }
745  i__1 = *q;
746  for (i__ = r__ + 1;
747  i__ <= i__1;
748  ++i__)
749  {
750  iwork[i__] = i__ - r__;
751  }
752  if (wantu1)
753  {
754  dlapmt_(&c_false, p, q, &u1[u1_offset], ldu1, &iwork[1]);
755  }
756  if (wantv1t)
757  {
758  dlapmr_(&c_false, q, q, &v1t[v1t_offset], ldv1t, &iwork[1]);
759  }
760  }
761  }
762  else
763  {
764  /* Case 4: R = M-Q */
765  /* Simultaneously bidiagonalize X11 and X21 */
766  i__1 = lorbdb - *m;
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) ;
768  /* Accumulate Householder reflectors */
769  if (wantu1 && *p > 0)
770  {
771  dcopy_(p, &work[iorbdb], &c__1, &u1[u1_offset], &c__1);
772  i__1 = *p;
773  for (j = 2;
774  j <= i__1;
775  ++j)
776  {
777  u1[j * u1_dim1 + 1] = 0.;
778  }
779  i__1 = *p - 1;
780  i__2 = *m - *q - 1;
781  dlacpy_("L", &i__1, &i__2, &x11[x11_dim1 + 2], ldx11, &u1[( u1_dim1 << 1) + 2], ldu1);
782  i__1 = *m - *q;
783  dorgqr_fla(p, p, &i__1, &u1[u1_offset], ldu1, &work[itaup1], &work[ iorgqr], &lorgqr, &childinfo);
784  }
785  if (wantu2 && *m - *p > 0)
786  {
787  i__1 = *m - *p;
788  dcopy_(&i__1, &work[iorbdb + *p], &c__1, &u2[u2_offset], &c__1);
789  i__1 = *m - *p;
790  for (j = 2;
791  j <= i__1;
792  ++j)
793  {
794  u2[j * u2_dim1 + 1] = 0.;
795  }
796  i__1 = *m - *p - 1;
797  i__2 = *m - *q - 1;
798  dlacpy_("L", &i__1, &i__2, &x21[x21_dim1 + 2], ldx21, &u2[( u2_dim1 << 1) + 2], ldu2);
799  i__1 = *m - *p;
800  i__2 = *m - *p;
801  i__3 = *m - *q;
802  dorgqr_fla(&i__1, &i__2, &i__3, &u2[u2_offset], ldu2, &work[itaup2], &work[iorgqr], &lorgqr, &childinfo);
803  }
804  if (wantv1t && *q > 0)
805  {
806  i__1 = *m - *q;
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);
811  i__1 = -(*p) + *q;
812  i__2 = *q - *p;
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);
815  }
816  /* Simultaneously diagonalize X11 and X21. */
817  i__1 = *m - *p;
818  i__2 = *m - *q;
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);
820  /* Permute rows and columns to place identity submatrices in */
821  /* preferred positions */
822  if (*p > r__)
823  {
824  i__1 = r__;
825  for (i__ = 1;
826  i__ <= i__1;
827  ++i__)
828  {
829  iwork[i__] = *p - r__ + i__;
830  }
831  i__1 = *p;
832  for (i__ = r__ + 1;
833  i__ <= i__1;
834  ++i__)
835  {
836  iwork[i__] = i__ - r__;
837  }
838  if (wantu1)
839  {
840  dlapmt_(&c_false, p, p, &u1[u1_offset], ldu1, &iwork[1]);
841  }
842  if (wantv1t)
843  {
844  dlapmr_(&c_false, p, q, &v1t[v1t_offset], ldv1t, &iwork[1]);
845  }
846  }
847  }
848  return 0;
849  /* End of DORCSD2BY1 */
850 }
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
double doublereal
Definition: FLA_f2c.h:31
int logical
Definition: FLA_f2c.h:36
int integer
Definition: FLA_f2c.h:25
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