SRC/dgebrd.f

   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