Annotation of rpl/lapack/lapack/dgebrd.f, revision 1.6

1.1       bertrand    1:       SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
                      2:      $                   INFO )
                      3: *
                      4: *  -- LAPACK routine (version 3.2) --
                      5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      7: *     November 2006
                      8: *
                      9: *     .. Scalar Arguments ..
                     10:       INTEGER            INFO, LDA, LWORK, M, N
                     11: *     ..
                     12: *     .. Array Arguments ..
                     13:       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAUP( * ),
                     14:      $                   TAUQ( * ), WORK( * )
                     15: *     ..
                     16: *
                     17: *  Purpose
                     18: *  =======
                     19: *
                     20: *  DGEBRD reduces a general real M-by-N matrix A to upper or lower
                     21: *  bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
                     22: *
                     23: *  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
                     24: *
                     25: *  Arguments
                     26: *  =========
                     27: *
                     28: *  M       (input) INTEGER
                     29: *          The number of rows in the matrix A.  M >= 0.
                     30: *
                     31: *  N       (input) INTEGER
                     32: *          The number of columns in the matrix A.  N >= 0.
                     33: *
                     34: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
                     35: *          On entry, the M-by-N general matrix to be reduced.
                     36: *          On exit,
                     37: *          if m >= n, the diagonal and the first superdiagonal are
                     38: *            overwritten with the upper bidiagonal matrix B; the
                     39: *            elements below the diagonal, with the array TAUQ, represent
                     40: *            the orthogonal matrix Q as a product of elementary
                     41: *            reflectors, and the elements above the first superdiagonal,
                     42: *            with the array TAUP, represent the orthogonal matrix P as
                     43: *            a product of elementary reflectors;
                     44: *          if m < n, the diagonal and the first subdiagonal are
                     45: *            overwritten with the lower bidiagonal matrix B; the
                     46: *            elements below the first subdiagonal, with the array TAUQ,
                     47: *            represent the orthogonal matrix Q as a product of
                     48: *            elementary reflectors, and the elements above the diagonal,
                     49: *            with the array TAUP, represent the orthogonal matrix P as
                     50: *            a product of elementary reflectors.
                     51: *          See Further Details.
                     52: *
                     53: *  LDA     (input) INTEGER
                     54: *          The leading dimension of the array A.  LDA >= max(1,M).
                     55: *
                     56: *  D       (output) DOUBLE PRECISION array, dimension (min(M,N))
                     57: *          The diagonal elements of the bidiagonal matrix B:
                     58: *          D(i) = A(i,i).
                     59: *
                     60: *  E       (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
                     61: *          The off-diagonal elements of the bidiagonal matrix B:
                     62: *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
                     63: *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
                     64: *
                     65: *  TAUQ    (output) DOUBLE PRECISION array dimension (min(M,N))
                     66: *          The scalar factors of the elementary reflectors which
                     67: *          represent the orthogonal matrix Q. See Further Details.
                     68: *
                     69: *  TAUP    (output) DOUBLE PRECISION array, dimension (min(M,N))
                     70: *          The scalar factors of the elementary reflectors which
                     71: *          represent the orthogonal matrix P. See Further Details.
                     72: *
                     73: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                     74: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                     75: *
                     76: *  LWORK   (input) INTEGER
                     77: *          The length of the array WORK.  LWORK >= max(1,M,N).
                     78: *          For optimum performance LWORK >= (M+N)*NB, where NB
                     79: *          is the optimal blocksize.
                     80: *
                     81: *          If LWORK = -1, then a workspace query is assumed; the routine
                     82: *          only calculates the optimal size of the WORK array, returns
                     83: *          this value as the first entry of the WORK array, and no error
                     84: *          message related to LWORK is issued by XERBLA.
                     85: *
                     86: *  INFO    (output) INTEGER
                     87: *          = 0:  successful exit
                     88: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
                     89: *
                     90: *  Further Details
                     91: *  ===============
                     92: *
                     93: *  The matrices Q and P are represented as products of elementary
                     94: *  reflectors:
                     95: *
                     96: *  If m >= n,
                     97: *
                     98: *     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
                     99: *
                    100: *  Each H(i) and G(i) has the form:
                    101: *
                    102: *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
                    103: *
                    104: *  where tauq and taup are real scalars, and v and u are real vectors;
                    105: *  v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
                    106: *  u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
                    107: *  tauq is stored in TAUQ(i) and taup in TAUP(i).
                    108: *
                    109: *  If m < n,
                    110: *
                    111: *     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
                    112: *
                    113: *  Each H(i) and G(i) has the form:
                    114: *
                    115: *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
                    116: *
                    117: *  where tauq and taup are real scalars, and v and u are real vectors;
                    118: *  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
                    119: *  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
                    120: *  tauq is stored in TAUQ(i) and taup in TAUP(i).
                    121: *
                    122: *  The contents of A on exit are illustrated by the following examples:
                    123: *
                    124: *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
                    125: *
                    126: *    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
                    127: *    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
                    128: *    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
                    129: *    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
                    130: *    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
                    131: *    (  v1  v2  v3  v4  v5 )
                    132: *
                    133: *  where d and e denote diagonal and off-diagonal elements of B, vi
                    134: *  denotes an element of the vector defining H(i), and ui an element of
                    135: *  the vector defining G(i).
                    136: *
                    137: *  =====================================================================
                    138: *
                    139: *     .. Parameters ..
                    140:       DOUBLE PRECISION   ONE
                    141:       PARAMETER          ( ONE = 1.0D+0 )
                    142: *     ..
                    143: *     .. Local Scalars ..
                    144:       LOGICAL            LQUERY
                    145:       INTEGER            I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
                    146:      $                   NBMIN, NX
                    147:       DOUBLE PRECISION   WS
                    148: *     ..
                    149: *     .. External Subroutines ..
                    150:       EXTERNAL           DGEBD2, DGEMM, DLABRD, XERBLA
                    151: *     ..
                    152: *     .. Intrinsic Functions ..
                    153:       INTRINSIC          DBLE, MAX, MIN
                    154: *     ..
                    155: *     .. External Functions ..
                    156:       INTEGER            ILAENV
                    157:       EXTERNAL           ILAENV
                    158: *     ..
                    159: *     .. Executable Statements ..
                    160: *
                    161: *     Test the input parameters
                    162: *
                    163:       INFO = 0
                    164:       NB = MAX( 1, ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) )
                    165:       LWKOPT = ( M+N )*NB
                    166:       WORK( 1 ) = DBLE( LWKOPT )
                    167:       LQUERY = ( LWORK.EQ.-1 )
                    168:       IF( M.LT.0 ) THEN
                    169:          INFO = -1
                    170:       ELSE IF( N.LT.0 ) THEN
                    171:          INFO = -2
                    172:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
                    173:          INFO = -4
                    174:       ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
                    175:          INFO = -10
                    176:       END IF
                    177:       IF( INFO.LT.0 ) THEN
                    178:          CALL XERBLA( 'DGEBRD', -INFO )
                    179:          RETURN
                    180:       ELSE IF( LQUERY ) THEN
                    181:          RETURN
                    182:       END IF
                    183: *
                    184: *     Quick return if possible
                    185: *
                    186:       MINMN = MIN( M, N )
                    187:       IF( MINMN.EQ.0 ) THEN
                    188:          WORK( 1 ) = 1
                    189:          RETURN
                    190:       END IF
                    191: *
                    192:       WS = MAX( M, N )
                    193:       LDWRKX = M
                    194:       LDWRKY = N
                    195: *
                    196:       IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
                    197: *
                    198: *        Set the crossover point NX.
                    199: *
                    200:          NX = MAX( NB, ILAENV( 3, 'DGEBRD', ' ', M, N, -1, -1 ) )
                    201: *
                    202: *        Determine when to switch from blocked to unblocked code.
                    203: *
                    204:          IF( NX.LT.MINMN ) THEN
                    205:             WS = ( M+N )*NB
                    206:             IF( LWORK.LT.WS ) THEN
                    207: *
                    208: *              Not enough work space for the optimal NB, consider using
                    209: *              a smaller block size.
                    210: *
                    211:                NBMIN = ILAENV( 2, 'DGEBRD', ' ', M, N, -1, -1 )
                    212:                IF( LWORK.GE.( M+N )*NBMIN ) THEN
                    213:                   NB = LWORK / ( M+N )
                    214:                ELSE
                    215:                   NB = 1
                    216:                   NX = MINMN
                    217:                END IF
                    218:             END IF
                    219:          END IF
                    220:       ELSE
                    221:          NX = MINMN
                    222:       END IF
                    223: *
                    224:       DO 30 I = 1, MINMN - NX, NB
                    225: *
                    226: *        Reduce rows and columns i:i+nb-1 to bidiagonal form and return
                    227: *        the matrices X and Y which are needed to update the unreduced
                    228: *        part of the matrix
                    229: *
                    230:          CALL DLABRD( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
                    231:      $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
                    232:      $                WORK( LDWRKX*NB+1 ), LDWRKY )
                    233: *
                    234: *        Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
                    235: *        of the form  A := A - V*Y' - X*U'
                    236: *
                    237:          CALL DGEMM( 'No transpose', 'Transpose', M-I-NB+1, N-I-NB+1,
                    238:      $               NB, -ONE, A( I+NB, I ), LDA,
                    239:      $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
                    240:      $               A( I+NB, I+NB ), LDA )
                    241:          CALL DGEMM( 'No transpose', 'No transpose', M-I-NB+1, N-I-NB+1,
                    242:      $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
                    243:      $               ONE, A( I+NB, I+NB ), LDA )
                    244: *
                    245: *        Copy diagonal and off-diagonal elements of B back into A
                    246: *
                    247:          IF( M.GE.N ) THEN
                    248:             DO 10 J = I, I + NB - 1
                    249:                A( J, J ) = D( J )
                    250:                A( J, J+1 ) = E( J )
                    251:    10       CONTINUE
                    252:          ELSE
                    253:             DO 20 J = I, I + NB - 1
                    254:                A( J, J ) = D( J )
                    255:                A( J+1, J ) = E( J )
                    256:    20       CONTINUE
                    257:          END IF
                    258:    30 CONTINUE
                    259: *
                    260: *     Use unblocked code to reduce the remainder of the matrix
                    261: *
                    262:       CALL DGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
                    263:      $             TAUQ( I ), TAUP( I ), WORK, IINFO )
                    264:       WORK( 1 ) = WS
                    265:       RETURN
                    266: *
                    267: *     End of DGEBRD
                    268: *
                    269:       END

CVSweb interface <joel.bertrand@systella.fr>