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
logical function lsame(CA, CB)
LSAME
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
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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.