253 SUBROUTINE cuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
254 $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
255 $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK,
263 CHARACTER JOBU1, JOBU2, JOBV1T
264 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
266 INTEGER LRWORK, LRWORKMIN, LRWORKOPT
271 COMPLEX U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
272 $ x11(ldx11,*), x21(ldx21,*)
280 PARAMETER ( ONE = (1.0e0,0.0e0), zero = (0.0e0,0.0e0) )
283 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
284 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
285 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
286 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
287 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
288 $ lworkmin, lworkopt, r
289 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
305 INTRINSIC int, max, min
312 wantu1 = lsame( jobu1,
'Y' )
313 wantu2 = lsame( jobu2,
'Y' )
314 wantv1t = lsame( jobv1t,
'Y' )
315 lquery = ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 )
319 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
321 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
323 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
325 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
327 ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) )
THEN
329 ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) )
THEN
331 ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) )
THEN
335 r = min( p, m-p, q, m-q )
370 IF( info .EQ. 0 )
THEN
372 ib11d = iphi + max( 1, r-1 )
373 ib11e = ib11d + max( 1, r )
374 ib12d = ib11e + max( 1, r - 1 )
375 ib12e = ib12d + max( 1, r )
376 ib21d = ib12e + max( 1, r - 1 )
377 ib21e = ib21d + max( 1, r )
378 ib22d = ib21e + max( 1, r - 1 )
379 ib22e = ib22d + max( 1, r )
380 ibbcsd = ib22e + max( 1, r - 1 )
382 itaup2 = itaup1 + max( 1, p )
383 itauq1 = itaup2 + max( 1, m-p )
384 iorbdb = itauq1 + max( 1, q )
385 iorgqr = itauq1 + max( 1, q )
386 iorglq = itauq1 + max( 1, q )
392 CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
393 $ dum, cdum, cdum, cdum, work, -1,
395 lorbdb = int( work(1) )
396 IF( wantu1 .AND. p .GT. 0 )
THEN
397 CALL cungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
399 lorgqrmin = max( lorgqrmin, p )
400 lorgqropt = max( lorgqropt, int( work(1) ) )
402 IF( wantu2 .AND. m-p .GT. 0 )
THEN
403 CALL cungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
405 lorgqrmin = max( lorgqrmin, m-p )
406 lorgqropt = max( lorgqropt, int( work(1) ) )
408 IF( wantv1t .AND. q .GT. 0 )
THEN
409 CALL cunglq( q-1, q-1, q-1, v1t, ldv1t,
410 $ cdum, work(1), -1, childinfo )
411 lorglqmin = max( lorglqmin, q-1 )
412 lorglqopt = max( lorglqopt, int( work(1) ) )
414 CALL cbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
415 $ dum(1), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
416 $ 1, dum, dum, dum, dum, dum, dum, dum, dum,
417 $ rwork(1), -1, childinfo )
418 lbbcsd = int( rwork(1) )
419 ELSE IF( r .EQ. p )
THEN
420 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
421 $ cdum, cdum, cdum, work(1), -1, childinfo )
422 lorbdb = int( work(1) )
423 IF( wantu1 .AND. p .GT. 0 )
THEN
424 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum, work(1),
426 lorgqrmin = max( lorgqrmin, p-1 )
427 lorgqropt = max( lorgqropt, int( work(1) ) )
429 IF( wantu2 .AND. m-p .GT. 0 )
THEN
430 CALL cungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
432 lorgqrmin = max( lorgqrmin, m-p )
433 lorgqropt = max( lorgqropt, int( work(1) ) )
435 IF( wantv1t .AND. q .GT. 0 )
THEN
436 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
438 lorglqmin = max( lorglqmin, q )
439 lorglqopt = max( lorglqopt, int( work(1) ) )
441 CALL cbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
442 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
443 $ dum, dum, dum, dum, dum, dum, dum, dum,
444 $ rwork(1), -1, childinfo )
445 lbbcsd = int( rwork(1) )
446 ELSE IF( r .EQ. m-p )
THEN
447 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
448 $ cdum, cdum, cdum, work(1), -1, childinfo )
449 lorbdb = int( work(1) )
450 IF( wantu1 .AND. p .GT. 0 )
THEN
451 CALL cungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
453 lorgqrmin = max( lorgqrmin, p )
454 lorgqropt = max( lorgqropt, int( work(1) ) )
456 IF( wantu2 .AND. m-p .GT. 0 )
THEN
457 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
458 $ work(1), -1, childinfo )
459 lorgqrmin = max( lorgqrmin, m-p-1 )
460 lorgqropt = max( lorgqropt, int( work(1) ) )
462 IF( wantv1t .AND. q .GT. 0 )
THEN
463 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
465 lorglqmin = max( lorglqmin, q )
466 lorglqopt = max( lorglqopt, int( work(1) ) )
468 CALL cbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
469 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
470 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
471 $ rwork(1), -1, childinfo )
472 lbbcsd = int( rwork(1) )
474 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
475 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
477 lorbdb = m + int( work(1) )
478 IF( wantu1 .AND. p .GT. 0 )
THEN
479 CALL cungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
481 lorgqrmin = max( lorgqrmin, p )
482 lorgqropt = max( lorgqropt, int( work(1) ) )
484 IF( wantu2 .AND. m-p .GT. 0 )
THEN
485 CALL cungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1), -1,
487 lorgqrmin = max( lorgqrmin, m-p )
488 lorgqropt = max( lorgqropt, int( work(1) ) )
490 IF( wantv1t .AND. q .GT. 0 )
THEN
491 CALL cunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
493 lorglqmin = max( lorglqmin, q )
494 lorglqopt = max( lorglqopt, int( work(1) ) )
496 CALL cbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
497 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
498 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
499 $ rwork(1), -1, childinfo )
500 lbbcsd = int( rwork(1) )
502 lrworkmin = ibbcsd+lbbcsd-1
503 lrworkopt = lrworkmin
505 lworkmin = max( iorbdb+lorbdb-1,
506 $ iorgqr+lorgqrmin-1,
507 $ iorglq+lorglqmin-1 )
508 lworkopt = max( iorbdb+lorbdb-1,
509 $ iorgqr+lorgqropt-1,
510 $ iorglq+lorglqopt-1 )
512 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
515 IF( lrwork .LT. lrworkmin .AND. .NOT.lquery )
THEN
519 IF( info .NE. 0 )
THEN
520 CALL xerbla(
'CUNCSD2BY1', -info )
522 ELSE IF( lquery )
THEN
525 lorgqr = lwork-iorgqr+1
526 lorglq = lwork-iorglq+1
537 CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
538 $ rwork(iphi), work(itaup1), work(itaup2),
539 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
543 IF( wantu1 .AND. p .GT. 0 )
THEN
544 CALL clacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
545 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
546 $ lorgqr, childinfo )
548 IF( wantu2 .AND. m-p .GT. 0 )
THEN
549 CALL clacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
550 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
551 $ work(iorgqr), lorgqr, childinfo )
553 IF( wantv1t .AND. q .GT. 0 )
THEN
559 CALL clacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
561 CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
562 $ work(iorglq), lorglq, childinfo )
567 CALL cbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
568 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
569 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
570 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
571 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
572 $ lrwork-ibbcsd+1, childinfo )
577 IF( q .GT. 0 .AND. wantu2 )
THEN
579 iwork(i) = m - p - q + i
584 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
586 ELSE IF( r .EQ. p )
THEN
592 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
593 $ rwork(iphi), work(itaup1), work(itaup2),
594 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
598 IF( wantu1 .AND. p .GT. 0 )
THEN
604 CALL clacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
605 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
606 $ work(iorgqr), lorgqr, childinfo )
608 IF( wantu2 .AND. m-p .GT. 0 )
THEN
609 CALL clacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
610 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
611 $ work(iorgqr), lorgqr, childinfo )
613 IF( wantv1t .AND. q .GT. 0 )
THEN
614 CALL clacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
615 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
616 $ work(iorglq), lorglq, childinfo )
621 CALL cbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
622 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
623 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
624 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
625 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
631 IF( q .GT. 0 .AND. wantu2 )
THEN
633 iwork(i) = m - p - q + i
638 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
640 ELSE IF( r .EQ. m-p )
THEN
646 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
647 $ rwork(iphi), work(itaup1), work(itaup2),
648 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
652 IF( wantu1 .AND. p .GT. 0 )
THEN
653 CALL clacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
654 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
655 $ lorgqr, childinfo )
657 IF( wantu2 .AND. m-p .GT. 0 )
THEN
663 CALL clacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
665 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
666 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
668 IF( wantv1t .AND. q .GT. 0 )
THEN
669 CALL clacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
670 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
671 $ work(iorglq), lorglq, childinfo )
676 CALL cbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
677 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
678 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
679 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
680 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
681 $ rwork(ibbcsd), lbbcsd, childinfo )
694 CALL clapmt( .false., p, q, u1, ldu1, iwork )
697 CALL clapmr( .false., q, q, v1t, ldv1t, iwork )
706 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
707 $ rwork(iphi), work(itaup1), work(itaup2),
708 $ work(itauq1), work(iorbdb), work(iorbdb+m),
709 $ lorbdb-m, childinfo )
714 IF( wantu2 .AND. m-p .GT. 0 )
THEN
715 CALL ccopy( m-p, work(iorbdb+p), 1, u2, 1 )
717 IF( wantu1 .AND. p .GT. 0 )
THEN
718 CALL ccopy( p, work(iorbdb), 1, u1, 1 )
722 CALL clacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
724 CALL cungqr( p, p, m-q, u1, ldu1, work(itaup1),
725 $ work(iorgqr), lorgqr, childinfo )
727 IF( wantu2 .AND. m-p .GT. 0 )
THEN
731 CALL clacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
733 CALL cungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
734 $ work(iorgqr), lorgqr, childinfo )
736 IF( wantv1t .AND. q .GT. 0 )
THEN
737 CALL clacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
738 CALL clacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
739 $ v1t(m-q+1,m-q+1), ldv1t )
740 CALL clacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
741 $ v1t(p+1,p+1), ldv1t )
742 CALL cunglq( q, q, q, v1t, ldv1t, work(itauq1),
743 $ work(iorglq), lorglq, childinfo )
748 CALL cbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
749 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
750 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
751 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
752 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
753 $ rwork(ibbcsd), lbbcsd, childinfo )
766 CALL clapmt( .false., p, p, u1, ldu1, iwork )
769 CALL clapmr( .false., p, q, v1t, ldv1t, iwork )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cunbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB2
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
subroutine cunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
CUNBDB4
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 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 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
subroutine cunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB1