260 SUBROUTINE cuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
261 $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
262 $ ldv1t, work, lwork, rwork, lrwork, iwork,
271 CHARACTER jobu1, jobu2, jobv1t
272 INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
274 INTEGER lrwork, lrworkmin, lrworkopt
279 COMPLEX u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
280 $ x11(ldx11,*), x21(ldx21,*)
288 parameter( one = (1.0e0,0.0e0), zero = (0.0e0,0.0e0) )
291 INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
292 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
293 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
294 $
j, lbbcsd, lorbdb, lorglq, lorglqmin,
295 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
296 $ lworkmin, lworkopt, r
297 LOGICAL lquery, wantu1, wantu2, wantv1t
309 INTRINSIC int, max, min
316 wantu1 =
lsame( jobu1,
'Y' )
317 wantu2 =
lsame( jobu2,
'Y' )
318 wantv1t =
lsame( jobv1t,
'Y' )
319 lquery = lwork .EQ. -1
323 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
325 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
327 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
329 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
331 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
333 ELSE IF( wantu2 .AND. ldu2 .LT. m - p )
THEN
335 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
339 r = min( p, m-p, q, m-q )
374 IF( info .EQ. 0 )
THEN
376 ib11d = iphi + max( 1, r-1 )
377 ib11e = ib11d + max( 1, r )
378 ib12d = ib11e + max( 1, r - 1 )
379 ib12e = ib12d + max( 1, r )
380 ib21d = ib12e + max( 1, r - 1 )
381 ib21e = ib21d + max( 1, r )
382 ib22d = ib21e + max( 1, r - 1 )
383 ib22e = ib22d + max( 1, r )
384 ibbcsd = ib22e + max( 1, r - 1 )
386 itaup2 = itaup1 + max( 1, p )
387 itauq1 = itaup2 + max( 1, m-p )
388 iorbdb = itauq1 + max( 1, q )
389 iorgqr = itauq1 + max( 1, q )
390 iorglq = itauq1 + max( 1, q )
392 CALL
cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
393 $ 0, 0, work, -1, childinfo )
394 lorbdb = int( work(1) )
395 IF( p .GE. m-p )
THEN
396 CALL
cungqr( p, p, q, u1, ldu1, 0, work(1), -1,
398 lorgqrmin = max( 1, p )
399 lorgqropt = int( work(1) )
401 CALL
cungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
403 lorgqrmin = max( 1, m-p )
404 lorgqropt = int( work(1) )
406 CALL
cunglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
407 $ 0, work(1), -1, childinfo )
408 lorglqmin = max( 1, q-1 )
409 lorglqopt = int( work(1) )
410 CALL
cbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
411 $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
412 $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
413 lbbcsd = int( rwork(1) )
414 ELSE IF( r .EQ. p )
THEN
415 CALL
cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
416 $ 0, 0, work(1), -1, childinfo )
417 lorbdb = int( work(1) )
418 IF( p-1 .GE. m-p )
THEN
419 CALL
cungqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
421 lorgqrmin = max( 1, p-1 )
422 lorgqropt = int( work(1) )
424 CALL
cungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
426 lorgqrmin = max( 1, m-p )
427 lorgqropt = int( work(1) )
429 CALL
cunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
431 lorglqmin = max( 1, q )
432 lorglqopt = int( work(1) )
433 CALL
cbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
434 $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
435 $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
436 lbbcsd = int( rwork(1) )
437 ELSE IF( r .EQ. m-p )
THEN
438 CALL
cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
439 $ 0, 0, work(1), -1, childinfo )
440 lorbdb = int( work(1) )
441 IF( p .GE. m-p-1 )
THEN
442 CALL
cungqr( p, p, q, u1, ldu1, 0, work(1), -1,
444 lorgqrmin = max( 1, p )
445 lorgqropt = int( work(1) )
447 CALL
cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
448 $ work(1), -1, childinfo )
449 lorgqrmin = max( 1, m-p-1 )
450 lorgqropt = int( work(1) )
452 CALL
cunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
454 lorglqmin = max( 1, q )
455 lorglqopt = int( work(1) )
456 CALL
cbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
457 $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
458 $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
460 lbbcsd = int( rwork(1) )
462 CALL
cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
463 $ 0, 0, 0, work(1), -1, childinfo )
464 lorbdb = m + int( work(1) )
465 IF( p .GE. m-p )
THEN
466 CALL
cungqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
468 lorgqrmin = max( 1, p )
469 lorgqropt = int( work(1) )
471 CALL
cungqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
473 lorgqrmin = max( 1, m-p )
474 lorgqropt = int( work(1) )
476 CALL
cunglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
478 lorglqmin = max( 1, q )
479 lorglqopt = int( work(1) )
480 CALL
cbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
481 $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
482 $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
484 lbbcsd = int( rwork(1) )
486 lrworkmin = ibbcsd+lbbcsd-1
487 lrworkopt = lrworkmin
489 lworkmin = max( iorbdb+lorbdb-1,
490 $ iorgqr+lorgqrmin-1,
491 $ iorglq+lorglqmin-1 )
492 lworkopt = max( iorbdb+lorbdb-1,
493 $ iorgqr+lorgqropt-1,
494 $ iorglq+lorglqopt-1 )
496 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
500 IF( info .NE. 0 )
THEN
501 CALL
xerbla(
'CUNCSD2BY1', -info )
503 ELSE IF( lquery )
THEN
506 lorgqr = lwork-iorgqr+1
507 lorglq = lwork-iorglq+1
518 CALL
cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
519 $ rwork(iphi), work(itaup1), work(itaup2),
520 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
524 IF( wantu1 .AND. p .GT. 0 )
THEN
525 CALL
clacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
526 CALL
cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
527 $ lorgqr, childinfo )
529 IF( wantu2 .AND. m-p .GT. 0 )
THEN
530 CALL
clacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
531 CALL
cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
532 $ work(iorgqr), lorgqr, childinfo )
534 IF( wantv1t .AND. q .GT. 0 )
THEN
540 CALL
clacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
542 CALL
cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
543 $ work(iorglq), lorglq, childinfo )
548 CALL
cbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
549 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
550 $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
551 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
552 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
558 IF( q .GT. 0 .AND. wantu2 )
THEN
560 iwork(i) = m - p - q + i
565 CALL
clapmt( .false., m-p, m-p, u2, ldu2, iwork )
567 ELSE IF( r .EQ. p )
THEN
573 CALL
cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
574 $ rwork(iphi), work(itaup1), work(itaup2),
575 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
579 IF( wantu1 .AND. p .GT. 0 )
THEN
585 CALL
clacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
586 CALL
cungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
587 $ work(iorgqr), lorgqr, childinfo )
589 IF( wantu2 .AND. m-p .GT. 0 )
THEN
590 CALL
clacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
591 CALL
cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
592 $ work(iorgqr), lorgqr, childinfo )
594 IF( wantv1t .AND. q .GT. 0 )
THEN
595 CALL
clacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
596 CALL
cunglq( q, q, r, v1t, ldv1t, work(itauq1),
597 $ work(iorglq), lorglq, childinfo )
602 CALL
cbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
603 $ rwork(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
604 $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
605 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
606 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
612 IF( q .GT. 0 .AND. wantu2 )
THEN
614 iwork(i) = m - p - q + i
619 CALL
clapmt( .false., m-p, m-p, u2, ldu2, iwork )
621 ELSE IF( r .EQ. m-p )
THEN
627 CALL
cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
628 $ rwork(iphi), work(itaup1), work(itaup2),
629 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
633 IF( wantu1 .AND. p .GT. 0 )
THEN
634 CALL
clacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
635 CALL
cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
636 $ lorgqr, childinfo )
638 IF( wantu2 .AND. m-p .GT. 0 )
THEN
644 CALL
clacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
646 CALL
cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
647 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
649 IF( wantv1t .AND. q .GT. 0 )
THEN
650 CALL
clacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
651 CALL
cunglq( q, q, r, v1t, ldv1t, work(itauq1),
652 $ work(iorglq), lorglq, childinfo )
657 CALL
cbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
658 $ theta, rwork(iphi), 0, 1, v1t, ldv1t, u2, ldu2,
659 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
660 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
661 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
662 $ rwork(ibbcsd), lbbcsd, childinfo )
675 CALL
clapmt( .false., p, q, u1, ldu1, iwork )
678 CALL
clapmr( .false., q, q, v1t, ldv1t, iwork )
687 CALL
cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
688 $ rwork(iphi), work(itaup1), work(itaup2),
689 $ work(itauq1), work(iorbdb), work(iorbdb+m),
690 $ lorbdb-m, childinfo )
694 IF( wantu1 .AND. p .GT. 0 )
THEN
695 CALL
ccopy( p, work(iorbdb), 1, u1, 1 )
699 CALL
clacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
701 CALL
cungqr( p, p, m-q, u1, ldu1, work(itaup1),
702 $ work(iorgqr), lorgqr, childinfo )
704 IF( wantu2 .AND. m-p .GT. 0 )
THEN
705 CALL
ccopy( m-p, work(iorbdb+p), 1, u2, 1 )
709 CALL
clacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
711 CALL
cungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
712 $ work(iorgqr), lorgqr, childinfo )
714 IF( wantv1t .AND. q .GT. 0 )
THEN
715 CALL
clacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
716 CALL
clacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
717 $ v1t(m-q+1,m-q+1), ldv1t )
718 CALL
clacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
719 $ v1t(p+1,p+1), ldv1t )
720 CALL
cunglq( q, q, q, v1t, ldv1t, work(itauq1),
721 $ work(iorglq), lorglq, childinfo )
726 CALL
cbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
727 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
728 $ ldv1t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
729 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
730 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
744 CALL
clapmt( .false., p, p, u1, ldu1, iwork )
747 CALL
clapmr( .false., p, q, v1t, ldv1t, iwork )