SRC/cgebrd.f

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