libflame  revision_anchor
Functions
sorcsd2by1.c File Reference

(r)

Functions

int sorcsd2by1_ (char *jobu1, char *jobu2, char *jobv1t, integer *m, integer *p, integer *q, real *x11, integer *ldx11, real *x21, integer *ldx21, real *theta, real *u1, integer *ldu1, real *u2, integer *ldu2, real *v1t, integer *ldv1t, real *work, integer *lwork, integer *iwork, integer *info)
 

Function Documentation

◆ sorcsd2by1_()

int sorcsd2by1_ ( char *  jobu1,
char *  jobu2,
char *  jobv1t,
integer m,
integer p,
integer q,
real x11,
integer ldx11,
real x21,
integer ldx21,
real theta,
real u1,
integer ldu1,
real u2,
integer ldu2,
real v1t,
integer ldv1t,
real work,
integer lwork,
integer iwork,
integer info 
)

References sorglq_fla(), and sorgqr_fla().

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