252 SUBROUTINE zuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
253 $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
254 $ 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
269 DOUBLE PRECISION RWORK(*)
270 DOUBLE PRECISION THETA(*)
271 COMPLEX*16 U1(ldu1,*), U2(ldu2,*), V1T(ldv1t,*), WORK(*),
272 $ x11(ldx11,*), x21(ldx21,*)
280 parameter ( one = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
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
292 DOUBLE PRECISION DUM( 1 )
293 COMPLEX*16 CDUM( 1, 1 )
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
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 zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
393 $ cdum, cdum, cdum, work, -1, childinfo )
394 lorbdb = int( work(1) )
395 IF( wantu1 .AND. p .GT. 0 )
THEN
396 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
398 lorgqrmin = max( lorgqrmin, p )
399 lorgqropt = max( lorgqropt, int( work(1) ) )
401 IF( wantu2 .AND. m-p .GT. 0 )
THEN
402 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
404 lorgqrmin = max( lorgqrmin, m-p )
405 lorgqropt = max( lorgqropt, int( work(1) ) )
407 IF( wantv1t .AND. q .GT. 0 )
THEN
408 CALL zunglq( q-1, q-1, q-1, v1t, ldv1t,
409 $ cdum, work(1), -1, childinfo )
410 lorglqmin = max( lorglqmin, q-1 )
411 lorglqopt = max( lorglqopt, int( work(1) ) )
413 CALL zbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
414 $ dum, u1, ldu1, u2, ldu2, v1t, ldv1t, cdum, 1,
415 $ dum, dum, dum, dum, dum, dum, dum, dum,
416 $ rwork(1), -1, childinfo )
417 lbbcsd = int( rwork(1) )
418 ELSE IF( r .EQ. p )
THEN
419 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
420 $ cdum, cdum, cdum, work(1), -1, childinfo )
421 lorbdb = int( work(1) )
422 IF( wantu1 .AND. p .GT. 0 )
THEN
423 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum, work(1),
425 lorgqrmin = max( lorgqrmin, p-1 )
426 lorgqropt = max( lorgqropt, int( work(1) ) )
428 IF( wantu2 .AND. m-p .GT. 0 )
THEN
429 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
431 lorgqrmin = max( lorgqrmin, m-p )
432 lorgqropt = max( lorgqropt, int( work(1) ) )
434 IF( wantv1t .AND. q .GT. 0 )
THEN
435 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
437 lorglqmin = max( lorglqmin, q )
438 lorglqopt = max( lorglqopt, int( work(1) ) )
440 CALL zbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
441 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
442 $ dum, dum, dum, dum, dum, dum, dum, dum,
443 $ rwork(1), -1, childinfo )
444 lbbcsd = int( rwork(1) )
445 ELSE IF( r .EQ. m-p )
THEN
446 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
447 $ cdum, cdum, cdum, work(1), -1, childinfo )
448 lorbdb = int( work(1) )
449 IF( wantu1 .AND. p .GT. 0 )
THEN
450 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
452 lorgqrmin = max( lorgqrmin, p )
453 lorgqropt = max( lorgqropt, int( work(1) ) )
455 IF( wantu2 .AND. m-p .GT. 0 )
THEN
456 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
457 $ work(1), -1, childinfo )
458 lorgqrmin = max( lorgqrmin, m-p-1 )
459 lorgqropt = max( lorgqropt, int( work(1) ) )
461 IF( wantv1t .AND. q .GT. 0 )
THEN
462 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
464 lorglqmin = max( lorglqmin, q )
465 lorglqopt = max( lorglqopt, int( work(1) ) )
467 CALL zbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
468 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
469 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
470 $ rwork(1), -1, childinfo )
471 lbbcsd = int( rwork(1) )
473 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
474 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
476 lorbdb = m + int( work(1) )
477 IF( wantu1 .AND. p .GT. 0 )
THEN
478 CALL zungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
480 lorgqrmin = max( lorgqrmin, p )
481 lorgqropt = max( lorgqropt, int( work(1) ) )
483 IF( wantu2 .AND. m-p .GT. 0 )
THEN
484 CALL zungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1), -1,
486 lorgqrmin = max( lorgqrmin, m-p )
487 lorgqropt = max( lorgqropt, int( work(1) ) )
489 IF( wantv1t .AND. q .GT. 0 )
THEN
490 CALL zunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
492 lorglqmin = max( lorglqmin, q )
493 lorglqopt = max( lorglqopt, int( work(1) ) )
495 CALL zbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
496 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
497 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
498 $ rwork(1), -1, childinfo )
499 lbbcsd = int( rwork(1) )
501 lrworkmin = ibbcsd+lbbcsd-1
502 lrworkopt = lrworkmin
504 lworkmin = max( iorbdb+lorbdb-1,
505 $ iorgqr+lorgqrmin-1,
506 $ iorglq+lorglqmin-1 )
507 lworkopt = max( iorbdb+lorbdb-1,
508 $ iorgqr+lorgqropt-1,
509 $ iorglq+lorglqopt-1 )
511 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
515 IF( info .NE. 0 )
THEN
516 CALL xerbla(
'ZUNCSD2BY1', -info )
518 ELSE IF( lquery )
THEN
521 lorgqr = lwork-iorgqr+1
522 lorglq = lwork-iorglq+1
533 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
534 $ rwork(iphi), work(itaup1), work(itaup2),
535 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
539 IF( wantu1 .AND. p .GT. 0 )
THEN
540 CALL zlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
541 CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
542 $ lorgqr, childinfo )
544 IF( wantu2 .AND. m-p .GT. 0 )
THEN
545 CALL zlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
546 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
547 $ work(iorgqr), lorgqr, childinfo )
549 IF( wantv1t .AND. q .GT. 0 )
THEN
555 CALL zlacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
557 CALL zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
558 $ work(iorglq), lorglq, childinfo )
563 CALL zbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
564 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
565 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
566 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
567 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
573 IF( q .GT. 0 .AND. wantu2 )
THEN
575 iwork(i) = m - p - q + i
580 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
582 ELSE IF( r .EQ. p )
THEN
588 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
589 $ rwork(iphi), work(itaup1), work(itaup2),
590 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
594 IF( wantu1 .AND. p .GT. 0 )
THEN
600 CALL zlacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
601 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
602 $ work(iorgqr), lorgqr, childinfo )
604 IF( wantu2 .AND. m-p .GT. 0 )
THEN
605 CALL zlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
606 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
607 $ work(iorgqr), lorgqr, childinfo )
609 IF( wantv1t .AND. q .GT. 0 )
THEN
610 CALL zlacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
611 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
612 $ work(iorglq), lorglq, childinfo )
617 CALL zbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
618 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
619 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
620 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
621 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
627 IF( q .GT. 0 .AND. wantu2 )
THEN
629 iwork(i) = m - p - q + i
634 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
636 ELSE IF( r .EQ. m-p )
THEN
642 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
643 $ rwork(iphi), work(itaup1), work(itaup2),
644 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
648 IF( wantu1 .AND. p .GT. 0 )
THEN
649 CALL zlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
650 CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
651 $ lorgqr, childinfo )
653 IF( wantu2 .AND. m-p .GT. 0 )
THEN
659 CALL zlacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
661 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
662 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
664 IF( wantv1t .AND. q .GT. 0 )
THEN
665 CALL zlacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
666 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
667 $ work(iorglq), lorglq, childinfo )
672 CALL zbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
673 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
674 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
675 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
676 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
677 $ rwork(ibbcsd), lbbcsd, childinfo )
690 CALL zlapmt( .false., p, q, u1, ldu1, iwork )
693 CALL zlapmr( .false., q, q, v1t, ldv1t, iwork )
702 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
703 $ rwork(iphi), work(itaup1), work(itaup2),
704 $ work(itauq1), work(iorbdb), work(iorbdb+m),
705 $ lorbdb-m, childinfo )
709 IF( wantu1 .AND. p .GT. 0 )
THEN
710 CALL zcopy( p, work(iorbdb), 1, u1, 1 )
714 CALL zlacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
716 CALL zungqr( p, p, m-q, u1, ldu1, work(itaup1),
717 $ work(iorgqr), lorgqr, childinfo )
719 IF( wantu2 .AND. m-p .GT. 0 )
THEN
720 CALL zcopy( m-p, work(iorbdb+p), 1, u2, 1 )
724 CALL zlacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
726 CALL zungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
727 $ work(iorgqr), lorgqr, childinfo )
729 IF( wantv1t .AND. q .GT. 0 )
THEN
730 CALL zlacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
731 CALL zlacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
732 $ v1t(m-q+1,m-q+1), ldv1t )
733 CALL zlacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
734 $ v1t(p+1,p+1), ldv1t )
735 CALL zunglq( q, q, q, v1t, ldv1t, work(itauq1),
736 $ work(iorglq), lorglq, childinfo )
741 CALL zbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
742 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
743 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
744 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
745 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
746 $ rwork(ibbcsd), lbbcsd, childinfo )
759 CALL zlapmt( .false., p, p, u1, ldu1, iwork )
762 CALL zlapmr( .false., p, q, v1t, ldv1t, iwork )
subroutine zunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB3
subroutine zbbcsd(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)
ZBBCSD
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB1
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
ZUNBDB4
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zuncsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZUNCSD2BY1
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlapmt(FORWRD, M, N, X, LDX, K)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine zlapmr(FORWRD, M, N, X, LDX, K)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
subroutine zunbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB2