LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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
                      "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.
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 284 of file dorbdb.f.

287 *
288 * -- LAPACK computational routine --
289 * -- LAPACK is a software package provided by Univ. of Tennessee, --
290 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
291 *
292 * .. Scalar Arguments ..
293  CHARACTER SIGNS, TRANS
294  INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
295  $ Q
296 * ..
297 * .. Array Arguments ..
298  DOUBLE PRECISION PHI( * ), THETA( * )
299  DOUBLE PRECISION TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
300  $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
301  $ X21( LDX21, * ), X22( LDX22, * )
302 * ..
303 *
304 * ====================================================================
305 *
306 * .. Parameters ..
307  DOUBLE PRECISION REALONE
308  parameter( realone = 1.0d0 )
309  DOUBLE PRECISION ONE
310  parameter( one = 1.0d0 )
311 * ..
312 * .. Local Scalars ..
313  LOGICAL COLMAJOR, LQUERY
314  INTEGER I, LWORKMIN, LWORKOPT
315  DOUBLE PRECISION Z1, Z2, Z3, Z4
316 * ..
317 * .. External Subroutines ..
318  EXTERNAL daxpy, dlarf, dlarfgp, dscal, 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,i-1),
402  $ 1, x11(i,i), 1 )
403  END IF
404  IF( i .EQ. 1 ) THEN
405  CALL dscal( m-p-i+1, z2, x21(i,i), 1 )
406  ELSE
407  CALL dscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i), 1 )
408  CALL daxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i,i-1),
409  $ 1, x21(i,i), 1 )
410  END IF
411 *
412  theta(i) = atan2( dnrm2( m-p-i+1, x21(i,i), 1 ),
413  $ dnrm2( p-i+1, x11(i,i), 1 ) )
414 *
415  IF( p .GT. i ) THEN
416  CALL dlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
417  ELSE IF( p .EQ. i ) THEN
418  CALL dlarfgp( p-i+1, x11(i,i), x11(i,i), 1, taup1(i) )
419  END IF
420  x11(i,i) = one
421  IF ( m-p .GT. i ) THEN
422  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
423  $ taup2(i) )
424  ELSE IF ( m-p .EQ. i ) THEN
425  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i), 1, taup2(i) )
426  END IF
427  x21(i,i) = one
428 *
429  IF ( q .GT. i ) THEN
430  CALL dlarf( 'L', p-i+1, q-i, x11(i,i), 1, taup1(i),
431  $ x11(i,i+1), ldx11, work )
432  END IF
433  IF ( m-q+1 .GT. i ) THEN
434  CALL dlarf( 'L', p-i+1, m-q-i+1, x11(i,i), 1, taup1(i),
435  $ x12(i,i), ldx12, work )
436  END IF
437  IF ( q .GT. i ) THEN
438  CALL dlarf( 'L', m-p-i+1, q-i, x21(i,i), 1, taup2(i),
439  $ x21(i,i+1), ldx21, work )
440  END IF
441  IF ( m-q+1 .GT. i ) THEN
442  CALL dlarf( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1, taup2(i),
443  $ x22(i,i), ldx22, work )
444  END IF
445 *
446  IF( i .LT. q ) THEN
447  CALL dscal( q-i, -z1*z3*sin(theta(i)), x11(i,i+1),
448  $ ldx11 )
449  CALL daxpy( q-i, z2*z3*cos(theta(i)), x21(i,i+1), ldx21,
450  $ x11(i,i+1), ldx11 )
451  END IF
452  CALL dscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i), ldx12 )
453  CALL daxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i), ldx22,
454  $ x12(i,i), ldx12 )
455 *
456  IF( i .LT. q )
457  $ phi(i) = atan2( dnrm2( q-i, x11(i,i+1), ldx11 ),
458  $ dnrm2( m-q-i+1, x12(i,i), ldx12 ) )
459 *
460  IF( i .LT. q ) THEN
461  IF ( q-i .EQ. 1 ) THEN
462  CALL dlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
463  $ tauq1(i) )
464  ELSE
465  CALL dlarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
466  $ tauq1(i) )
467  END IF
468  x11(i,i+1) = one
469  END IF
470  IF ( q+i-1 .LT. m ) THEN
471  IF ( m-q .EQ. i ) THEN
472  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
473  $ tauq2(i) )
474  ELSE
475  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
476  $ tauq2(i) )
477  END IF
478  END IF
479  x12(i,i) = one
480 *
481  IF( i .LT. q ) THEN
482  CALL dlarf( 'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
483  $ x11(i+1,i+1), ldx11, work )
484  CALL dlarf( 'R', m-p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
485  $ x21(i+1,i+1), ldx21, work )
486  END IF
487  IF ( p .GT. i ) THEN
488  CALL dlarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
489  $ x12(i+1,i), ldx12, work )
490  END IF
491  IF ( m-p .GT. i ) THEN
492  CALL dlarf( 'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
493  $ tauq2(i), x22(i+1,i), ldx22, work )
494  END IF
495 *
496  END DO
497 *
498 * Reduce columns Q + 1, ..., P of X12, X22
499 *
500  DO i = q + 1, p
501 *
502  CALL dscal( m-q-i+1, -z1*z4, x12(i,i), ldx12 )
503  IF ( i .GE. m-q ) THEN
504  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
505  $ tauq2(i) )
506  ELSE
507  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
508  $ tauq2(i) )
509  END IF
510  x12(i,i) = one
511 *
512  IF ( p .GT. i ) THEN
513  CALL dlarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
514  $ x12(i+1,i), ldx12, work )
515  END IF
516  IF( m-p-q .GE. 1 )
517  $ CALL dlarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
518  $ tauq2(i), x22(q+1,i), ldx22, work )
519 *
520  END DO
521 *
522 * Reduce columns P + 1, ..., M - Q of X12, X22
523 *
524  DO i = 1, m - p - q
525 *
526  CALL dscal( m-p-q-i+1, z2*z4, x22(q+i,p+i), ldx22 )
527  IF ( i .EQ. m-p-q ) THEN
528  CALL dlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i),
529  $ ldx22, tauq2(p+i) )
530  ELSE
531  CALL dlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
532  $ ldx22, tauq2(p+i) )
533  END IF
534  x22(q+i,p+i) = one
535  IF ( i .LT. m-p-q ) THEN
536  CALL dlarf( 'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i), ldx22,
537  $ tauq2(p+i), x22(q+i+1,p+i), ldx22, work )
538  END IF
539 *
540  END DO
541 *
542  ELSE
543 *
544 * Reduce columns 1, ..., Q of X11, X12, X21, X22
545 *
546  DO i = 1, q
547 *
548  IF( i .EQ. 1 ) THEN
549  CALL dscal( p-i+1, z1, x11(i,i), ldx11 )
550  ELSE
551  CALL dscal( p-i+1, z1*cos(phi(i-1)), x11(i,i), ldx11 )
552  CALL daxpy( p-i+1, -z1*z3*z4*sin(phi(i-1)), x12(i-1,i),
553  $ ldx12, x11(i,i), ldx11 )
554  END IF
555  IF( i .EQ. 1 ) THEN
556  CALL dscal( m-p-i+1, z2, x21(i,i), ldx21 )
557  ELSE
558  CALL dscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i), ldx21 )
559  CALL daxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i-1,i),
560  $ ldx22, x21(i,i), ldx21 )
561  END IF
562 *
563  theta(i) = atan2( dnrm2( m-p-i+1, x21(i,i), ldx21 ),
564  $ dnrm2( p-i+1, x11(i,i), ldx11 ) )
565 *
566  CALL dlarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11, taup1(i) )
567  x11(i,i) = one
568  IF ( i .EQ. m-p ) THEN
569  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
570  $ taup2(i) )
571  ELSE
572  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
573  $ taup2(i) )
574  END IF
575  x21(i,i) = one
576 *
577  IF ( q .GT. i ) THEN
578  CALL dlarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
579  $ x11(i+1,i), ldx11, work )
580  END IF
581  IF ( m-q+1 .GT. i ) THEN
582  CALL dlarf( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11,
583  $ taup1(i), x12(i,i), ldx12, work )
584  END IF
585  IF ( q .GT. i ) THEN
586  CALL dlarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
587  $ x21(i+1,i), ldx21, work )
588  END IF
589  IF ( m-q+1 .GT. i ) THEN
590  CALL dlarf( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
591  $ taup2(i), x22(i,i), ldx22, work )
592  END IF
593 *
594  IF( i .LT. q ) THEN
595  CALL dscal( q-i, -z1*z3*sin(theta(i)), x11(i+1,i), 1 )
596  CALL daxpy( q-i, z2*z3*cos(theta(i)), x21(i+1,i), 1,
597  $ x11(i+1,i), 1 )
598  END IF
599  CALL dscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i), 1 )
600  CALL daxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i), 1,
601  $ x12(i,i), 1 )
602 *
603  IF( i .LT. q )
604  $ phi(i) = atan2( dnrm2( q-i, x11(i+1,i), 1 ),
605  $ dnrm2( m-q-i+1, x12(i,i), 1 ) )
606 *
607  IF( i .LT. q ) THEN
608  IF ( q-i .EQ. 1) THEN
609  CALL dlarfgp( q-i, x11(i+1,i), x11(i+1,i), 1,
610  $ tauq1(i) )
611  ELSE
612  CALL dlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1,
613  $ tauq1(i) )
614  END IF
615  x11(i+1,i) = one
616  END IF
617  IF ( m-q .GT. i ) THEN
618  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
619  $ tauq2(i) )
620  ELSE
621  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i,i), 1,
622  $ tauq2(i) )
623  END IF
624  x12(i,i) = one
625 *
626  IF( i .LT. q ) THEN
627  CALL dlarf( 'L', q-i, p-i, x11(i+1,i), 1, tauq1(i),
628  $ x11(i+1,i+1), ldx11, work )
629  CALL dlarf( 'L', q-i, m-p-i, x11(i+1,i), 1, tauq1(i),
630  $ x21(i+1,i+1), ldx21, work )
631  END IF
632  CALL dlarf( 'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
633  $ x12(i,i+1), ldx12, work )
634  IF ( m-p-i .GT. 0 ) THEN
635  CALL dlarf( 'L', m-q-i+1, m-p-i, x12(i,i), 1, tauq2(i),
636  $ x22(i,i+1), ldx22, work )
637  END IF
638 *
639  END DO
640 *
641 * Reduce columns Q + 1, ..., P of X12, X22
642 *
643  DO i = q + 1, p
644 *
645  CALL dscal( m-q-i+1, -z1*z4, x12(i,i), 1 )
646  CALL dlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
647  x12(i,i) = one
648 *
649  IF ( p .GT. i ) THEN
650  CALL dlarf( 'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
651  $ x12(i,i+1), ldx12, work )
652  END IF
653  IF( m-p-q .GE. 1 )
654  $ CALL dlarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1, tauq2(i),
655  $ x22(i,q+1), ldx22, work )
656 *
657  END DO
658 *
659 * Reduce columns P + 1, ..., M - Q of X12, X22
660 *
661  DO i = 1, m - p - q
662 *
663  CALL dscal( m-p-q-i+1, z2*z4, x22(p+i,q+i), 1 )
664  IF ( m-p-q .EQ. i ) THEN
665  CALL dlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i,q+i), 1,
666  $ tauq2(p+i) )
667  ELSE
668  CALL dlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
669  $ tauq2(p+i) )
670  CALL dlarf( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
671  $ tauq2(p+i), x22(p+i,q+i+1), ldx22, work )
672  END IF
673  x22(p+i,q+i) = one
674 *
675  END DO
676 *
677  END IF
678 *
679  RETURN
680 *
681 * End of DORBDB
682 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:124
subroutine dlarfgp(N, ALPHA, X, INCX, TAU)
DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: dlarfgp.f:104
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: