286 SUBROUTINE dorbdb( 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 DOUBLE PRECISION 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 )
316 LOGICAL colmajor, lquery
317 INTEGER i, lworkmin, lworkopt
318 DOUBLE PRECISION z1, z2, z3, z4
324 DOUBLE PRECISION dnrm2
329 INTRINSIC atan2, cos, max, sin
336 colmajor = .NOT.
lsame( trans,
'T' )
337 IF( .NOT.
lsame( signs,
'O' ) )
THEN
348 lquery = lwork .EQ. -1
352 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
354 ELSE IF( q .LT. 0 .OR. q .GT. p .OR. q .GT. m-p .OR.
357 ELSE IF( colmajor .AND. ldx11 .LT. max( 1, p ) )
THEN
359 ELSE IF( .NOT.colmajor .AND. ldx11 .LT. max( 1, q ) )
THEN
361 ELSE IF( colmajor .AND. ldx12 .LT. max( 1, p ) )
THEN
363 ELSE IF( .NOT.colmajor .AND. ldx12 .LT. max( 1, m-q ) )
THEN
365 ELSE IF( colmajor .AND. ldx21 .LT. max( 1, m-p ) )
THEN
367 ELSE IF( .NOT.colmajor .AND. ldx21 .LT. max( 1, q ) )
THEN
369 ELSE IF( colmajor .AND. ldx22 .LT. max( 1, m-p ) )
THEN
371 ELSE IF( .NOT.colmajor .AND. ldx22 .LT. max( 1, m-q ) )
THEN
377 IF( info .EQ. 0 )
THEN
381 IF( lwork .LT. lworkmin .AND. .NOT. lquery )
THEN
385 IF( info .NE. 0 )
THEN
386 CALL
xerbla(
'xORBDB', -info )
388 ELSE IF( lquery )
THEN
401 CALL
dscal( p-i+1, z1, x11(i,i), 1 )
403 CALL
dscal( p-i+1, z1*cos(phi(i-1)), x11(i,i), 1 )
404 CALL
daxpy( p-i+1, -z1*z3*z4*sin(phi(i-1)), x12(i,i-1),
408 CALL
dscal( m-p-i+1, z2, x21(i,i), 1 )
410 CALL
dscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i), 1 )
411 CALL
daxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i,i-1),
415 theta(i) = atan2(
dnrm2( m-p-i+1, x21(i,i), 1 ),
416 $
dnrm2( p-i+1, x11(i,i), 1 ) )
419 CALL
dlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
420 ELSE IF( p .EQ. i )
THEN
421 CALL
dlarfgp( p-i+1, x11(i,i), x11(i,i), 1, taup1(i) )
424 IF ( m-p .GT. i )
THEN
425 CALL
dlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
427 ELSE IF ( m-p .EQ. i )
THEN
428 CALL
dlarfgp( m-p-i+1, x21(i,i), x21(i,i), 1, taup2(i) )
433 CALL
dlarf(
'L', p-i+1, q-i, x11(i,i), 1, taup1(i),
434 $ x11(i,i+1), ldx11, work )
436 IF ( m-q+1 .GT. i )
THEN
437 CALL
dlarf(
'L', p-i+1, m-q-i+1, x11(i,i), 1, taup1(i),
438 $ x12(i,i), ldx12, work )
441 CALL
dlarf(
'L', m-p-i+1, q-i, x21(i,i), 1, taup2(i),
442 $ x21(i,i+1), ldx21, work )
444 IF ( m-q+1 .GT. i )
THEN
445 CALL
dlarf(
'L', m-p-i+1, m-q-i+1, x21(i,i), 1, taup2(i),
446 $ x22(i,i), ldx22, work )
450 CALL
dscal( q-i, -z1*z3*sin(theta(i)), x11(i,i+1),
452 CALL
daxpy( q-i, z2*z3*cos(theta(i)), x21(i,i+1), ldx21,
453 $ x11(i,i+1), ldx11 )
455 CALL
dscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i), ldx12 )
456 CALL
daxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i), ldx22,
460 $ phi(i) = atan2(
dnrm2( q-i, x11(i,i+1), ldx11 ),
461 $
dnrm2( m-q-i+1, x12(i,i), ldx12 ) )
464 IF ( q-i .EQ. 1 )
THEN
465 CALL
dlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
468 CALL
dlarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
473 IF ( q+i-1 .LT. m )
THEN
474 IF ( m-q .EQ. i )
THEN
475 CALL
dlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
478 CALL
dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
485 CALL
dlarf(
'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
486 $ x11(i+1,i+1), ldx11, work )
487 CALL
dlarf(
'R', m-p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
488 $ x21(i+1,i+1), ldx21, work )
491 CALL
dlarf(
'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
492 $ x12(i+1,i), ldx12, work )
494 IF ( m-p .GT. i )
THEN
495 CALL
dlarf(
'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
496 $ tauq2(i), x22(i+1,i), ldx22, work )
505 CALL
dscal( m-q-i+1, -z1*z4, x12(i,i), ldx12 )
506 IF ( i .GE. m-q )
THEN
507 CALL
dlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
510 CALL
dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
516 CALL
dlarf(
'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
517 $ x12(i+1,i), ldx12, work )
520 $ CALL
dlarf(
'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
521 $ tauq2(i), x22(q+1,i), ldx22, work )
529 CALL
dscal( m-p-q-i+1, z2*z4, x22(q+i,p+i), ldx22 )
530 IF ( i .EQ. m-p-q )
THEN
531 CALL
dlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i),
532 $ ldx22, tauq2(p+i) )
534 CALL
dlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
535 $ ldx22, tauq2(p+i) )
538 IF ( i .LT. m-p-q )
THEN
539 CALL
dlarf(
'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i), ldx22,
540 $ tauq2(p+i), x22(q+i+1,p+i), ldx22, work )
552 CALL
dscal( p-i+1, z1, x11(i,i), ldx11 )
554 CALL
dscal( p-i+1, z1*cos(phi(i-1)), x11(i,i), ldx11 )
555 CALL
daxpy( p-i+1, -z1*z3*z4*sin(phi(i-1)), x12(i-1,i),
556 $ ldx12, x11(i,i), ldx11 )
559 CALL
dscal( m-p-i+1, z2, x21(i,i), ldx21 )
561 CALL
dscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i), ldx21 )
562 CALL
daxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i-1,i),
563 $ ldx22, x21(i,i), ldx21 )
566 theta(i) = atan2(
dnrm2( m-p-i+1, x21(i,i), ldx21 ),
567 $
dnrm2( p-i+1, x11(i,i), ldx11 ) )
569 CALL
dlarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11, taup1(i) )
571 IF ( i .EQ. m-p )
THEN
572 CALL
dlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
575 CALL
dlarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
581 CALL
dlarf(
'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
582 $ x11(i+1,i), ldx11, work )
584 IF ( m-q+1 .GT. i )
THEN
585 CALL
dlarf(
'R', m-q-i+1, p-i+1, x11(i,i), ldx11,
586 $ taup1(i), x12(i,i), ldx12, work )
589 CALL
dlarf(
'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
590 $ x21(i+1,i), ldx21, work )
592 IF ( m-q+1 .GT. i )
THEN
593 CALL
dlarf(
'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
594 $ taup2(i), x22(i,i), ldx22, work )
598 CALL
dscal( q-i, -z1*z3*sin(theta(i)), x11(i+1,i), 1 )
599 CALL
daxpy( q-i, z2*z3*cos(theta(i)), x21(i+1,i), 1,
602 CALL
dscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i), 1 )
603 CALL
daxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i), 1,
607 $ phi(i) = atan2(
dnrm2( q-i, x11(i+1,i), 1 ),
608 $
dnrm2( m-q-i+1, x12(i,i), 1 ) )
611 IF ( q-i .EQ. 1)
THEN
612 CALL
dlarfgp( q-i, x11(i+1,i), x11(i+1,i), 1,
615 CALL
dlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1,
620 IF ( m-q .GT. i )
THEN
621 CALL
dlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
624 CALL
dlarfgp( m-q-i+1, x12(i,i), x12(i,i), 1,
630 CALL
dlarf(
'L', q-i, p-i, x11(i+1,i), 1, tauq1(i),
631 $ x11(i+1,i+1), ldx11, work )
632 CALL
dlarf(
'L', q-i, m-p-i, x11(i+1,i), 1, tauq1(i),
633 $ x21(i+1,i+1), ldx21, work )
635 CALL
dlarf(
'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
636 $ x12(i,i+1), ldx12, work )
637 IF ( m-p-i .GT. 0 )
THEN
638 CALL
dlarf(
'L', m-q-i+1, m-p-i, x12(i,i), 1, tauq2(i),
639 $ x22(i,i+1), ldx22, work )
648 CALL
dscal( m-q-i+1, -z1*z4, x12(i,i), 1 )
649 CALL
dlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
653 CALL
dlarf(
'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
654 $ x12(i,i+1), ldx12, work )
657 $ CALL
dlarf(
'L', m-q-i+1, m-p-q, x12(i,i), 1, tauq2(i),
658 $ x22(i,q+1), ldx22, work )
666 CALL
dscal( m-p-q-i+1, z2*z4, x22(p+i,q+i), 1 )
667 IF ( m-p-q .EQ. i )
THEN
668 CALL
dlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i,q+i), 1,
671 CALL
dlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
673 CALL
dlarf(
'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
674 $ tauq2(p+i), x22(p+i,q+i+1), ldx22, work )
LOGICAL function lsame(CA, CB)
LSAME
subroutine dlarfgp(N, ALPHA, X, INCX, TAU)
DLARFGP generates an elementary reflector (Householder matrix) with non-negatibe beta.
subroutine dorbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
DORBDB
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
DOUBLE PRECISION function dnrm2(N, X, INCX)
DNRM2