LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
                      "other" convention);
          otherwise:  The upper-right block is made nonpositive (the
                      "default" 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.
Date
November 2015
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 289 of file dorbdb.f.

289 *
290 * -- LAPACK computational routine (version 3.6.0) --
291 * -- LAPACK is a software package provided by Univ. of Tennessee, --
292 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
293 * November 2015
294 *
295 * .. Scalar Arguments ..
296  CHARACTER signs, trans
297  INTEGER info, ldx11, ldx12, ldx21, ldx22, lwork, m, p,
298  $ q
299 * ..
300 * .. Array Arguments ..
301  DOUBLE PRECISION phi( * ), theta( * )
302  DOUBLE PRECISION taup1( * ), taup2( * ), tauq1( * ), tauq2( * ),
303  $ work( * ), x11( ldx11, * ), x12( ldx12, * ),
304  $ x21( ldx21, * ), x22( ldx22, * )
305 * ..
306 *
307 * ====================================================================
308 *
309 * .. Parameters ..
310  DOUBLE PRECISION realone
311  parameter ( realone = 1.0d0 )
312  DOUBLE PRECISION one
313  parameter ( one = 1.0d0 )
314 * ..
315 * .. Local Scalars ..
316  LOGICAL colmajor, lquery
317  INTEGER i, lworkmin, lworkopt
318  DOUBLE PRECISION z1, z2, z3, z4
319 * ..
320 * .. External Subroutines ..
321  EXTERNAL daxpy, dlarf, dlarfgp, dscal, xerbla
322 * ..
323 * .. External Functions ..
324  DOUBLE PRECISION dnrm2
325  LOGICAL lsame
326  EXTERNAL dnrm2, lsame
327 * ..
328 * .. Intrinsic Functions
329  INTRINSIC atan2, cos, max, sin
330 * ..
331 * .. Executable Statements ..
332 *
333 * Test input arguments
334 *
335  info = 0
336  colmajor = .NOT. lsame( trans, 'T' )
337  IF( .NOT. lsame( signs, 'O' ) ) THEN
338  z1 = realone
339  z2 = realone
340  z3 = realone
341  z4 = realone
342  ELSE
343  z1 = realone
344  z2 = -realone
345  z3 = realone
346  z4 = -realone
347  END IF
348  lquery = lwork .EQ. -1
349 *
350  IF( m .LT. 0 ) THEN
351  info = -3
352  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
353  info = -4
354  ELSE IF( q .LT. 0 .OR. q .GT. p .OR. q .GT. m-p .OR.
355  $ q .GT. m-q ) THEN
356  info = -5
357  ELSE IF( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
358  info = -7
359  ELSE IF( .NOT.colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
360  info = -7
361  ELSE IF( colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
362  info = -9
363  ELSE IF( .NOT.colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
364  info = -9
365  ELSE IF( colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
366  info = -11
367  ELSE IF( .NOT.colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
368  info = -11
369  ELSE IF( colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
370  info = -13
371  ELSE IF( .NOT.colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
372  info = -13
373  END IF
374 *
375 * Compute workspace
376 *
377  IF( info .EQ. 0 ) THEN
378  lworkopt = m - q
379  lworkmin = m - q
380  work(1) = lworkopt
381  IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
382  info = -21
383  END IF
384  END IF
385  IF( info .NE. 0 ) THEN
386  CALL xerbla( 'xORBDB', -info )
387  RETURN
388  ELSE IF( lquery ) THEN
389  RETURN
390  END IF
391 *
392 * Handle column-major and row-major separately
393 *
394  IF( colmajor ) THEN
395 *
396 * Reduce columns 1, ..., Q of X11, X12, X21, and X22
397 *
398  DO i = 1, q
399 *
400  IF( i .EQ. 1 ) THEN
401  CALL dscal( p-i+1, z1, x11(i,i), 1 )
402  ELSE
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),
405  $ 1, x11(i,i), 1 )
406  END IF
407  IF( i .EQ. 1 ) THEN
408  CALL dscal( m-p-i+1, z2, x21(i,i), 1 )
409  ELSE
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),
412  $ 1, x21(i,i), 1 )
413  END IF
414 *
415  theta(i) = atan2( dnrm2( m-p-i+1, x21(i,i), 1 ),
416  $ dnrm2( p-i+1, x11(i,i), 1 ) )
417 *
418  IF( p .GT. i ) THEN
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) )
422  END IF
423  x11(i,i) = one
424  IF ( m-p .GT. i ) THEN
425  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
426  $ taup2(i) )
427  ELSE IF ( m-p .EQ. i ) THEN
428  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i), 1, taup2(i) )
429  END IF
430  x21(i,i) = one
431 *
432  IF ( q .GT. i ) THEN
433  CALL dlarf( 'L', p-i+1, q-i, x11(i,i), 1, taup1(i),
434  $ x11(i,i+1), ldx11, work )
435  END IF
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 )
439  END IF
440  IF ( q .GT. i ) THEN
441  CALL dlarf( 'L', m-p-i+1, q-i, x21(i,i), 1, taup2(i),
442  $ x21(i,i+1), ldx21, work )
443  END IF
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 )
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), ldx21,
453  $ x11(i,i+1), ldx11 )
454  END IF
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,
457  $ x12(i,i), ldx12 )
458 *
459  IF( i .LT. q )
460  $ phi(i) = atan2( dnrm2( q-i, x11(i,i+1), ldx11 ),
461  $ dnrm2( m-q-i+1, x12(i,i), ldx12 ) )
462 *
463  IF( i .LT. q ) THEN
464  IF ( q-i .EQ. 1 ) THEN
465  CALL dlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
466  $ tauq1(i) )
467  ELSE
468  CALL dlarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
469  $ tauq1(i) )
470  END IF
471  x11(i,i+1) = one
472  END IF
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,
476  $ tauq2(i) )
477  ELSE
478  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
479  $ tauq2(i) )
480  END IF
481  END IF
482  x12(i,i) = one
483 *
484  IF( i .LT. q ) THEN
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 )
489  END IF
490  IF ( p .GT. i ) THEN
491  CALL dlarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
492  $ x12(i+1,i), ldx12, work )
493  END IF
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 )
497  END IF
498 *
499  END DO
500 *
501 * Reduce columns Q + 1, ..., P of X12, X22
502 *
503  DO i = q + 1, p
504 *
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,
508  $ tauq2(i) )
509  ELSE
510  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
511  $ tauq2(i) )
512  END IF
513  x12(i,i) = one
514 *
515  IF ( p .GT. i ) THEN
516  CALL dlarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
517  $ x12(i+1,i), ldx12, work )
518  END IF
519  IF( m-p-q .GE. 1 )
520  $ CALL dlarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
521  $ tauq2(i), x22(q+1,i), ldx22, work )
522 *
523  END DO
524 *
525 * Reduce columns P + 1, ..., M - Q of X12, X22
526 *
527  DO i = 1, m - p - q
528 *
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) )
533  ELSE
534  CALL dlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
535  $ ldx22, tauq2(p+i) )
536  END IF
537  x22(q+i,p+i) = one
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 )
541  END IF
542 *
543  END DO
544 *
545  ELSE
546 *
547 * Reduce columns 1, ..., Q of X11, X12, X21, X22
548 *
549  DO i = 1, q
550 *
551  IF( i .EQ. 1 ) THEN
552  CALL dscal( p-i+1, z1, x11(i,i), ldx11 )
553  ELSE
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 )
557  END IF
558  IF( i .EQ. 1 ) THEN
559  CALL dscal( m-p-i+1, z2, x21(i,i), ldx21 )
560  ELSE
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 )
564  END IF
565 *
566  theta(i) = atan2( dnrm2( m-p-i+1, x21(i,i), ldx21 ),
567  $ dnrm2( p-i+1, x11(i,i), ldx11 ) )
568 *
569  CALL dlarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11, taup1(i) )
570  x11(i,i) = one
571  IF ( i .EQ. m-p ) THEN
572  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
573  $ taup2(i) )
574  ELSE
575  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
576  $ taup2(i) )
577  END IF
578  x21(i,i) = one
579 *
580  IF ( q .GT. i ) THEN
581  CALL dlarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
582  $ x11(i+1,i), ldx11, work )
583  END IF
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 )
587  END IF
588  IF ( q .GT. i ) THEN
589  CALL dlarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
590  $ x21(i+1,i), ldx21, work )
591  END IF
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 )
595  END IF
596 *
597  IF( i .LT. q ) THEN
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,
600  $ x11(i+1,i), 1 )
601  END IF
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,
604  $ x12(i,i), 1 )
605 *
606  IF( i .LT. q )
607  $ phi(i) = atan2( dnrm2( q-i, x11(i+1,i), 1 ),
608  $ dnrm2( m-q-i+1, x12(i,i), 1 ) )
609 *
610  IF( i .LT. q ) THEN
611  IF ( q-i .EQ. 1) THEN
612  CALL dlarfgp( q-i, x11(i+1,i), x11(i+1,i), 1,
613  $ tauq1(i) )
614  ELSE
615  CALL dlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1,
616  $ tauq1(i) )
617  END IF
618  x11(i+1,i) = one
619  END IF
620  IF ( m-q .GT. i ) THEN
621  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
622  $ tauq2(i) )
623  ELSE
624  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i), 1,
625  $ tauq2(i) )
626  END IF
627  x12(i,i) = one
628 *
629  IF( i .LT. q ) THEN
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 )
634  END IF
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 )
640  END IF
641 *
642  END DO
643 *
644 * Reduce columns Q + 1, ..., P of X12, X22
645 *
646  DO i = q + 1, p
647 *
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) )
650  x12(i,i) = one
651 *
652  IF ( p .GT. i ) THEN
653  CALL dlarf( 'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
654  $ x12(i,i+1), ldx12, work )
655  END IF
656  IF( m-p-q .GE. 1 )
657  $ CALL dlarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1, tauq2(i),
658  $ x22(i,q+1), ldx22, work )
659 *
660  END DO
661 *
662 * Reduce columns P + 1, ..., M - Q of X12, X22
663 *
664  DO i = 1, m - p - q
665 *
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,
669  $ tauq2(p+i) )
670  ELSE
671  CALL dlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
672  $ tauq2(p+i) )
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 )
675  END IF
676  x22(p+i,q+i) = one
677 *
678  END DO
679 *
680  END IF
681 *
682  RETURN
683 *
684 * End of DORBDB
685 *
dlarfgp
subroutine dlarfgp(N, ALPHA, X, INCX, TAU)
DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: dlarfgp.f:106
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
dlarf
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dnrm2
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: