259 SUBROUTINE zuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
260 $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
261 $ ldv1t, work, lwork, rwork, lrwork, iwork,
270 CHARACTER JOBU1, JOBU2, JOBV1T
271 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
273 INTEGER LRWORK, LRWORKMIN, LRWORKOPT
276 DOUBLE PRECISION RWORK(*)
277 DOUBLE PRECISION THETA(*)
278 COMPLEX*16 U1(ldu1,*), U2(ldu2,*), V1T(ldv1t,*), WORK(*),
279 $ x11(ldx11,*), x21(ldx21,*)
287 parameter( one = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
290 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
291 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
292 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
293 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
294 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
295 $ lworkmin, lworkopt, r
296 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
308 INTRINSIC int, max, min
315 wantu1 = lsame( jobu1,
'Y' )
316 wantu2 = lsame( jobu2,
'Y' )
317 wantv1t = lsame( jobv1t,
'Y' )
318 lquery = lwork .EQ. -1
322 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
324 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
326 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
328 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
330 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
332 ELSE IF( wantu2 .AND. ldu2 .LT. m - p )
THEN
334 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
338 r = min( p, m-p, q, m-q )
373 IF( info .EQ. 0 )
THEN
375 ib11d = iphi + max( 1, r-1 )
376 ib11e = ib11d + max( 1, r )
377 ib12d = ib11e + max( 1, r - 1 )
378 ib12e = ib12d + max( 1, r )
379 ib21d = ib12e + max( 1, r - 1 )
380 ib21e = ib21d + max( 1, r )
381 ib22d = ib21e + max( 1, r - 1 )
382 ib22e = ib22d + max( 1, r )
383 ibbcsd = ib22e + max( 1, r - 1 )
385 itaup2 = itaup1 + max( 1, p )
386 itauq1 = itaup2 + max( 1, m-p )
387 iorbdb = itauq1 + max( 1, q )
388 iorgqr = itauq1 + max( 1, q )
389 iorglq = itauq1 + max( 1, q )
391 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
392 $ 0, 0, work, -1, childinfo )
393 lorbdb = int( work(1) )
394 IF( p .GE. m-p )
THEN
395 CALL zungqr( p, p, q, u1, ldu1, 0, work(1), -1,
397 lorgqrmin = max( 1, p )
398 lorgqropt = int( work(1) )
400 CALL zungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
402 lorgqrmin = max( 1, m-p )
403 lorgqropt = int( work(1) )
405 CALL zunglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
406 $ 0, work(1), -1, childinfo )
407 lorglqmin = max( 1, q-1 )
408 lorglqopt = int( work(1) )
409 CALL zbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
410 $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
411 $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
412 lbbcsd = int( rwork(1) )
413 ELSE IF( r .EQ. p )
THEN
414 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
415 $ 0, 0, work(1), -1, childinfo )
416 lorbdb = int( work(1) )
417 IF( p-1 .GE. m-p )
THEN
418 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
420 lorgqrmin = max( 1, p-1 )
421 lorgqropt = int( work(1) )
423 CALL zungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
425 lorgqrmin = max( 1, m-p )
426 lorgqropt = int( work(1) )
428 CALL zunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
430 lorglqmin = max( 1, q )
431 lorglqopt = int( work(1) )
432 CALL zbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
433 $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
434 $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
435 lbbcsd = int( rwork(1) )
436 ELSE IF( r .EQ. m-p )
THEN
437 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
438 $ 0, 0, work(1), -1, childinfo )
439 lorbdb = int( work(1) )
440 IF( p .GE. m-p-1 )
THEN
441 CALL zungqr( p, p, q, u1, ldu1, 0, work(1), -1,
443 lorgqrmin = max( 1, p )
444 lorgqropt = int( work(1) )
446 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
447 $ work(1), -1, childinfo )
448 lorgqrmin = max( 1, m-p-1 )
449 lorgqropt = int( work(1) )
451 CALL zunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
453 lorglqmin = max( 1, q )
454 lorglqopt = int( work(1) )
455 CALL zbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
456 $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
457 $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
459 lbbcsd = int( rwork(1) )
461 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
462 $ 0, 0, 0, work(1), -1, childinfo )
463 lorbdb = m + int( work(1) )
464 IF( p .GE. m-p )
THEN
465 CALL zungqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
467 lorgqrmin = max( 1, p )
468 lorgqropt = int( work(1) )
470 CALL zungqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
472 lorgqrmin = max( 1, m-p )
473 lorgqropt = int( work(1) )
475 CALL zunglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
477 lorglqmin = max( 1, q )
478 lorglqopt = int( work(1) )
479 CALL zbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
480 $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
481 $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
483 lbbcsd = int( rwork(1) )
485 lrworkmin = ibbcsd+lbbcsd-1
486 lrworkopt = lrworkmin
488 lworkmin = max( iorbdb+lorbdb-1,
489 $ iorgqr+lorgqrmin-1,
490 $ iorglq+lorglqmin-1 )
491 lworkopt = max( iorbdb+lorbdb-1,
492 $ iorgqr+lorgqropt-1,
493 $ iorglq+lorglqopt-1 )
495 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
499 IF( info .NE. 0 )
THEN
500 CALL xerbla(
'ZUNCSD2BY1', -info )
502 ELSE IF( lquery )
THEN
505 lorgqr = lwork-iorgqr+1
506 lorglq = lwork-iorglq+1
517 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
518 $ rwork(iphi), work(itaup1), work(itaup2),
519 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
523 IF( wantu1 .AND. p .GT. 0 )
THEN
524 CALL zlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
525 CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
526 $ lorgqr, childinfo )
528 IF( wantu2 .AND. m-p .GT. 0 )
THEN
529 CALL zlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
530 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
531 $ work(iorgqr), lorgqr, childinfo )
533 IF( wantv1t .AND. q .GT. 0 )
THEN
539 CALL zlacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
541 CALL zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
542 $ work(iorglq), lorglq, childinfo )
547 CALL zbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
548 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
549 $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
550 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
551 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
557 IF( q .GT. 0 .AND. wantu2 )
THEN
559 iwork(i) = m - p - q + i
564 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
566 ELSE IF( r .EQ. p )
THEN
572 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
573 $ rwork(iphi), work(itaup1), work(itaup2),
574 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
578 IF( wantu1 .AND. p .GT. 0 )
THEN
584 CALL zlacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
585 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
586 $ work(iorgqr), lorgqr, childinfo )
588 IF( wantu2 .AND. m-p .GT. 0 )
THEN
589 CALL zlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
590 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
591 $ work(iorgqr), lorgqr, childinfo )
593 IF( wantv1t .AND. q .GT. 0 )
THEN
594 CALL zlacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
595 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
596 $ work(iorglq), lorglq, childinfo )
601 CALL zbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
602 $ rwork(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
603 $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
604 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
605 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
611 IF( q .GT. 0 .AND. wantu2 )
THEN
613 iwork(i) = m - p - q + i
618 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
620 ELSE IF( r .EQ. m-p )
THEN
626 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
627 $ rwork(iphi), work(itaup1), work(itaup2),
628 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
632 IF( wantu1 .AND. p .GT. 0 )
THEN
633 CALL zlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
634 CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
635 $ lorgqr, childinfo )
637 IF( wantu2 .AND. m-p .GT. 0 )
THEN
643 CALL zlacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
645 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
646 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
648 IF( wantv1t .AND. q .GT. 0 )
THEN
649 CALL zlacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
650 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
651 $ work(iorglq), lorglq, childinfo )
656 CALL zbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
657 $ theta, rwork(iphi), 0, 1, v1t, ldv1t, u2, ldu2,
658 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
659 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
660 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
661 $ rwork(ibbcsd), lbbcsd, childinfo )
674 CALL zlapmt( .false., p, q, u1, ldu1, iwork )
677 CALL zlapmr( .false., q, q, v1t, ldv1t, iwork )
686 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
687 $ rwork(iphi), work(itaup1), work(itaup2),
688 $ work(itauq1), work(iorbdb), work(iorbdb+m),
689 $ lorbdb-m, childinfo )
693 IF( wantu1 .AND. p .GT. 0 )
THEN
694 CALL zcopy( p, work(iorbdb), 1, u1, 1 )
698 CALL zlacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
700 CALL zungqr( p, p, m-q, u1, ldu1, work(itaup1),
701 $ work(iorgqr), lorgqr, childinfo )
703 IF( wantu2 .AND. m-p .GT. 0 )
THEN
704 CALL zcopy( m-p, work(iorbdb+p), 1, u2, 1 )
708 CALL zlacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
710 CALL zungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
711 $ work(iorgqr), lorgqr, childinfo )
713 IF( wantv1t .AND. q .GT. 0 )
THEN
714 CALL zlacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
715 CALL zlacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
716 $ v1t(m-q+1,m-q+1), ldv1t )
717 CALL zlacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
718 $ v1t(p+1,p+1), ldv1t )
719 CALL zunglq( q, q, q, v1t, ldv1t, work(itauq1),
720 $ work(iorglq), lorglq, childinfo )
725 CALL zbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
726 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
727 $ ldv1t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
728 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
729 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
743 CALL zlapmt( .false., p, p, u1, ldu1, iwork )
746 CALL zlapmr( .false., p, q, v1t, ldv1t, iwork )
subroutine zunbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
ZUNBDB4
subroutine zlapmt(FORWRD, M, N, X, LDX, K)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zunbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB1
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zunbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB2
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB3
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
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 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 zlapmr(FORWRD, M, N, X, LDX, K)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.