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
326 DOUBLE PRECISION DZNRM2
328 EXTERNAL dznrm2, lsame
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 xerbla(SRNAME, INFO)
XERBLA
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-negative beta.
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.