Return to dgebrd.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief \b DGEBRD
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DGEBRD + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgebrd.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgebrd.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgebrd.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
   22: *                          INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       INTEGER            INFO, LDA, LWORK, M, N
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAUP( * ),
   29: *      $                   TAUQ( * ), WORK( * )
   30: *       ..
   31: *
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> DGEBRD reduces a general real M-by-N matrix A to upper or lower
   39: *> bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
   40: *>
   41: *> If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
   42: *> \endverbatim
   43: *
   44: *  Arguments:
   45: *  ==========
   46: *
   47: *> \param[in] M
   48: *> \verbatim
   49: *>          M is INTEGER
   50: *>          The number of rows in the matrix A.  M >= 0.
   51: *> \endverbatim
   52: *>
   53: *> \param[in] N
   54: *> \verbatim
   55: *>          N is INTEGER
   56: *>          The number of columns in the matrix A.  N >= 0.
   57: *> \endverbatim
   58: *>
   59: *> \param[in,out] A
   60: *> \verbatim
   61: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   62: *>          On entry, the M-by-N general matrix to be reduced.
   63: *>          On exit,
   64: *>          if m >= n, the diagonal and the first superdiagonal are
   65: *>            overwritten with the upper bidiagonal matrix B; the
   66: *>            elements below the diagonal, with the array TAUQ, represent
   67: *>            the orthogonal matrix Q as a product of elementary
   68: *>            reflectors, and the elements above the first superdiagonal,
   69: *>            with the array TAUP, represent the orthogonal matrix P as
   70: *>            a product of elementary reflectors;
   71: *>          if m < n, the diagonal and the first subdiagonal are
   72: *>            overwritten with the lower bidiagonal matrix B; the
   73: *>            elements below the first subdiagonal, with the array TAUQ,
   74: *>            represent the orthogonal matrix Q as a product of
   75: *>            elementary reflectors, and the elements above the diagonal,
   76: *>            with the array TAUP, represent the orthogonal matrix P as
   77: *>            a product of elementary reflectors.
   78: *>          See Further Details.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] LDA
   82: *> \verbatim
   83: *>          LDA is INTEGER
   84: *>          The leading dimension of the array A.  LDA >= max(1,M).
   85: *> \endverbatim
   86: *>
   87: *> \param[out] D
   88: *> \verbatim
   89: *>          D is DOUBLE PRECISION array, dimension (min(M,N))
   90: *>          The diagonal elements of the bidiagonal matrix B:
   91: *>          D(i) = A(i,i).
   92: *> \endverbatim
   93: *>
   94: *> \param[out] E
   95: *> \verbatim
   96: *>          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
   97: *>          The off-diagonal elements of the bidiagonal matrix B:
   98: *>          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
   99: *>          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
  100: *> \endverbatim
  101: *>
  102: *> \param[out] TAUQ
  103: *> \verbatim
  104: *>          TAUQ is DOUBLE PRECISION array, dimension (min(M,N))
  105: *>          The scalar factors of the elementary reflectors which
  106: *>          represent the orthogonal matrix Q. See Further Details.
  107: *> \endverbatim
  108: *>
  109: *> \param[out] TAUP
  110: *> \verbatim
  111: *>          TAUP is DOUBLE PRECISION array, dimension (min(M,N))
  112: *>          The scalar factors of the elementary reflectors which
  113: *>          represent the orthogonal matrix P. See Further Details.
  114: *> \endverbatim
  115: *>
  116: *> \param[out] WORK
  117: *> \verbatim
  118: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  119: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  120: *> \endverbatim
  121: *>
  122: *> \param[in] LWORK
  123: *> \verbatim
  124: *>          LWORK is INTEGER
  125: *>          The length of the array WORK.  LWORK >= max(1,M,N).
  126: *>          For optimum performance LWORK >= (M+N)*NB, where NB
  127: *>          is the optimal blocksize.
  128: *>
  129: *>          If LWORK = -1, then a workspace query is assumed; the routine
  130: *>          only calculates the optimal size of the WORK array, returns
  131: *>          this value as the first entry of the WORK array, and no error
  132: *>          message related to LWORK is issued by XERBLA.
  133: *> \endverbatim
  134: *>
  135: *> \param[out] INFO
  136: *> \verbatim
  137: *>          INFO is INTEGER
  138: *>          = 0:  successful exit
  139: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  140: *> \endverbatim
  141: *
  142: *  Authors:
  143: *  ========
  144: *
  145: *> \author Univ. of Tennessee
  146: *> \author Univ. of California Berkeley
  147: *> \author Univ. of Colorado Denver
  148: *> \author NAG Ltd.
  149: *
  150: *> \date November 2017
  151: *
  152: *> \ingroup doubleGEcomputational
  153: *
  154: *> \par Further Details:
  155: *  =====================
  156: *>
  157: *> \verbatim
  158: *>
  159: *>  The matrices Q and P are represented as products of elementary
  160: *>  reflectors:
  161: *>
  162: *>  If m >= n,
  163: *>
  164: *>     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
  165: *>
  166: *>  Each H(i) and G(i) has the form:
  167: *>
  168: *>     H(i) = I - tauq * v * v**T  and G(i) = I - taup * u * u**T
  169: *>
  170: *>  where tauq and taup are real scalars, and v and u are real vectors;
  171: *>  v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
  172: *>  u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
  173: *>  tauq is stored in TAUQ(i) and taup in TAUP(i).
  174: *>
  175: *>  If m < n,
  176: *>
  177: *>     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
  178: *>
  179: *>  Each H(i) and G(i) has the form:
  180: *>
  181: *>     H(i) = I - tauq * v * v**T  and G(i) = I - taup * u * u**T
  182: *>
  183: *>  where tauq and taup are real scalars, and v and u are real vectors;
  184: *>  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
  185: *>  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
  186: *>  tauq is stored in TAUQ(i) and taup in TAUP(i).
  187: *>
  188: *>  The contents of A on exit are illustrated by the following examples:
  189: *>
  190: *>  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
  191: *>
  192: *>    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
  193: *>    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
  194: *>    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
  195: *>    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
  196: *>    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
  197: *>    (  v1  v2  v3  v4  v5 )
  198: *>
  199: *>  where d and e denote diagonal and off-diagonal elements of B, vi
  200: *>  denotes an element of the vector defining H(i), and ui an element of
  201: *>  the vector defining G(i).
  202: *> \endverbatim
  203: *>
  204: *  =====================================================================
  205:       SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
  206:      $                   INFO )
  207: *
  208: *  -- LAPACK computational routine (version 3.8.0) --
  209: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  210: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  211: *     November 2017
  212: *
  213: *     .. Scalar Arguments ..
  214:       INTEGER            INFO, LDA, LWORK, M, N
  215: *     ..
  216: *     .. Array Arguments ..
  217:       DOUBLE PRECISION   A( LDA, * ), D( * ), E( * ), TAUP( * ),
  218:      $                   TAUQ( * ), WORK( * )
  219: *     ..
  220: *
  221: *  =====================================================================
  222: *
  223: *     .. Parameters ..
  224:       DOUBLE PRECISION   ONE
  225:       PARAMETER          ( ONE = 1.0D+0 )
  226: *     ..
  227: *     .. Local Scalars ..
  228:       LOGICAL            LQUERY
  229:       INTEGER            I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
  230:      $                   NBMIN, NX, WS
  231: *     ..
  232: *     .. External Subroutines ..
  233:       EXTERNAL           DGEBD2, DGEMM, DLABRD, XERBLA
  234: *     ..
  235: *     .. Intrinsic Functions ..
  236:       INTRINSIC          DBLE, MAX, MIN
  237: *     ..
  238: *     .. External Functions ..
  239:       INTEGER            ILAENV
  240:       EXTERNAL           ILAENV
  241: *     ..
  242: *     .. Executable Statements ..
  243: *
  244: *     Test the input parameters
  245: *
  246:       INFO = 0
  247:       NB = MAX( 1, ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) )
  248:       LWKOPT = ( M+N )*NB
  249:       WORK( 1 ) = DBLE( LWKOPT )
  250:       LQUERY = ( LWORK.EQ.-1 )
  251:       IF( M.LT.0 ) THEN
  252:          INFO = -1
  253:       ELSE IF( N.LT.0 ) THEN
  254:          INFO = -2
  255:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  256:          INFO = -4
  257:       ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
  258:          INFO = -10
  259:       END IF
  260:       IF( INFO.LT.0 ) THEN
  261:          CALL XERBLA( 'DGEBRD', -INFO )
  262:          RETURN
  263:       ELSE IF( LQUERY ) THEN
  264:          RETURN
  265:       END IF
  266: *
  267: *     Quick return if possible
  268: *
  269:       MINMN = MIN( M, N )
  270:       IF( MINMN.EQ.0 ) THEN
  271:          WORK( 1 ) = 1
  272:          RETURN
  273:       END IF
  274: *
  275:       WS = MAX( M, N )
  276:       LDWRKX = M
  277:       LDWRKY = N
  278: *
  279:       IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
  280: *
  281: *        Set the crossover point NX.
  282: *
  283:          NX = MAX( NB, ILAENV( 3, 'DGEBRD', ' ', M, N, -1, -1 ) )
  284: *
  285: *        Determine when to switch from blocked to unblocked code.
  286: *
  287:          IF( NX.LT.MINMN ) THEN
  288:             WS = ( M+N )*NB
  289:             IF( LWORK.LT.WS ) THEN
  290: *
  291: *              Not enough work space for the optimal NB, consider using
  292: *              a smaller block size.
  293: *
  294:                NBMIN = ILAENV( 2, 'DGEBRD', ' ', M, N, -1, -1 )
  295:                IF( LWORK.GE.( M+N )*NBMIN ) THEN
  296:                   NB = LWORK / ( M+N )
  297:                ELSE
  298:                   NB = 1
  299:                   NX = MINMN
  300:                END IF
  301:             END IF
  302:          END IF
  303:       ELSE
  304:          NX = MINMN
  305:       END IF
  306: *
  307:       DO 30 I = 1, MINMN - NX, NB
  308: *
  309: *        Reduce rows and columns i:i+nb-1 to bidiagonal form and return
  310: *        the matrices X and Y which are needed to update the unreduced
  311: *        part of the matrix
  312: *
  313:          CALL DLABRD( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
  314:      $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
  315:      $                WORK( LDWRKX*NB+1 ), LDWRKY )
  316: *
  317: *        Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
  318: *        of the form  A := A - V*Y**T - X*U**T
  319: *
  320:          CALL DGEMM( 'No transpose', 'Transpose', M-I-NB+1, N-I-NB+1,
  321:      $               NB, -ONE, A( I+NB, I ), LDA,
  322:      $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
  323:      $               A( I+NB, I+NB ), LDA )
  324:          CALL DGEMM( 'No transpose', 'No transpose', M-I-NB+1, N-I-NB+1,
  325:      $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
  326:      $               ONE, A( I+NB, I+NB ), LDA )
  327: *
  328: *        Copy diagonal and off-diagonal elements of B back into A
  329: *
  330:          IF( M.GE.N ) THEN
  331:             DO 10 J = I, I + NB - 1
  332:                A( J, J ) = D( J )
  333:                A( J, J+1 ) = E( J )
  334:    10       CONTINUE
  335:          ELSE
  336:             DO 20 J = I, I + NB - 1
  337:                A( J, J ) = D( J )
  338:                A( J+1, J ) = E( J )
  339:    20       CONTINUE
  340:          END IF
  341:    30 CONTINUE
  342: *
  343: *     Use unblocked code to reduce the remainder of the matrix
  344: *
  345:       CALL DGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
  346:      $             TAUQ( I ), TAUP( I ), WORK, IINFO )
  347:       WORK( 1 ) = WS
  348:       RETURN
  349: *
  350: *     End of DGEBRD
  351: *
  352:       END