LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dorbdb()

subroutine dorbdb ( character trans,
character signs,
integer m,
integer p,
integer q,
double precision, dimension( ldx11, * ) x11,
integer ldx11,
double precision, dimension( ldx12, * ) x12,
integer ldx12,
double precision, dimension( ldx21, * ) x21,
integer ldx21,
double precision, dimension( ldx22, * ) x22,
integer ldx22,
double precision, dimension( * ) theta,
double precision, dimension( * ) phi,
double precision, dimension( * ) taup1,
double precision, dimension( * ) taup2,
double precision, dimension( * ) tauq1,
double precision, dimension( * ) tauq2,
double precision, dimension( * ) work,
integer lwork,
integer info )

DORBDB

Download DORBDB + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DORBDB simultaneously bidiagonalizes the blocks of an M-by-M
!> partitioned orthogonal matrix X:
!>
!>                                 [ B11 | B12 0  0 ]
!>     [ X11 | X12 ]   [ P1 |    ] [  0  |  0 -I  0 ] [ Q1 |    ]**T
!> X = [-----------] = [---------] [----------------] [---------]   .
!>     [ X21 | X22 ]   [    | P2 ] [ B21 | B22 0  0 ] [    | Q2 ]
!>                                 [  0  |  0  0  I ]
!>
!> X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
!> not the case, then X must be transposed and/or permuted. This can be
!> done in constant time using the TRANS and SIGNS options. See DORCSD
!> for details.)
!>
!> The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
!> (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
!> represented implicitly by Householder vectors.
!>
!> B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
!> implicitly by angles THETA, PHI.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER
!>          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
!>                      order;
!>          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
!>                      major order.
!>
[in]SIGNS
!>          SIGNS is CHARACTER
!>          = 'O':      The lower-left block is made nonpositive (the
!>                       convention);
!>          otherwise:  The upper-right block is made nonpositive (the
!>                       convention).
!>
[in]M
!>          M is INTEGER
!>          The number of rows and columns in X.
!>
[in]P
!>          P is INTEGER
!>          The number of rows in X11 and X12. 0 <= P <= M.
!>
[in]Q
!>          Q is INTEGER
!>          The number of columns in X11 and X21. 0 <= Q <=
!>          MIN(P,M-P,M-Q).
!>
[in,out]X11
!>          X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
!>          On entry, the top-left block of the orthogonal matrix to be
!>          reduced. On exit, the form depends on TRANS:
!>          If TRANS = 'N', then
!>             the columns of tril(X11) specify reflectors for P1,
!>             the rows of triu(X11,1) specify reflectors for Q1;
!>          else TRANS = 'T', and
!>             the rows of triu(X11) specify reflectors for P1,
!>             the columns of tril(X11,-1) specify reflectors for Q1.
!>
[in]LDX11
!>          LDX11 is INTEGER
!>          The leading dimension of X11. If TRANS = 'N', then LDX11 >=
!>          P; else LDX11 >= Q.
!>
[in,out]X12
!>          X12 is DOUBLE PRECISION array, dimension (LDX12,M-Q)
!>          On entry, the top-right block of the orthogonal matrix to
!>          be reduced. On exit, the form depends on TRANS:
!>          If TRANS = 'N', then
!>             the rows of triu(X12) specify the first P reflectors for
!>             Q2;
!>          else TRANS = 'T', and
!>             the columns of tril(X12) specify the first P reflectors
!>             for Q2.
!>
[in]LDX12
!>          LDX12 is INTEGER
!>          The leading dimension of X12. If TRANS = 'N', then LDX12 >=
!>          P; else LDX11 >= M-Q.
!>
[in,out]X21
!>          X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
!>          On entry, the bottom-left block of the orthogonal matrix to
!>          be reduced. On exit, the form depends on TRANS:
!>          If TRANS = 'N', then
!>             the columns of tril(X21) specify reflectors for P2;
!>          else TRANS = 'T', and
!>             the rows of triu(X21) specify reflectors for P2.
!>
[in]LDX21
!>          LDX21 is INTEGER
!>          The leading dimension of X21. If TRANS = 'N', then LDX21 >=
!>          M-P; else LDX21 >= Q.
!>
[in,out]X22
!>          X22 is DOUBLE PRECISION array, dimension (LDX22,M-Q)
!>          On entry, the bottom-right block of the orthogonal matrix to
!>          be reduced. On exit, the form depends on TRANS:
!>          If TRANS = 'N', then
!>             the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
!>             M-P-Q reflectors for Q2,
!>          else TRANS = 'T', and
!>             the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
!>             M-P-Q reflectors for P2.
!>
[in]LDX22
!>          LDX22 is INTEGER
!>          The leading dimension of X22. If TRANS = 'N', then LDX22 >=
!>          M-P; else LDX22 >= M-Q.
!>
[out]THETA
!>          THETA is DOUBLE PRECISION array, dimension (Q)
!>          The entries of the bidiagonal blocks B11, B12, B21, B22 can
!>          be computed from the angles THETA and PHI. See Further
!>          Details.
!>
[out]PHI
!>          PHI is DOUBLE PRECISION array, dimension (Q-1)
!>          The entries of the bidiagonal blocks B11, B12, B21, B22 can
!>          be computed from the angles THETA and PHI. See Further
!>          Details.
!>
[out]TAUP1
!>          TAUP1 is DOUBLE PRECISION array, dimension (P)
!>          The scalar factors of the elementary reflectors that define
!>          P1.
!>
[out]TAUP2
!>          TAUP2 is DOUBLE PRECISION array, dimension (M-P)
!>          The scalar factors of the elementary reflectors that define
!>          P2.
!>
[out]TAUQ1
!>          TAUQ1 is DOUBLE PRECISION array, dimension (Q)
!>          The scalar factors of the elementary reflectors that define
!>          Q1.
!>
[out]TAUQ2
!>          TAUQ2 is DOUBLE PRECISION array, dimension (M-Q)
!>          The scalar factors of the elementary reflectors that define
!>          Q2.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (LWORK)
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK. LWORK >= M-Q.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The bidiagonal blocks B11, B12, B21, and B22 are represented
!>  implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
!>  PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
!>  lower bidiagonal. Every entry in each bidiagonal band is a product
!>  of a sine or cosine of a THETA with a sine or cosine of a PHI. See
!>  [1] or DORCSD for details.
!>
!>  P1, P2, Q1, and Q2 are represented as products of elementary
!>  reflectors. See DORCSD for details on generating P1, P2, Q1, and Q2
!>  using DORGQR and DORGLQ.
!>
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

Definition at line 282 of file dorbdb.f.

286*
287* -- LAPACK computational routine --
288* -- LAPACK is a software package provided by Univ. of Tennessee, --
289* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290*
291* .. Scalar Arguments ..
292 CHARACTER SIGNS, TRANS
293 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
294 $ Q
295* ..
296* .. Array Arguments ..
297 DOUBLE PRECISION PHI( * ), THETA( * )
298 DOUBLE PRECISION TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
299 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
300 $ X21( LDX21, * ), X22( LDX22, * )
301* ..
302*
303* ====================================================================
304*
305* .. Parameters ..
306 DOUBLE PRECISION REALONE
307 parameter( realone = 1.0d0 )
308 DOUBLE PRECISION ONE
309 parameter( one = 1.0d0 )
310* ..
311* .. Local Scalars ..
312 LOGICAL COLMAJOR, LQUERY
313 INTEGER I, LWORKMIN, LWORKOPT
314 DOUBLE PRECISION Z1, Z2, Z3, Z4
315* ..
316* .. External Subroutines ..
317 EXTERNAL daxpy, dlarf1f, dlarfgp, dscal,
318 $ xerbla
319* ..
320* .. External Functions ..
321 DOUBLE PRECISION DNRM2
322 LOGICAL LSAME
323 EXTERNAL dnrm2, lsame
324* ..
325* .. Intrinsic Functions
326 INTRINSIC atan2, cos, max, sin
327* ..
328* .. Executable Statements ..
329*
330* Test input arguments
331*
332 info = 0
333 colmajor = .NOT. lsame( trans, 'T' )
334 IF( .NOT. lsame( signs, 'O' ) ) THEN
335 z1 = realone
336 z2 = realone
337 z3 = realone
338 z4 = realone
339 ELSE
340 z1 = realone
341 z2 = -realone
342 z3 = realone
343 z4 = -realone
344 END IF
345 lquery = lwork .EQ. -1
346*
347 IF( m .LT. 0 ) THEN
348 info = -3
349 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
350 info = -4
351 ELSE IF( q .LT. 0 .OR. q .GT. p .OR. q .GT. m-p .OR.
352 $ q .GT. m-q ) THEN
353 info = -5
354 ELSE IF( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
355 info = -7
356 ELSE IF( .NOT.colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
357 info = -7
358 ELSE IF( colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
359 info = -9
360 ELSE IF( .NOT.colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
361 info = -9
362 ELSE IF( colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
363 info = -11
364 ELSE IF( .NOT.colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
365 info = -11
366 ELSE IF( colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
367 info = -13
368 ELSE IF( .NOT.colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
369 info = -13
370 END IF
371*
372* Compute workspace
373*
374 IF( info .EQ. 0 ) THEN
375 lworkopt = m - q
376 lworkmin = m - q
377 work(1) = lworkopt
378 IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
379 info = -21
380 END IF
381 END IF
382 IF( info .NE. 0 ) THEN
383 CALL xerbla( 'xORBDB', -info )
384 RETURN
385 ELSE IF( lquery ) THEN
386 RETURN
387 END IF
388*
389* Handle column-major and row-major separately
390*
391 IF( colmajor ) THEN
392*
393* Reduce columns 1, ..., Q of X11, X12, X21, and X22
394*
395 DO i = 1, q
396*
397 IF( i .EQ. 1 ) THEN
398 CALL dscal( p-i+1, z1, x11(i,i), 1 )
399 ELSE
400 CALL dscal( p-i+1, z1*cos(phi(i-1)), x11(i,i), 1 )
401 CALL daxpy( p-i+1, -z1*z3*z4*sin(phi(i-1)), x12(i,
402 $ i-1),
403 $ 1, x11(i,i), 1 )
404 END IF
405 IF( i .EQ. 1 ) THEN
406 CALL dscal( m-p-i+1, z2, x21(i,i), 1 )
407 ELSE
408 CALL dscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i), 1 )
409 CALL daxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i,
410 $ i-1),
411 $ 1, x21(i,i), 1 )
412 END IF
413*
414 theta(i) = atan2( dnrm2( m-p-i+1, x21(i,i), 1 ),
415 $ dnrm2( p-i+1, x11(i,i), 1 ) )
416*
417 IF( p .GT. i ) THEN
418 CALL dlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1,
419 $ taup1(i) )
420 ELSE IF( p .EQ. i ) THEN
421 CALL dlarfgp( p-i+1, x11(i,i), x11(i,i), 1, taup1(i) )
422 END IF
423 IF ( m-p .GT. i ) THEN
424 CALL dlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
425 $ taup2(i) )
426 ELSE IF ( m-p .EQ. i ) THEN
427 CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i), 1,
428 $ taup2(i) )
429 END IF
430*
431 IF ( q .GT. i ) THEN
432 CALL dlarf1f( 'L', p-i+1, q-i, x11(i,i), 1, taup1(i),
433 $ x11(i,i+1), ldx11, work )
434 END IF
435 IF ( m-q+1 .GT. i ) THEN
436 CALL dlarf1f( 'L', p-i+1, m-q-i+1, x11(i,i), 1,
437 $ taup1(i),
438 $ x12(i,i), ldx12, work )
439 END IF
440 IF ( q .GT. i ) THEN
441 CALL dlarf1f( 'L', m-p-i+1, q-i, x21(i,i), 1,
442 $ taup2(i), x21(i,i+1), ldx21, work )
443 END IF
444 IF ( m-q+1 .GT. i ) THEN
445 CALL dlarf1f( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1,
446 $ taup2(i), x22(i,i), ldx22, work )
447 END IF
448*
449 IF( i .LT. q ) THEN
450 CALL dscal( q-i, -z1*z3*sin(theta(i)), x11(i,i+1),
451 $ ldx11 )
452 CALL daxpy( q-i, z2*z3*cos(theta(i)), x21(i,i+1),
453 $ ldx21,
454 $ x11(i,i+1), ldx11 )
455 END IF
456 CALL dscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i),
457 $ ldx12 )
458 CALL daxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i),
459 $ ldx22,
460 $ x12(i,i), ldx12 )
461*
462 IF( i .LT. q )
463 $ phi(i) = atan2( dnrm2( q-i, x11(i,i+1), ldx11 ),
464 $ dnrm2( m-q-i+1, x12(i,i), ldx12 ) )
465*
466 IF( i .LT. q ) THEN
467 IF ( q-i .EQ. 1 ) THEN
468 CALL dlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
469 $ tauq1(i) )
470 ELSE
471 CALL dlarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
472 $ tauq1(i) )
473 END IF
474 END IF
475 IF ( q+i-1 .LT. m ) THEN
476 IF ( m-q .EQ. i ) THEN
477 CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
478 $ tauq2(i) )
479 ELSE
480 CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
481 $ tauq2(i) )
482 END IF
483 END IF
484*
485 IF( i .LT. q ) THEN
486 CALL dlarf1f( 'R', p-i, q-i, x11(i,i+1), ldx11,
487 $ tauq1(i),
488 $ x11(i+1,i+1), ldx11, work )
489 CALL dlarf1f( 'R', m-p-i, q-i, x11(i,i+1), ldx11,
490 $ tauq1(i),
491 $ x21(i+1,i+1), ldx21, work )
492 END IF
493 IF ( p .GT. i ) THEN
494 CALL dlarf1f( 'R', p-i, m-q-i+1, x12(i,i), ldx12,
495 $ tauq2(i),
496 $ x12(i+1,i), ldx12, work )
497 END IF
498 IF ( m-p .GT. i ) THEN
499 CALL dlarf1f( 'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
500 $ tauq2(i), x22(i+1,i), ldx22, work )
501 END IF
502*
503 END DO
504*
505* Reduce columns Q + 1, ..., P of X12, X22
506*
507 DO i = q + 1, p
508*
509 CALL dscal( m-q-i+1, -z1*z4, x12(i,i), ldx12 )
510 IF ( i .GE. m-q ) THEN
511 CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
512 $ tauq2(i) )
513 ELSE
514 CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
515 $ tauq2(i) )
516 END IF
517*
518 IF ( p .GT. i ) THEN
519 CALL dlarf1f( 'R', p-i, m-q-i+1, x12(i,i), ldx12,
520 $ tauq2(i),
521 $ x12(i+1,i), ldx12, work )
522 END IF
523 IF( m-p-q .GE. 1 )
524 $ CALL dlarf1f( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
525 $ tauq2(i), x22(q+1,i), ldx22, work )
526*
527 END DO
528*
529* Reduce columns P + 1, ..., M - Q of X12, X22
530*
531 DO i = 1, m - p - q
532*
533 CALL dscal( m-p-q-i+1, z2*z4, x22(q+i,p+i), ldx22 )
534 IF ( i .EQ. m-p-q ) THEN
535 CALL dlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i),
536 $ ldx22, tauq2(p+i) )
537 ELSE
538 CALL dlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
539 $ ldx22, tauq2(p+i) )
540 END IF
541 IF ( i .LT. m-p-q ) THEN
542 CALL dlarf1f( 'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i),
543 $ ldx22,
544 $ tauq2(p+i), x22(q+i+1,p+i), ldx22, work )
545 END IF
546*
547 END DO
548*
549 ELSE
550*
551* Reduce columns 1, ..., Q of X11, X12, X21, X22
552*
553 DO i = 1, q
554*
555 IF( i .EQ. 1 ) THEN
556 CALL dscal( p-i+1, z1, x11(i,i), ldx11 )
557 ELSE
558 CALL dscal( p-i+1, z1*cos(phi(i-1)), x11(i,i), ldx11 )
559 CALL daxpy( p-i+1, -z1*z3*z4*sin(phi(i-1)), x12(i-1,
560 $ i),
561 $ ldx12, x11(i,i), ldx11 )
562 END IF
563 IF( i .EQ. 1 ) THEN
564 CALL dscal( m-p-i+1, z2, x21(i,i), ldx21 )
565 ELSE
566 CALL dscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i),
567 $ ldx21 )
568 CALL daxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i-1,
569 $ i),
570 $ ldx22, x21(i,i), ldx21 )
571 END IF
572*
573 theta(i) = atan2( dnrm2( m-p-i+1, x21(i,i), ldx21 ),
574 $ dnrm2( p-i+1, x11(i,i), ldx11 ) )
575*
576 CALL dlarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11,
577 $ taup1(i) )
578 IF ( i .EQ. m-p ) THEN
579 CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
580 $ taup2(i) )
581 ELSE
582 CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
583 $ taup2(i) )
584 END IF
585*
586 IF ( q .GT. i ) THEN
587 CALL dlarf1f( 'R', q-i, p-i+1, x11(i,i), ldx11,
588 $ taup1(i),
589 $ x11(i+1,i), ldx11, work )
590 END IF
591 IF ( m-q+1 .GT. i ) THEN
592 CALL dlarf1f( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11,
593 $ taup1(i), x12(i,i), ldx12, work )
594 END IF
595 IF ( q .GT. i ) THEN
596 CALL dlarf1f( 'R', q-i, m-p-i+1, x21(i,i), ldx21,
597 $ taup2(i),
598 $ x21(i+1,i), ldx21, work )
599 END IF
600 IF ( m-q+1 .GT. i ) THEN
601 CALL dlarf1f( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
602 $ taup2(i), x22(i,i), ldx22, work )
603 END IF
604*
605 IF( i .LT. q ) THEN
606 CALL dscal( q-i, -z1*z3*sin(theta(i)), x11(i+1,i), 1 )
607 CALL daxpy( q-i, z2*z3*cos(theta(i)), x21(i+1,i), 1,
608 $ x11(i+1,i), 1 )
609 END IF
610 CALL dscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i), 1 )
611 CALL daxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i), 1,
612 $ x12(i,i), 1 )
613*
614 IF( i .LT. q )
615 $ phi(i) = atan2( dnrm2( q-i, x11(i+1,i), 1 ),
616 $ dnrm2( m-q-i+1, x12(i,i), 1 ) )
617*
618 IF( i .LT. q ) THEN
619 IF ( q-i .EQ. 1) THEN
620 CALL dlarfgp( q-i, x11(i+1,i), x11(i+1,i), 1,
621 $ tauq1(i) )
622 ELSE
623 CALL dlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1,
624 $ tauq1(i) )
625 END IF
626 END IF
627 IF ( m-q .GT. i ) THEN
628 CALL dlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
629 $ tauq2(i) )
630 ELSE
631 CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i), 1,
632 $ tauq2(i) )
633 END IF
634*
635 IF( i .LT. q ) THEN
636 CALL dlarf1f( 'L', q-i, p-i, x11(i+1,i), 1, tauq1(i),
637 $ x11(i+1,i+1), ldx11, work )
638 CALL dlarf1f( 'L', q-i, m-p-i, x11(i+1,i), 1,
639 $ tauq1(i), x21(i+1,i+1), ldx21, work )
640 END IF
641 CALL dlarf1f( 'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
642 $ x12(i,i+1), ldx12, work )
643 IF ( m-p-i .GT. 0 ) THEN
644 CALL dlarf1f( 'L', m-q-i+1, m-p-i, x12(i,i), 1,
645 $ tauq2(i), x22(i,i+1), ldx22, work )
646 END IF
647*
648 END DO
649*
650* Reduce columns Q + 1, ..., P of X12, X22
651*
652 DO i = q + 1, p
653*
654 CALL dscal( m-q-i+1, -z1*z4, x12(i,i), 1 )
655 CALL dlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
656 $ tauq2(i) )
657*
658 IF ( p .GT. i ) THEN
659 CALL dlarf1f( 'L', m-q-i+1, p-i, x12(i,i), 1,
660 $ tauq2(i), x12(i,i+1), ldx12, work )
661 END IF
662 IF( m-p-q .GE. 1 )
663 $ CALL dlarf1f( 'L', m-q-i+1, m-p-q, x12(i,i), 1,
664 $ tauq2(i), x22(i,q+1), ldx22, work )
665*
666 END DO
667*
668* Reduce columns P + 1, ..., M - Q of X12, X22
669*
670 DO i = 1, m - p - q
671*
672 CALL dscal( m-p-q-i+1, z2*z4, x22(p+i,q+i), 1 )
673 IF ( m-p-q .EQ. i ) THEN
674 CALL dlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i,q+i),
675 $ 1,
676 $ tauq2(p+i) )
677 ELSE
678 CALL dlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i),
679 $ 1,
680 $ tauq2(p+i) )
681 CALL dlarf1f( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i),
682 $ 1, tauq2(p+i), x22(p+i,q+i+1), ldx22,
683 $ work )
684 END IF
685*
686 END DO
687*
688 END IF
689*
690 RETURN
691*
692* End of DORBDB
693*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dlarf1f
subroutine dlarf1f(side, m, n, v, incv, tau, c, ldc, work)
DLARF1F applies an elementary reflector to a general rectangular
Definition dlarf1f.f:157
dlarfgp
subroutine dlarfgp(n, alpha, x, incx, tau)
DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition dlarfgp.f:102
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: