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 )