SRC/clarft.f

   1 *> \brief \b CLARFT 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 CLARFT + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLARFT( 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            T( LDT, * ), TAU( * ), V( LDV, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> CLARFT forms the triangular factor T of a complex block reflector H
  38 *> of order n, which is defined as a product of k elementary reflectors.
  39 *>
  40 *> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
  41 *>
  42 *> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
  43 *>
  44 *> If STOREV = 'C', the vector which defines the elementary reflector
  45 *> H(i) is stored in the i-th column of the array V, and
  46 *>
  47 *>    H  =  I - V * T * V**H
  48 *>
  49 *> If STOREV = 'R', the vector which defines the elementary reflector
  50 *> H(i) is stored in the i-th row of the array V, and
  51 *>
  52 *>    H  =  I - V**H * T * V
  53 *> \endverbatim
  54 *
  55 *  Arguments:
  56 *  ==========
  57 *
  58 *> \param[in] DIRECT
  59 *> \verbatim
  60 *>          DIRECT is CHARACTER*1
  61 *>          Specifies the order in which the elementary reflectors are
  62 *>          multiplied to form the block reflector:
  63 *>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
  64 *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
  65 *> \endverbatim
  66 *>
  67 *> \param[in] STOREV
  68 *> \verbatim
  69 *>          STOREV is CHARACTER*1
  70 *>          Specifies how the vectors which define the elementary
  71 *>          reflectors are stored (see also Further Details):
  72 *>          = 'C': columnwise
  73 *>          = 'R': rowwise
  74 *> \endverbatim
  75 *>
  76 *> \param[in] N
  77 *> \verbatim
  78 *>          N is INTEGER
  79 *>          The order of the block reflector H. N >= 0.
  80 *> \endverbatim
  81 *>
  82 *> \param[in] K
  83 *> \verbatim
  84 *>          K is INTEGER
  85 *>          The order of the triangular factor T (= the number of
  86 *>          elementary reflectors). K >= 1.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] V
  90 *> \verbatim
  91 *>          V is COMPLEX array, dimension
  92 *>                               (LDV,K) if STOREV = 'C'
  93 *>                               (LDV,N) if STOREV = 'R'
  94 *>          The matrix V. See further details.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] LDV
  98 *> \verbatim
  99 *>          LDV is INTEGER
 100 *>          The leading dimension of the array V.
 101 *>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] TAU
 105 *> \verbatim
 106 *>          TAU is COMPLEX array, dimension (K)
 107 *>          TAU(i) must contain the scalar factor of the elementary
 108 *>          reflector H(i).
 109 *> \endverbatim
 110 *>
 111 *> \param[out] T
 112 *> \verbatim
 113 *>          T is COMPLEX array, dimension (LDT,K)
 114 *>          The k by k triangular factor T of the block reflector.
 115 *>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
 116 *>          lower triangular. The rest of the array is not used.
 117 *> \endverbatim
 118 *>
 119 *> \param[in] LDT
 120 *> \verbatim
 121 *>          LDT is INTEGER
 122 *>          The leading dimension of the array T. LDT >= K.
 123 *> \endverbatim
 124 *
 125 *  Authors:
 126 *  ========
 127 *
 128 *> \author Univ. of Tennessee
 129 *> \author Univ. of California Berkeley
 130 *> \author Univ. of Colorado Denver
 131 *> \author NAG Ltd.
 132 *
 133 *> \date September 2012
 134 *
 135 *> \ingroup complexOTHERauxiliary
 136 *
 137 *> \par Further Details:
 138 *  =====================
 139 *>
 140 *> \verbatim
 141 *>
 142 *>  The shape of the matrix V and the storage of the vectors which define
 143 *>  the H(i) is best illustrated by the following example with n = 5 and
 144 *>  k = 3. The elements equal to 1 are not stored.
 145 *>
 146 *>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
 147 *>
 148 *>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
 149 *>                   ( v1  1    )                     (     1 v2 v2 v2 )
 150 *>                   ( v1 v2  1 )                     (        1 v3 v3 )
 151 *>                   ( v1 v2 v3 )
 152 *>                   ( v1 v2 v3 )
 153 *>
 154 *>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
 155 *>
 156 *>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
 157 *>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
 158 *>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
 159 *>                   (     1 v3 )
 160 *>                   (        1 )
 161 *> \endverbatim
 162 *>
 163 *  =====================================================================
 164       SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
 165 *
 166 *  -- LAPACK auxiliary routine (version 3.4.2) --
 167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 169 *     September 2012
 170 *
 171 *     .. Scalar Arguments ..
 172       CHARACTER          DIRECT, STOREV
 173       INTEGER            K, LDT, LDV, N
 174 *     ..
 175 *     .. Array Arguments ..
 176       COMPLEX            T( LDT, * ), TAU( * ), V( LDV, * )
 177 *     ..
 178 *
 179 *  =====================================================================
 180 *
 181 *     .. Parameters ..
 182       COMPLEX            ONE, ZERO
 183       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
 184      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
 185 *     ..
 186 *     .. Local Scalars ..
 187       INTEGER            I, J, PREVLASTV, LASTV
 188 *     ..
 189 *     .. External Subroutines ..
 190       EXTERNAL           CGEMV, CTRMV
 191 *     ..
 192 *     .. External Functions ..
 193       LOGICAL            LSAME
 194       EXTERNAL           LSAME
 195 *     ..
 196 *     .. Executable Statements ..
 197 *
 198 *     Quick return if possible
 199 *
 200       IF( N.EQ.0 )
 201      $   RETURN
 202 *
 203       IF( LSAME( DIRECT, 'F' ) ) THEN
 204          PREVLASTV = N
 205          DO I = 1, K
 206             PREVLASTV = MAX( PREVLASTV, I )
 207             IF( TAU( I ).EQ.ZERO ) THEN
 208 *
 209 *              H(i)  =  I
 210 *
 211                DO J = 1, I
 212                   T( J, I ) = ZERO
 213                END DO
 214             ELSE
 215 *
 216 *              general case
 217 *
 218                IF( LSAME( STOREV, 'C' ) ) THEN
 219 *                 Skip any trailing zeros.
 220                   DO LASTV = N, I+1, -1
 221                      IF( V( LASTV, I ).NE.ZERO ) EXIT
 222                   END DO
 223                   DO J = 1, I-1
 224                      T( J, I ) = -TAU( I ) * CONJG( V( I , J ) )
 225                   END DO
 226                   J = MIN( LASTV, PREVLASTV )
 227 *
 228 *                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
 229 *
 230                   CALL CGEMV( 'Conjugate transpose', J-I, I-1,
 231      $                        -TAU( I ), V( I+1, 1 ), LDV,
 232      $                        V( I+1, I ), 1,
 233      $                        ONE, T( 1, I ), 1 )
 234                ELSE
 235 *                 Skip any trailing zeros.
 236                   DO LASTV = N, I+1, -1
 237                      IF( V( I, LASTV ).NE.ZERO ) EXIT
 238                   END DO
 239                   DO J = 1, I-1
 240                      T( J, I ) = -TAU( I ) * V( J , I )
 241                   END DO
 242                   J = MIN( LASTV, PREVLASTV )
 243 *
 244 *                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
 245 *
 246                   CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
 247      $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
 248      $                        ONE, T( 1, I ), LDT )
 249                END IF
 250 *
 251 *              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
 252 *
 253                CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
 254      $                     LDT, T( 1, I ), 1 )
 255                T( I, I ) = TAU( I )
 256                IF( I.GT.1 ) THEN
 257                   PREVLASTV = MAX( PREVLASTV, LASTV )
 258                ELSE
 259                   PREVLASTV = LASTV
 260                END IF
 261             END IF
 262          END DO
 263       ELSE
 264          PREVLASTV = 1
 265          DO I = K, 1, -1
 266             IF( TAU( I ).EQ.ZERO ) THEN
 267 *
 268 *              H(i)  =  I
 269 *
 270                DO J = I, K
 271                   T( J, I ) = ZERO
 272                END DO
 273             ELSE
 274 *
 275 *              general case
 276 *
 277                IF( I.LT.K ) THEN
 278                   IF( LSAME( STOREV, 'C' ) ) THEN
 279 *                    Skip any leading zeros.
 280                      DO LASTV = 1, I-1
 281                         IF( V( LASTV, I ).NE.ZERO ) EXIT
 282                      END DO
 283                      DO J = I+1, K
 284                         T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
 285                      END DO
 286                      J = MAX( LASTV, PREVLASTV )
 287 *
 288 *                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
 289 *
 290                      CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I,
 291      $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
 292      $                           1, ONE, T( I+1, I ), 1 )
 293                   ELSE
 294 *                    Skip any leading zeros.
 295                      DO LASTV = 1, I-1
 296                         IF( V( I, LASTV ).NE.ZERO ) EXIT
 297                      END DO
 298                      DO J = I+1, K
 299                         T( J, I ) = -TAU( I ) * V( J, N-K+I )
 300                      END DO
 301                      J = MAX( LASTV, PREVLASTV )
 302 *
 303 *                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
 304 *
 305                      CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
 306      $                           V( I+1, J ), LDV, V( I, J ), LDV,
 307      $                           ONE, T( I+1, I ), LDT )
 308                   END IF
 309 *
 310 *                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
 311 *
 312                   CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
 313      $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
 314                   IF( I.GT.1 ) THEN
 315                      PREVLASTV = MIN( PREVLASTV, LASTV )
 316                   ELSE
 317                      PREVLASTV = LASTV
 318                   END IF
 319                END IF
 320                T( I, I ) = TAU( I )
 321             END IF
 322          END DO
 323       END IF
 324       RETURN
 325 *
 326 *     End of CLARFT
 327 *
 328       END