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 )
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}
integer
int integer
Definition FLA_f2c.h:25
logical
int logical
Definition FLA_f2c.h:36
real
float real
Definition FLA_f2c.h:30
sorglq_fla
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
sorgqr_fla
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

References sorglq_fla(), and sorgqr_fla().

Generated on Thu Sep 23 2021 17:46:17 for libflame by doxygen 1.12.0