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 )
subroutine cuncsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
CUNCSD2BY1
subroutine clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB1
subroutine cbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, RWORK, LRWORK, INFO)
CBBCSD
subroutine cunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
CUNBDB4
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cunbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB2
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
subroutine cunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB3