ground/gcs/src/libs/eigen/lapack/dlarft.f

   1 *> \brief \b DLARFT
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLARFT + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarft.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarft.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarft.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLARFT( 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 *       DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DLARFT forms the triangular factor T of a real 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**T
  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**T * 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 April 2012
 134 *
 135 *> \ingroup doubleOTHERauxiliary
 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 DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
 165 *
 166 *  -- LAPACK auxiliary routine (version 3.4.1) --
 167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 169 *     April 2012
 170 *
 171 *     .. Scalar Arguments ..
 172       CHARACTER          DIRECT, STOREV
 173       INTEGER            K, LDT, LDV, N
 174 *     ..
 175 *     .. Array Arguments ..
 176       DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
 177 *     ..
 178 *
 179 *  =====================================================================
 180 *
 181 *     .. Parameters ..
 182       DOUBLE PRECISION   ONE, ZERO
 183       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 184 *     ..
 185 *     .. Local Scalars ..
 186       INTEGER            I, J, PREVLASTV, LASTV
 187 *     ..
 188 *     .. External Subroutines ..
 189       EXTERNAL           DGEMV, DTRMV
 190 *     ..
 191 *     .. External Functions ..
 192       LOGICAL            LSAME
 193       EXTERNAL           LSAME
 194 *     ..
 195 *     .. Executable Statements ..
 196 *
 197 *     Quick return if possible
 198 *
 199       IF( N.EQ.0 )
 200      $   RETURN
 201 *
 202       IF( LSAME( DIRECT, 'F' ) ) THEN
 203          PREVLASTV = N
 204          DO I = 1, K
 205             PREVLASTV = MAX( I, PREVLASTV )
 206             IF( TAU( I ).EQ.ZERO ) THEN
 207 *
 208 *              H(i)  =  I
 209 *
 210                DO J = 1, I
 211                   T( J, I ) = ZERO
 212                END DO
 213             ELSE
 214 *
 215 *              general case
 216 *
 217                IF( LSAME( STOREV, 'C' ) ) THEN
 218 *                 Skip any trailing zeros.
 219                   DO LASTV = N, I+1, -1
 220                      IF( V( LASTV, I ).NE.ZERO ) EXIT
 221                   END DO
 222                   DO J = 1, I-1
 223                      T( J, I ) = -TAU( I ) * V( I , J )
 224                   END DO
 225                   J = MIN( LASTV, PREVLASTV )
 226 *
 227 *                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
 228 *
 229                   CALL DGEMV( 'Transpose', J-I, I-1, -TAU( I ),
 230      $                        V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE,
 231      $                        T( 1, I ), 1 )
 232                ELSE
 233 *                 Skip any trailing zeros.
 234                   DO LASTV = N, I+1, -1
 235                      IF( V( I, LASTV ).NE.ZERO ) EXIT
 236                   END DO
 237                   DO J = 1, I-1
 238                      T( J, I ) = -TAU( I ) * V( J , I )
 239                   END DO
 240                   J = MIN( LASTV, PREVLASTV )
 241 *
 242 *                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
 243 *
 244                   CALL DGEMV( 'No transpose', I-1, J-I, -TAU( I ),
 245      $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV, ONE,
 246      $                        T( 1, I ), 1 )
 247                END IF
 248 *
 249 *              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
 250 *
 251                CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
 252      $                     LDT, T( 1, I ), 1 )
 253                T( I, I ) = TAU( I )
 254                IF( I.GT.1 ) THEN
 255                   PREVLASTV = MAX( PREVLASTV, LASTV )
 256                ELSE
 257                   PREVLASTV = LASTV
 258                END IF
 259             END IF
 260          END DO
 261       ELSE
 262          PREVLASTV = 1
 263          DO I = K, 1, -1
 264             IF( TAU( I ).EQ.ZERO ) THEN
 265 *
 266 *              H(i)  =  I
 267 *
 268                DO J = I, K
 269                   T( J, I ) = ZERO
 270                END DO
 271             ELSE
 272 *
 273 *              general case
 274 *
 275                IF( I.LT.K ) THEN
 276                   IF( LSAME( STOREV, 'C' ) ) THEN
 277 *                    Skip any leading zeros.
 278                      DO LASTV = 1, I-1
 279                         IF( V( LASTV, I ).NE.ZERO ) EXIT
 280                      END DO
 281                      DO J = I+1, K
 282                         T( J, I ) = -TAU( I ) * V( N-K+I , J )
 283                      END DO
 284                      J = MAX( LASTV, PREVLASTV )
 285 *
 286 *                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
 287 *
 288                      CALL DGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ),
 289      $                           V( J, I+1 ), LDV, V( J, I ), 1, ONE,
 290      $                           T( I+1, I ), 1 )
 291                   ELSE
 292 *                    Skip any leading zeros.
 293                      DO LASTV = 1, I-1
 294                         IF( V( I, LASTV ).NE.ZERO ) EXIT
 295                      END DO
 296                      DO J = I+1, K
 297                         T( J, I ) = -TAU( I ) * V( J, N-K+I )
 298                      END DO
 299                      J = MAX( LASTV, PREVLASTV )
 300 *
 301 *                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
 302 *
 303                      CALL DGEMV( 'No transpose', K-I, N-K+I-J,
 304      $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
 305      $                    ONE, T( I+1, I ), 1 )
 306                   END IF
 307 *
 308 *                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
 309 *
 310                   CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
 311      $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
 312                   IF( I.GT.1 ) THEN
 313                      PREVLASTV = MIN( PREVLASTV, LASTV )
 314                   ELSE
 315                      PREVLASTV = LASTV
 316                   END IF
 317                END IF
 318                T( I, I ) = TAU( I )
 319             END IF
 320          END DO
 321       END IF
 322       RETURN
 323 *
 324 *     End of DLARFT
 325 *
 326       END