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