SRC/zlahr2.f

   1 *> \brief \b ZLAHR2 reduces the specified number of first columns of a general rectangular matrix A so that elements below the specified subdiagonal are zero, and returns auxiliary matrices which are needed to apply the transformation to the unreduced part of A.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLAHR2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahr2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahr2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahr2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            K, LDA, LDT, LDY, N, NB
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       COMPLEX*16        A( LDA, * ), T( LDT, NB ), TAU( NB ),
  28 *      $                   Y( LDY, NB )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1)
  38 *> matrix A so that elements below the k-th subdiagonal are zero. The
  39 *> reduction is performed by an unitary similarity transformation
  40 *> Q**H * A * Q. The routine returns the matrices V and T which determine
  41 *> Q as a block reflector I - V*T*V**H, and also the matrix Y = A * V * T.
  42 *>
  43 *> This is an auxiliary routine called by ZGEHRD.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] N
  50 *> \verbatim
  51 *>          N is INTEGER
  52 *>          The order of the matrix A.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] K
  56 *> \verbatim
  57 *>          K is INTEGER
  58 *>          The offset for the reduction. Elements below the k-th
  59 *>          subdiagonal in the first NB columns are reduced to zero.
  60 *>          K < N.
  61 *> \endverbatim
  62 *>
  63 *> \param[in] NB
  64 *> \verbatim
  65 *>          NB is INTEGER
  66 *>          The number of columns to be reduced.
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] A
  70 *> \verbatim
  71 *>          A is COMPLEX*16 array, dimension (LDA,N-K+1)
  72 *>          On entry, the n-by-(n-k+1) general matrix A.
  73 *>          On exit, the elements on and above the k-th subdiagonal in
  74 *>          the first NB columns are overwritten with the corresponding
  75 *>          elements of the reduced matrix; the elements below the k-th
  76 *>          subdiagonal, with the array TAU, represent the matrix Q as a
  77 *>          product of elementary reflectors. The other columns of A are
  78 *>          unchanged. 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,N).
  85 *> \endverbatim
  86 *>
  87 *> \param[out] TAU
  88 *> \verbatim
  89 *>          TAU is COMPLEX*16 array, dimension (NB)
  90 *>          The scalar factors of the elementary reflectors. See Further
  91 *>          Details.
  92 *> \endverbatim
  93 *>
  94 *> \param[out] T
  95 *> \verbatim
  96 *>          T is COMPLEX*16 array, dimension (LDT,NB)
  97 *>          The upper triangular matrix T.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDT
 101 *> \verbatim
 102 *>          LDT is INTEGER
 103 *>          The leading dimension of the array T.  LDT >= NB.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] Y
 107 *> \verbatim
 108 *>          Y is COMPLEX*16 array, dimension (LDY,NB)
 109 *>          The n-by-nb matrix Y.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] LDY
 113 *> \verbatim
 114 *>          LDY is INTEGER
 115 *>          The leading dimension of the array Y. LDY >= N.
 116 *> \endverbatim
 117 *
 118 *  Authors:
 119 *  ========
 120 *
 121 *> \author Univ. of Tennessee
 122 *> \author Univ. of California Berkeley
 123 *> \author Univ. of Colorado Denver
 124 *> \author NAG Ltd.
 125 *
 126 *> \date September 2012
 127 *
 128 *> \ingroup complex16OTHERauxiliary
 129 *
 130 *> \par Further Details:
 131 *  =====================
 132 *>
 133 *> \verbatim
 134 *>
 135 *>  The matrix Q is represented as a product of nb elementary reflectors
 136 *>
 137 *>     Q = H(1) H(2) . . . H(nb).
 138 *>
 139 *>  Each H(i) has the form
 140 *>
 141 *>     H(i) = I - tau * v * v**H
 142 *>
 143 *>  where tau is a complex scalar, and v is a complex vector with
 144 *>  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
 145 *>  A(i+k+1:n,i), and tau in TAU(i).
 146 *>
 147 *>  The elements of the vectors v together form the (n-k+1)-by-nb matrix
 148 *>  V which is needed, with T and Y, to apply the transformation to the
 149 *>  unreduced part of the matrix, using an update of the form:
 150 *>  A := (I - V*T*V**H) * (A - Y*V**H).
 151 *>
 152 *>  The contents of A on exit are illustrated by the following example
 153 *>  with n = 7, k = 3 and nb = 2:
 154 *>
 155 *>     ( a   a   a   a   a )
 156 *>     ( a   a   a   a   a )
 157 *>     ( a   a   a   a   a )
 158 *>     ( h   h   a   a   a )
 159 *>     ( v1  h   a   a   a )
 160 *>     ( v1  v2  a   a   a )
 161 *>     ( v1  v2  a   a   a )
 162 *>
 163 *>  where a denotes an element of the original matrix A, h denotes a
 164 *>  modified element of the upper Hessenberg matrix H, and vi denotes an
 165 *>  element of the vector defining H(i).
 166 *>
 167 *>  This subroutine is a slight modification of LAPACK-3.0's DLAHRD
 168 *>  incorporating improvements proposed by Quintana-Orti and Van de
 169 *>  Gejin. Note that the entries of A(1:K,2:NB) differ from those
 170 *>  returned by the original LAPACK-3.0's DLAHRD routine. (This
 171 *>  subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)
 172 *> \endverbatim
 173 *
 174 *> \par References:
 175 *  ================
 176 *>
 177 *>  Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
 178 *>  performance of reduction to Hessenberg form," ACM Transactions on
 179 *>  Mathematical Software, 32(2):180-194, June 2006.
 180 *>
 181 *  =====================================================================
 182       SUBROUTINE ZLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
 183 *
 184 *  -- LAPACK auxiliary routine (version 3.4.2) --
 185 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 186 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 187 *     September 2012
 188 *
 189 *     .. Scalar Arguments ..
 190       INTEGER            K, LDA, LDT, LDY, N, NB
 191 *     ..
 192 *     .. Array Arguments ..
 193       COMPLEX*16        A( LDA, * ), T( LDT, NB ), TAU( NB ),
 194      $                   Y( LDY, NB )
 195 *     ..
 196 *
 197 *  =====================================================================
 198 *
 199 *     .. Parameters ..
 200       COMPLEX*16        ZERO, ONE
 201       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
 202      $                     ONE = ( 1.0D+0, 0.0D+0 ) )
 203 *     ..
 204 *     .. Local Scalars ..
 205       INTEGER            I
 206       COMPLEX*16        EI
 207 *     ..
 208 *     .. External Subroutines ..
 209       EXTERNAL           ZAXPY, ZCOPY, ZGEMM, ZGEMV, ZLACPY,
 210      $                   ZLARFG, ZSCAL, ZTRMM, ZTRMV, ZLACGV
 211 *     ..
 212 *     .. Intrinsic Functions ..
 213       INTRINSIC          MIN
 214 *     ..
 215 *     .. Executable Statements ..
 216 *
 217 *     Quick return if possible
 218 *
 219       IF( N.LE.1 )
 220      $   RETURN
 221 *
 222       DO 10 I = 1, NB
 223          IF( I.GT.1 ) THEN
 224 *
 225 *           Update A(K+1:N,I)
 226 *
 227 *           Update I-th column of A - Y * V**H
 228 *
 229             CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
 230             CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY,
 231      $                  A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
 232             CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
 233 *
 234 *           Apply I - V * T**H * V**H to this column (call it b) from the
 235 *           left, using the last column of T as workspace
 236 *
 237 *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
 238 *                    ( V2 )             ( b2 )
 239 *
 240 *           where V1 is unit lower triangular
 241 *
 242 *           w := V1**H * b1
 243 *
 244             CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
 245             CALL ZTRMV( 'Lower', 'Conjugate transpose', 'UNIT',
 246      $                  I-1, A( K+1, 1 ),
 247      $                  LDA, T( 1, NB ), 1 )
 248 *
 249 *           w := w + V2**H * b2
 250 *
 251             CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1,
 252      $                  ONE, A( K+I, 1 ),
 253      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
 254 *
 255 *           w := T**H * w
 256 *
 257             CALL ZTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT',
 258      $                  I-1, T, LDT,
 259      $                  T( 1, NB ), 1 )
 260 *
 261 *           b2 := b2 - V2*w
 262 *
 263             CALL ZGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE,
 264      $                  A( K+I, 1 ),
 265      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
 266 *
 267 *           b1 := b1 - V1*w
 268 *
 269             CALL ZTRMV( 'Lower', 'NO TRANSPOSE',
 270      $                  'UNIT', I-1,
 271      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
 272             CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
 273 *
 274             A( K+I-1, I-1 ) = EI
 275          END IF
 276 *
 277 *        Generate the elementary reflector H(I) to annihilate
 278 *        A(K+I+1:N,I)
 279 *
 280          CALL ZLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
 281      $                TAU( I ) )
 282          EI = A( K+I, I )
 283          A( K+I, I ) = ONE
 284 *
 285 *        Compute  Y(K+1:N,I)
 286 *
 287          CALL ZGEMV( 'NO TRANSPOSE', N-K, N-K-I+1,
 288      $               ONE, A( K+1, I+1 ),
 289      $               LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 )
 290          CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1,
 291      $               ONE, A( K+I, 1 ), LDA,
 292      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
 293          CALL ZGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE,
 294      $               Y( K+1, 1 ), LDY,
 295      $               T( 1, I ), 1, ONE, Y( K+1, I ), 1 )
 296          CALL ZSCAL( N-K, TAU( I ), Y( K+1, I ), 1 )
 297 *
 298 *        Compute T(1:I,I)
 299 *
 300          CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
 301          CALL ZTRMV( 'Upper', 'No Transpose', 'NON-UNIT',
 302      $               I-1, T, LDT,
 303      $               T( 1, I ), 1 )
 304          T( I, I ) = TAU( I )
 305 *
 306    10 CONTINUE
 307       A( K+NB, NB ) = EI
 308 *
 309 *     Compute Y(1:K,1:NB)
 310 *
 311       CALL ZLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY )
 312       CALL ZTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE',
 313      $            'UNIT', K, NB,
 314      $            ONE, A( K+1, 1 ), LDA, Y, LDY )
 315       IF( N.GT.K+NB )
 316      $   CALL ZGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K,
 317      $               NB, N-K-NB, ONE,
 318      $               A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y,
 319      $               LDY )
 320       CALL ZTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE',
 321      $            'NON-UNIT', K, NB,
 322      $            ONE, T, LDT, Y, LDY )
 323 *
 324       RETURN
 325 *
 326 *     End of ZLAHR2
 327 *
 328       END