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 
)
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 */
251 dorglq_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *),
252 dorgqr_fla(integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *);
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}
doublereal
double doublereal
Definition FLA_f2c.h:31
integer
int integer
Definition FLA_f2c.h:25
logical
int logical
Definition FLA_f2c.h:36
i
int i
Definition bl1_axmyv2.c:145
dorglq_fla
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
dorgqr_fla
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

References dorglq_fla(), dorgqr_fla(), and i.

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