Return to dgebrd.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence 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 2011
  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.4.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 2011
  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
  231:       DOUBLE PRECISION   WS
  232: *     ..
  233: *     .. External Subroutines ..
  234:       EXTERNAL           DGEBD2, DGEMM, DLABRD, XERBLA
  235: *     ..
  236: *     .. Intrinsic Functions ..
  237:       INTRINSIC          DBLE, MAX, MIN
  238: *     ..
  239: *     .. External Functions ..
  240:       INTEGER            ILAENV
  241:       EXTERNAL           ILAENV
  242: *     ..
  243: *     .. Executable Statements ..
  244: *
  245: *     Test the input parameters
  246: *
  247:       INFO = 0
  248:       NB = MAX( 1, ILAENV( 1, 'DGEBRD', ' ', M, N, -1, -1 ) )
  249:       LWKOPT = ( M+N )*NB
  250:       WORK( 1 ) = DBLE( LWKOPT )
  251:       LQUERY = ( LWORK.EQ.-1 )
  252:       IF( M.LT.0 ) THEN
  253:          INFO = -1
  254:       ELSE IF( N.LT.0 ) THEN
  255:          INFO = -2
  256:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  257:          INFO = -4
  258:       ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
  259:          INFO = -10
  260:       END IF
  261:       IF( INFO.LT.0 ) THEN
  262:          CALL XERBLA( 'DGEBRD', -INFO )
  263:          RETURN
  264:       ELSE IF( LQUERY ) THEN
  265:          RETURN
  266:       END IF
  267: *
  268: *     Quick return if possible
  269: *
  270:       MINMN = MIN( M, N )
  271:       IF( MINMN.EQ.0 ) THEN
  272:          WORK( 1 ) = 1
  273:          RETURN
  274:       END IF
  275: *
  276:       WS = MAX( M, N )
  277:       LDWRKX = M
  278:       LDWRKY = N
  279: *
  280:       IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
  281: *
  282: *        Set the crossover point NX.
  283: *
  284:          NX = MAX( NB, ILAENV( 3, 'DGEBRD', ' ', M, N, -1, -1 ) )
  285: *
  286: *        Determine when to switch from blocked to unblocked code.
  287: *
  288:          IF( NX.LT.MINMN ) THEN
  289:             WS = ( M+N )*NB
  290:             IF( LWORK.LT.WS ) THEN
  291: *
  292: *              Not enough work space for the optimal NB, consider using
  293: *              a smaller block size.
  294: *
  295:                NBMIN = ILAENV( 2, 'DGEBRD', ' ', M, N, -1, -1 )
  296:                IF( LWORK.GE.( M+N )*NBMIN ) THEN
  297:                   NB = LWORK / ( M+N )
  298:                ELSE
  299:                   NB = 1
  300:                   NX = MINMN
  301:                END IF
  302:             END IF
  303:          END IF
  304:       ELSE
  305:          NX = MINMN
  306:       END IF
  307: *
  308:       DO 30 I = 1, MINMN - NX, NB
  309: *
  310: *        Reduce rows and columns i:i+nb-1 to bidiagonal form and return
  311: *        the matrices X and Y which are needed to update the unreduced
  312: *        part of the matrix
  313: *
  314:          CALL DLABRD( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
  315:      $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
  316:      $                WORK( LDWRKX*NB+1 ), LDWRKY )
  317: *
  318: *        Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
  319: *        of the form  A := A - V*Y**T - X*U**T
  320: *
  321:          CALL DGEMM( 'No transpose', 'Transpose', M-I-NB+1, N-I-NB+1,
  322:      $               NB, -ONE, A( I+NB, I ), LDA,
  323:      $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
  324:      $               A( I+NB, I+NB ), LDA )
  325:          CALL DGEMM( 'No transpose', 'No transpose', M-I-NB+1, N-I-NB+1,
  326:      $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
  327:      $               ONE, A( I+NB, I+NB ), LDA )
  328: *
  329: *        Copy diagonal and off-diagonal elements of B back into A
  330: *
  331:          IF( M.GE.N ) THEN
  332:             DO 10 J = I, I + NB - 1
  333:                A( J, J ) = D( J )
  334:                A( J, J+1 ) = E( J )
  335:    10       CONTINUE
  336:          ELSE
  337:             DO 20 J = I, I + NB - 1
  338:                A( J, J ) = D( J )
  339:                A( J+1, J ) = E( J )
  340:    20       CONTINUE
  341:          END IF
  342:    30 CONTINUE
  343: *
  344: *     Use unblocked code to reduce the remainder of the matrix
  345: *
  346:       CALL DGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
  347:      $             TAUQ( I ), TAUP( I ), WORK, IINFO )
  348:       WORK( 1 ) = WS
  349:       RETURN
  350: *
  351: *     End of DGEBRD
  352: *
  353:       END