SRC/zlarzt.f

   1 *> \brief \b ZLARZT forms the triangular factor T of a block reflector H = I - vtvH.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLARZT + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarzt.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarzt.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarzt.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLARZT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          DIRECT, STOREV
  25 *       INTEGER            K, LDT, LDV, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       COMPLEX*16         T( LDT, * ), TAU( * ), V( LDV, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZLARZT forms the triangular factor T of a complex block reflector
  38 *> H of order > n, which is defined as a product of k elementary
  39 *> reflectors.
  40 *>
  41 *> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
  42 *>
  43 *> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
  44 *>
  45 *> If STOREV = 'C', the vector which defines the elementary reflector
  46 *> H(i) is stored in the i-th column of the array V, and
  47 *>
  48 *>    H  =  I - V * T * V**H
  49 *>
  50 *> If STOREV = 'R', the vector which defines the elementary reflector
  51 *> H(i) is stored in the i-th row of the array V, and
  52 *>
  53 *>    H  =  I - V**H * T * V
  54 *>
  55 *> Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
  56 *> \endverbatim
  57 *
  58 *  Arguments:
  59 *  ==========
  60 *
  61 *> \param[in] DIRECT
  62 *> \verbatim
  63 *>          DIRECT is CHARACTER*1
  64 *>          Specifies the order in which the elementary reflectors are
  65 *>          multiplied to form the block reflector:
  66 *>          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
  67 *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
  68 *> \endverbatim
  69 *>
  70 *> \param[in] STOREV
  71 *> \verbatim
  72 *>          STOREV is CHARACTER*1
  73 *>          Specifies how the vectors which define the elementary
  74 *>          reflectors are stored (see also Further Details):
  75 *>          = 'C': columnwise                        (not supported yet)
  76 *>          = 'R': rowwise
  77 *> \endverbatim
  78 *>
  79 *> \param[in] N
  80 *> \verbatim
  81 *>          N is INTEGER
  82 *>          The order of the block reflector H. N >= 0.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] K
  86 *> \verbatim
  87 *>          K is INTEGER
  88 *>          The order of the triangular factor T (= the number of
  89 *>          elementary reflectors). K >= 1.
  90 *> \endverbatim
  91 *>
  92 *> \param[in,out] V
  93 *> \verbatim
  94 *>          V is COMPLEX*16 array, dimension
  95 *>                               (LDV,K) if STOREV = 'C'
  96 *>                               (LDV,N) if STOREV = 'R'
  97 *>          The matrix V. See further details.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDV
 101 *> \verbatim
 102 *>          LDV is INTEGER
 103 *>          The leading dimension of the array V.
 104 *>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
 105 *> \endverbatim
 106 *>
 107 *> \param[in] TAU
 108 *> \verbatim
 109 *>          TAU is COMPLEX*16 array, dimension (K)
 110 *>          TAU(i) must contain the scalar factor of the elementary
 111 *>          reflector H(i).
 112 *> \endverbatim
 113 *>
 114 *> \param[out] T
 115 *> \verbatim
 116 *>          T is COMPLEX*16 array, dimension (LDT,K)
 117 *>          The k by k triangular factor T of the block reflector.
 118 *>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
 119 *>          lower triangular. The rest of the array is not used.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] LDT
 123 *> \verbatim
 124 *>          LDT is INTEGER
 125 *>          The leading dimension of the array T. LDT >= K.
 126 *> \endverbatim
 127 *
 128 *  Authors:
 129 *  ========
 130 *
 131 *> \author Univ. of Tennessee
 132 *> \author Univ. of California Berkeley
 133 *> \author Univ. of Colorado Denver
 134 *> \author NAG Ltd.
 135 *
 136 *> \date September 2012
 137 *
 138 *> \ingroup complex16OTHERcomputational
 139 *
 140 *> \par Contributors:
 141 *  ==================
 142 *>
 143 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 144 *
 145 *> \par Further Details:
 146 *  =====================
 147 *>
 148 *> \verbatim
 149 *>
 150 *>  The shape of the matrix V and the storage of the vectors which define
 151 *>  the H(i) is best illustrated by the following example with n = 5 and
 152 *>  k = 3. The elements equal to 1 are not stored; the corresponding
 153 *>  array elements are modified but restored on exit. The rest of the
 154 *>  array is not used.
 155 *>
 156 *>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
 157 *>
 158 *>                                              ______V_____
 159 *>         ( v1 v2 v3 )                        /            \
 160 *>         ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 )
 161 *>     V = ( v1 v2 v3 )                      ( v2 v2 v2 v2 v2 . . . 1   )
 162 *>         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     )
 163 *>         ( v1 v2 v3 )
 164 *>            .  .  .
 165 *>            .  .  .
 166 *>            1  .  .
 167 *>               1  .
 168 *>                  1
 169 *>
 170 *>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
 171 *>
 172 *>                                                        ______V_____
 173 *>            1                                          /            \
 174 *>            .  1                           ( 1 . . . . v1 v1 v1 v1 v1 )
 175 *>            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 )
 176 *>            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 )
 177 *>            .  .  .
 178 *>         ( v1 v2 v3 )
 179 *>         ( v1 v2 v3 )
 180 *>     V = ( v1 v2 v3 )
 181 *>         ( v1 v2 v3 )
 182 *>         ( v1 v2 v3 )
 183 *> \endverbatim
 184 *>
 185 *  =====================================================================
 186       SUBROUTINE ZLARZT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
 187 *
 188 *  -- LAPACK computational routine (version 3.4.2) --
 189 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 190 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 191 *     September 2012
 192 *
 193 *     .. Scalar Arguments ..
 194       CHARACTER          DIRECT, STOREV
 195       INTEGER            K, LDT, LDV, N
 196 *     ..
 197 *     .. Array Arguments ..
 198       COMPLEX*16         T( LDT, * ), TAU( * ), V( LDV, * )
 199 *     ..
 200 *
 201 *  =====================================================================
 202 *
 203 *     .. Parameters ..
 204       COMPLEX*16         ZERO
 205       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 206 *     ..
 207 *     .. Local Scalars ..
 208       INTEGER            I, INFO, J
 209 *     ..
 210 *     .. External Subroutines ..
 211       EXTERNAL           XERBLA, ZGEMV, ZLACGV, ZTRMV
 212 *     ..
 213 *     .. External Functions ..
 214       LOGICAL            LSAME
 215       EXTERNAL           LSAME
 216 *     ..
 217 *     .. Executable Statements ..
 218 *
 219 *     Check for currently supported options
 220 *
 221       INFO = 0
 222       IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
 223          INFO = -1
 224       ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
 225          INFO = -2
 226       END IF
 227       IF( INFO.NE.0 ) THEN
 228          CALL XERBLA( 'ZLARZT', -INFO )
 229          RETURN
 230       END IF
 231 *
 232       DO 20 I = K, 1, -1
 233          IF( TAU( I ).EQ.ZERO ) THEN
 234 *
 235 *           H(i)  =  I
 236 *
 237             DO 10 J = I, K
 238                T( J, I ) = ZERO
 239    10       CONTINUE
 240          ELSE
 241 *
 242 *           general case
 243 *
 244             IF( I.LT.K ) THEN
 245 *
 246 *              T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)**H
 247 *
 248                CALL ZLACGV( N, V( I, 1 ), LDV )
 249                CALL ZGEMV( 'No transpose', K-I, N, -TAU( I ),
 250      $                     V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO,
 251      $                     T( I+1, I ), 1 )
 252                CALL ZLACGV( N, V( I, 1 ), LDV )
 253 *
 254 *              T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i)
 255 *
 256                CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
 257      $                     T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
 258             END IF
 259             T( I, I ) = TAU( I )
 260          END IF
 261    20 CONTINUE
 262       RETURN
 263 *
 264 *     End of ZLARZT
 265 *
 266       END