286 SUBROUTINE zunbdb( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
287 $ x21, ldx21, x22, ldx22, theta, phi, taup1,
288 $ taup2, tauq1, tauq2, work, lwork, info )
296 CHARACTER signs, trans
297 INTEGER info, ldx11, ldx12, ldx21, ldx22, lwork, m, p,
301 DOUBLE PRECISION phi( * ), theta( * )
302 COMPLEX*16 taup1( * ), taup2( * ), tauq1( * ), tauq2( * ),
303 $ work( * ), x11( ldx11, * ), x12( ldx12, * ),
304 $ x21( ldx21, * ), x22( ldx22, * )
310 DOUBLE PRECISION realone
311 parameter( realone = 1.0d0 )
313 parameter( one = (1.0d0,0.0d0) )
316 LOGICAL colmajor, lquery
317 INTEGER i, lworkmin, lworkopt, pi1, qi1
318 DOUBLE PRECISION z1, z2, z3, z4
331 INTRINSIC atan2, cos, max, min, sin
332 INTRINSIC dcmplx, dconjg
339 colmajor = .NOT.
lsame( trans,
'T' )
340 IF( .NOT.
lsame( signs,
'O' ) )
THEN
351 lquery = lwork .EQ. -1
355 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
357 ELSE IF( q .LT. 0 .OR. q .GT. p .OR. q .GT. m-p .OR.
360 ELSE IF( colmajor .AND. ldx11 .LT. max( 1, p ) )
THEN
362 ELSE IF( .NOT.colmajor .AND. ldx11 .LT. max( 1, q ) )
THEN
364 ELSE IF( colmajor .AND. ldx12 .LT. max( 1, p ) )
THEN
366 ELSE IF( .NOT.colmajor .AND. ldx12 .LT. max( 1, m-q ) )
THEN
368 ELSE IF( colmajor .AND. ldx21 .LT. max( 1, m-p ) )
THEN
370 ELSE IF( .NOT.colmajor .AND. ldx21 .LT. max( 1, q ) )
THEN
372 ELSE IF( colmajor .AND. ldx22 .LT. max( 1, m-p ) )
THEN
374 ELSE IF( .NOT.colmajor .AND. ldx22 .LT. max( 1, m-q ) )
THEN
380 IF( info .EQ. 0 )
THEN
384 IF( lwork .LT. lworkmin .AND. .NOT. lquery )
THEN
388 IF( info .NE. 0 )
THEN
389 CALL
xerbla(
'xORBDB', -info )
391 ELSE IF( lquery )
THEN
404 CALL
zscal( p-i+1, dcmplx( z1, 0.0d0 ), x11(i,i), 1 )
406 CALL
zscal( p-i+1, dcmplx( z1*cos(phi(i-1)), 0.0d0 ),
408 CALL
zaxpy( p-i+1, dcmplx( -z1*z3*z4*sin(phi(i-1)),
409 $ 0.0d0 ), x12(i,i-1), 1, x11(i,i), 1 )
412 CALL
zscal( m-p-i+1, dcmplx( z2, 0.0d0 ), x21(i,i), 1 )
414 CALL
zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)), 0.0d0 ),
416 CALL
zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
417 $ 0.0d0 ), x22(i,i-1), 1, x21(i,i), 1 )
420 theta(i) = atan2(
dznrm2( m-p-i+1, x21(i,i), 1 ),
421 $
dznrm2( p-i+1, x11(i,i), 1 ) )
424 CALL
zlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
425 ELSE IF ( p .EQ. i )
THEN
426 CALL
zlarfgp( p-i+1, x11(i,i), x11(i,i), 1, taup1(i) )
429 IF ( m-p .GT. i )
THEN
430 CALL
zlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
432 ELSE IF ( m-p .EQ. i )
THEN
433 CALL
zlarfgp( m-p-i+1, x21(i,i), x21(i,i), 1,
439 CALL
zlarf(
'L', p-i+1, q-i, x11(i,i), 1,
440 $ dconjg(taup1(i)), x11(i,i+1), ldx11, work )
441 CALL
zlarf(
'L', m-p-i+1, q-i, x21(i,i), 1,
442 $ dconjg(taup2(i)), x21(i,i+1), ldx21, work )
444 IF ( m-q+1 .GT. i )
THEN
445 CALL
zlarf(
'L', p-i+1, m-q-i+1, x11(i,i), 1,
446 $ dconjg(taup1(i)), x12(i,i), ldx12, work )
447 CALL
zlarf(
'L', m-p-i+1, m-q-i+1, x21(i,i), 1,
448 $ dconjg(taup2(i)), x22(i,i), ldx22, work )
452 CALL
zscal( q-i, dcmplx( -z1*z3*sin(theta(i)), 0.0d0 ),
453 $ x11(i,i+1), ldx11 )
454 CALL
zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
455 $ x21(i,i+1), ldx21, x11(i,i+1), ldx11 )
457 CALL
zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)), 0.0d0 ),
459 CALL
zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)), 0.0d0 ),
460 $ x22(i,i), ldx22, x12(i,i), ldx12 )
463 $ phi(i) = atan2(
dznrm2( q-i, x11(i,i+1), ldx11 ),
464 $
dznrm2( m-q-i+1, x12(i,i), ldx12 ) )
467 CALL
zlacgv( q-i, x11(i,i+1), ldx11 )
468 IF ( i .EQ. q-1 )
THEN
469 CALL
zlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
472 CALL
zlarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
477 IF ( m-q+1 .GT. i )
THEN
478 CALL
zlacgv( m-q-i+1, x12(i,i), ldx12 )
479 IF ( m-q .EQ. i )
THEN
480 CALL
zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
483 CALL
zlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
490 CALL
zlarf(
'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
491 $ x11(i+1,i+1), ldx11, work )
492 CALL
zlarf(
'R', m-p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
493 $ x21(i+1,i+1), ldx21, work )
496 CALL
zlarf(
'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
497 $ x12(i+1,i), ldx12, work )
499 IF ( m-p .GT. i )
THEN
500 CALL
zlarf(
'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
501 $ tauq2(i), x22(i+1,i), ldx22, work )
505 $ CALL
zlacgv( q-i, x11(i,i+1), ldx11 )
506 CALL
zlacgv( m-q-i+1, x12(i,i), ldx12 )
514 CALL
zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i),
516 CALL
zlacgv( m-q-i+1, x12(i,i), ldx12 )
517 IF ( i .GE. m-q )
THEN
518 CALL
zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
521 CALL
zlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
527 CALL
zlarf(
'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
528 $ x12(i+1,i), ldx12, work )
531 $ CALL
zlarf(
'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
532 $ tauq2(i), x22(q+1,i), ldx22, work )
534 CALL
zlacgv( m-q-i+1, x12(i,i), ldx12 )
542 CALL
zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
543 $ x22(q+i,p+i), ldx22 )
544 CALL
zlacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
545 CALL
zlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
546 $ ldx22, tauq2(p+i) )
548 CALL
zlarf(
'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i), ldx22,
549 $ tauq2(p+i), x22(q+i+1,p+i), ldx22, work )
551 CALL
zlacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
562 CALL
zscal( p-i+1, dcmplx( z1, 0.0d0 ), x11(i,i),
565 CALL
zscal( p-i+1, dcmplx( z1*cos(phi(i-1)), 0.0d0 ),
567 CALL
zaxpy( p-i+1, dcmplx( -z1*z3*z4*sin(phi(i-1)),
568 $ 0.0d0 ), x12(i-1,i), ldx12, x11(i,i), ldx11 )
571 CALL
zscal( m-p-i+1, dcmplx( z2, 0.0d0 ), x21(i,i),
574 CALL
zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)), 0.0d0 ),
576 CALL
zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
577 $ 0.0d0 ), x22(i-1,i), ldx22, x21(i,i), ldx21 )
580 theta(i) = atan2(
dznrm2( m-p-i+1, x21(i,i), ldx21 ),
581 $
dznrm2( p-i+1, x11(i,i), ldx11 ) )
583 CALL
zlacgv( p-i+1, x11(i,i), ldx11 )
584 CALL
zlacgv( m-p-i+1, x21(i,i), ldx21 )
586 CALL
zlarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11, taup1(i) )
588 IF ( i .EQ. m-p )
THEN
589 CALL
zlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
592 CALL
zlarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
597 CALL
zlarf(
'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
598 $ x11(i+1,i), ldx11, work )
599 CALL
zlarf(
'R', m-q-i+1, p-i+1, x11(i,i), ldx11, taup1(i),
600 $ x12(i,i), ldx12, work )
601 CALL
zlarf(
'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
602 $ x21(i+1,i), ldx21, work )
603 CALL
zlarf(
'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
604 $ taup2(i), x22(i,i), ldx22, work )
606 CALL
zlacgv( p-i+1, x11(i,i), ldx11 )
607 CALL
zlacgv( m-p-i+1, x21(i,i), ldx21 )
610 CALL
zscal( q-i, dcmplx( -z1*z3*sin(theta(i)), 0.0d0 ),
612 CALL
zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
613 $ x21(i+1,i), 1, x11(i+1,i), 1 )
615 CALL
zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)), 0.0d0 ),
617 CALL
zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)), 0.0d0 ),
618 $ x22(i,i), 1, x12(i,i), 1 )
621 $ phi(i) = atan2(
dznrm2( q-i, x11(i+1,i), 1 ),
622 $
dznrm2( m-q-i+1, x12(i,i), 1 ) )
625 CALL
zlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1, tauq1(i) )
628 CALL
zlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
632 CALL
zlarf(
'L', q-i, p-i, x11(i+1,i), 1,
633 $ dconjg(tauq1(i)), x11(i+1,i+1), ldx11, work )
634 CALL
zlarf(
'L', q-i, m-p-i, x11(i+1,i), 1,
635 $ dconjg(tauq1(i)), x21(i+1,i+1), ldx21, work )
637 CALL
zlarf(
'L', m-q-i+1, p-i, x12(i,i), 1,
638 $ dconjg(tauq2(i)), x12(i,i+1), ldx12, work )
639 IF ( m-p .GT. i )
THEN
640 CALL
zlarf(
'L', m-q-i+1, m-p-i, x12(i,i), 1,
641 $ dconjg(tauq2(i)), x22(i,i+1), ldx22, work )
650 CALL
zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i), 1 )
651 CALL
zlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
655 CALL
zlarf(
'L', m-q-i+1, p-i, x12(i,i), 1,
656 $ dconjg(tauq2(i)), x12(i,i+1), ldx12, work )
659 $ CALL
zlarf(
'L', m-q-i+1, m-p-q, x12(i,i), 1,
660 $ dconjg(tauq2(i)), x22(i,q+1), ldx22, work )
668 CALL
zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
670 CALL
zlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
674 IF ( m-p-q .NE. i )
THEN
675 CALL
zlarf(
'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
676 $ dconjg(tauq2(p+i)), x22(p+i,q+i+1), ldx22,
subroutine zunbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
ZUNBDB
subroutine zlarfgp(N, ALPHA, X, INCX, TAU)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negatibe beta.
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
logical function lsame(CA, CB)
LSAME
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
double precision function dznrm2(N, X, INCX)
DZNRM2
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL