Home | History | Annotate | Download | only in lapack
      1 *> \brief \b SLARFT
      2 *
      3 *  =========== DOCUMENTATION ===========
      4 *
      5 * Online html documentation available at 
      6 *            http://www.netlib.org/lapack/explore-html/ 
      7 *
      8 *> \htmlonly
      9 *> Download SLARFT + dependencies 
     10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f"> 
     11 *> [TGZ]</a> 
     12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f"> 
     13 *> [ZIP]</a> 
     14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarft.f"> 
     15 *> [TXT]</a>
     16 *> \endhtmlonly 
     17 *
     18 *  Definition:
     19 *  ===========
     20 *
     21 *       SUBROUTINE SLARFT( 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 *       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
     29 *       ..
     30 *  
     31 *
     32 *> \par Purpose:
     33 *  =============
     34 *>
     35 *> \verbatim
     36 *>
     37 *> SLARFT 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 REAL 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 REAL 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 REAL 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 realOTHERauxiliary
    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 SLARFT( 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       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
    177 *     ..
    178 *
    179 *  =====================================================================
    180 *
    181 *     .. Parameters ..
    182       REAL               ONE, ZERO
    183       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
    184 *     ..
    185 *     .. Local Scalars ..
    186       INTEGER            I, J, PREVLASTV, LASTV
    187 *     ..
    188 *     .. External Subroutines ..
    189       EXTERNAL           SGEMV, STRMV
    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 SGEMV( '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 SGEMV( 'No transpose', I-1, J-I, -TAU( I ),
    245      $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 
    246      $                        ONE, 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 STRMV( '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 SGEMV( '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 SGEMV( '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 STRMV( '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 SLARFT
    325 *
    326       END
    327