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,
264 CHARACTER JOBU1, JOBU2, JOBV1T
265 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
267 INTEGER LRWORK, LRWORKMIN, LRWORKOPT
272 COMPLEX U1(ldu1,*), U2(ldu2,*), V1T(ldv1t,*), WORK(*),
273 $ x11(ldx11,*), x21(ldx21,*)
281 parameter ( one = (1.0e0,0.0e0), zero = (0.0e0,0.0e0) )
284 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
285 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
286 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
287 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
288 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
289 $ lworkmin, lworkopt, r
290 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
306 INTRINSIC int, max, min
313 wantu1 = lsame( jobu1,
'Y' )
314 wantu2 = lsame( jobu2,
'Y' )
315 wantv1t = lsame( jobv1t,
'Y' )
316 lquery = lwork .EQ. -1
320 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
322 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
324 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
326 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
328 ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) )
THEN
330 ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) )
THEN
332 ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) )
THEN
336 r = min( p, m-p, q, m-q )
371 IF( info .EQ. 0 )
THEN
373 ib11d = iphi + max( 1, r-1 )
374 ib11e = ib11d + max( 1, r )
375 ib12d = ib11e + max( 1, r - 1 )
376 ib12e = ib12d + max( 1, r )
377 ib21d = ib12e + max( 1, r - 1 )
378 ib21e = ib21d + max( 1, r )
379 ib22d = ib21e + max( 1, r - 1 )
380 ib22e = ib22d + max( 1, r )
381 ibbcsd = ib22e + max( 1, r - 1 )
383 itaup2 = itaup1 + max( 1, p )
384 itauq1 = itaup2 + max( 1, m-p )
385 iorbdb = itauq1 + max( 1, q )
386 iorgqr = itauq1 + max( 1, q )
387 iorglq = itauq1 + max( 1, q )
393 CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
394 $ dum, cdum, cdum, cdum, work, -1,
396 lorbdb = int( work(1) )
397 IF( wantu1 .AND. p .GT. 0 )
THEN
398 CALL cungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
400 lorgqrmin = max( lorgqrmin, p )
401 lorgqropt = max( lorgqropt, int( work(1) ) )
403 IF( wantu2 .AND. m-p .GT. 0 )
THEN
404 CALL cungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
406 lorgqrmin = max( lorgqrmin, m-p )
407 lorgqropt = max( lorgqropt, int( work(1) ) )
409 IF( wantv1t .AND. q .GT. 0 )
THEN
410 CALL cunglq( q-1, q-1, q-1, v1t, ldv1t,
411 $ cdum, work(1), -1, childinfo )
412 lorglqmin = max( lorglqmin, q-1 )
413 lorglqopt = max( lorglqopt, int( work(1) ) )
415 CALL cbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
416 $ dum(1), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
417 $ 1, dum, dum, dum, dum, dum, dum, dum, dum,
418 $ rwork(1), -1, childinfo )
419 lbbcsd = int( rwork(1) )
420 ELSE IF( r .EQ. p )
THEN
421 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
422 $ cdum, cdum, cdum, work(1), -1, childinfo )
423 lorbdb = int( work(1) )
424 IF( wantu1 .AND. p .GT. 0 )
THEN
425 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum, work(1),
427 lorgqrmin = max( lorgqrmin, p-1 )
428 lorgqropt = max( lorgqropt, int( work(1) ) )
430 IF( wantu2 .AND. m-p .GT. 0 )
THEN
431 CALL cungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
433 lorgqrmin = max( lorgqrmin, m-p )
434 lorgqropt = max( lorgqropt, int( work(1) ) )
436 IF( wantv1t .AND. q .GT. 0 )
THEN
437 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
439 lorglqmin = max( lorglqmin, q )
440 lorglqopt = max( lorglqopt, int( work(1) ) )
442 CALL cbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
443 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
444 $ dum, dum, dum, dum, dum, dum, dum, dum,
445 $ rwork(1), -1, childinfo )
446 lbbcsd = int( rwork(1) )
447 ELSE IF( r .EQ. m-p )
THEN
448 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
449 $ cdum, cdum, cdum, work(1), -1, childinfo )
450 lorbdb = int( work(1) )
451 IF( wantu1 .AND. p .GT. 0 )
THEN
452 CALL cungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
454 lorgqrmin = max( lorgqrmin, p )
455 lorgqropt = max( lorgqropt, int( work(1) ) )
457 IF( wantu2 .AND. m-p .GT. 0 )
THEN
458 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
459 $ work(1), -1, childinfo )
460 lorgqrmin = max( lorgqrmin, m-p-1 )
461 lorgqropt = max( lorgqropt, int( work(1) ) )
463 IF( wantv1t .AND. q .GT. 0 )
THEN
464 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
466 lorglqmin = max( lorglqmin, q )
467 lorglqopt = max( lorglqopt, int( work(1) ) )
469 CALL cbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
470 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
471 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
472 $ rwork(1), -1, childinfo )
473 lbbcsd = int( rwork(1) )
475 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
476 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
478 lorbdb = m + int( work(1) )
479 IF( wantu1 .AND. p .GT. 0 )
THEN
480 CALL cungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
482 lorgqrmin = max( lorgqrmin, p )
483 lorgqropt = max( lorgqropt, int( work(1) ) )
485 IF( wantu2 .AND. m-p .GT. 0 )
THEN
486 CALL cungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1), -1,
488 lorgqrmin = max( lorgqrmin, m-p )
489 lorgqropt = max( lorgqropt, int( work(1) ) )
491 IF( wantv1t .AND. q .GT. 0 )
THEN
492 CALL cunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
494 lorglqmin = max( lorglqmin, q )
495 lorglqopt = max( lorglqopt, int( work(1) ) )
497 CALL cbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
498 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
499 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
500 $ rwork(1), -1, childinfo )
501 lbbcsd = int( rwork(1) )
503 lrworkmin = ibbcsd+lbbcsd-1
504 lrworkopt = lrworkmin
506 lworkmin = max( iorbdb+lorbdb-1,
507 $ iorgqr+lorgqrmin-1,
508 $ iorglq+lorglqmin-1 )
509 lworkopt = max( iorbdb+lorbdb-1,
510 $ iorgqr+lorgqropt-1,
511 $ iorglq+lorglqopt-1 )
513 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
517 IF( info .NE. 0 )
THEN
518 CALL xerbla(
'CUNCSD2BY1', -info )
520 ELSE IF( lquery )
THEN
523 lorgqr = lwork-iorgqr+1
524 lorglq = lwork-iorglq+1
535 CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
536 $ rwork(iphi), work(itaup1), work(itaup2),
537 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
541 IF( wantu1 .AND. p .GT. 0 )
THEN
542 CALL clacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
543 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
544 $ lorgqr, childinfo )
546 IF( wantu2 .AND. m-p .GT. 0 )
THEN
547 CALL clacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
548 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
549 $ work(iorgqr), lorgqr, childinfo )
551 IF( wantv1t .AND. q .GT. 0 )
THEN
557 CALL clacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
559 CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
560 $ work(iorglq), lorglq, childinfo )
565 CALL cbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
566 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
567 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
568 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
569 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
575 IF( q .GT. 0 .AND. wantu2 )
THEN
577 iwork(i) = m - p - q + i
582 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
584 ELSE IF( r .EQ. p )
THEN
590 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
591 $ rwork(iphi), work(itaup1), work(itaup2),
592 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
596 IF( wantu1 .AND. p .GT. 0 )
THEN
602 CALL clacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
603 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
604 $ work(iorgqr), lorgqr, childinfo )
606 IF( wantu2 .AND. m-p .GT. 0 )
THEN
607 CALL clacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
608 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
609 $ work(iorgqr), lorgqr, childinfo )
611 IF( wantv1t .AND. q .GT. 0 )
THEN
612 CALL clacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
613 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
614 $ work(iorglq), lorglq, childinfo )
619 CALL cbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
620 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
621 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
622 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
623 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
629 IF( q .GT. 0 .AND. wantu2 )
THEN
631 iwork(i) = m - p - q + i
636 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
638 ELSE IF( r .EQ. m-p )
THEN
644 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
645 $ rwork(iphi), work(itaup1), work(itaup2),
646 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
650 IF( wantu1 .AND. p .GT. 0 )
THEN
651 CALL clacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
652 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
653 $ lorgqr, childinfo )
655 IF( wantu2 .AND. m-p .GT. 0 )
THEN
661 CALL clacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
663 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
664 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
666 IF( wantv1t .AND. q .GT. 0 )
THEN
667 CALL clacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
668 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
669 $ work(iorglq), lorglq, childinfo )
674 CALL cbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
675 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
676 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
677 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
678 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
679 $ rwork(ibbcsd), lbbcsd, childinfo )
692 CALL clapmt( .false., p, q, u1, ldu1, iwork )
695 CALL clapmr( .false., q, q, v1t, ldv1t, iwork )
704 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
705 $ rwork(iphi), work(itaup1), work(itaup2),
706 $ work(itauq1), work(iorbdb), work(iorbdb+m),
707 $ lorbdb-m, childinfo )
711 IF( wantu1 .AND. p .GT. 0 )
THEN
712 CALL ccopy( p, work(iorbdb), 1, u1, 1 )
716 CALL clacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
718 CALL cungqr( p, p, m-q, u1, ldu1, work(itaup1),
719 $ work(iorgqr), lorgqr, childinfo )
721 IF( wantu2 .AND. m-p .GT. 0 )
THEN
722 CALL ccopy( m-p, work(iorbdb+p), 1, u2, 1 )
726 CALL clacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
728 CALL cungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
729 $ work(iorgqr), lorgqr, childinfo )
731 IF( wantv1t .AND. q .GT. 0 )
THEN
732 CALL clacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
733 CALL clacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
734 $ v1t(m-q+1,m-q+1), ldv1t )
735 CALL clacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
736 $ v1t(p+1,p+1), ldv1t )
737 CALL cunglq( q, q, q, v1t, ldv1t, work(itauq1),
738 $ work(iorglq), lorglq, childinfo )
743 CALL cbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
744 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
745 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
746 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
747 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
748 $ rwork(ibbcsd), lbbcsd, childinfo )
761 CALL clapmt( .false., p, p, u1, ldu1, iwork )
764 CALL clapmr( .false., p, q, v1t, ldv1t, iwork )
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
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 cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
subroutine clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine cunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB1
subroutine cunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
CUNBDB4
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
CUNBDB3
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 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 cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ