Go to the documentation of this file.
00001 *> \brief \b SLARFT
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLARFT + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarft.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          DIRECT, STOREV
00025 *       INTEGER            K, LDT, LDV, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> SLARFT forms the triangular factor T of a real block reflector H
00038 *> of order n, which is defined as a product of k elementary reflectors.
00039 *>
00040 *> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
00041 *>
00042 *> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
00043 *>
00044 *> If STOREV = 'C', the vector which defines the elementary reflector
00045 *> H(i) is stored in the i-th column of the array V, and
00046 *>
00047 *>    H  =  I - V * T * V**T
00048 *>
00049 *> If STOREV = 'R', the vector which defines the elementary reflector
00050 *> H(i) is stored in the i-th row of the array V, and
00051 *>
00052 *>    H  =  I - V**T * T * V
00053 *> \endverbatim
00054 *
00055 *  Arguments:
00056 *  ==========
00057 *
00058 *> \param[in] DIRECT
00059 *> \verbatim
00060 *>          DIRECT is CHARACTER*1
00061 *>          Specifies the order in which the elementary reflectors are
00062 *>          multiplied to form the block reflector:
00063 *>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
00064 *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
00065 *> \endverbatim
00066 *>
00067 *> \param[in] STOREV
00068 *> \verbatim
00069 *>          STOREV is CHARACTER*1
00070 *>          Specifies how the vectors which define the elementary
00071 *>          reflectors are stored (see also Further Details):
00072 *>          = 'C': columnwise
00073 *>          = 'R': rowwise
00074 *> \endverbatim
00075 *>
00076 *> \param[in] N
00077 *> \verbatim
00078 *>          N is INTEGER
00079 *>          The order of the block reflector H. N >= 0.
00080 *> \endverbatim
00081 *>
00082 *> \param[in] K
00083 *> \verbatim
00084 *>          K is INTEGER
00085 *>          The order of the triangular factor T (= the number of
00086 *>          elementary reflectors). K >= 1.
00087 *> \endverbatim
00088 *>
00089 *> \param[in,out] V
00090 *> \verbatim
00091 *>          V is REAL array, dimension
00092 *>                               (LDV,K) if STOREV = 'C'
00093 *>                               (LDV,N) if STOREV = 'R'
00094 *>          The matrix V. See further details.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] LDV
00098 *> \verbatim
00099 *>          LDV is INTEGER
00100 *>          The leading dimension of the array V.
00101 *>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] TAU
00105 *> \verbatim
00106 *>          TAU is REAL array, dimension (K)
00107 *>          TAU(i) must contain the scalar factor of the elementary
00108 *>          reflector H(i).
00109 *> \endverbatim
00110 *>
00111 *> \param[out] T
00112 *> \verbatim
00113 *>          T is REAL array, dimension (LDT,K)
00114 *>          The k by k triangular factor T of the block reflector.
00115 *>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
00116 *>          lower triangular. The rest of the array is not used.
00117 *> \endverbatim
00118 *>
00119 *> \param[in] LDT
00120 *> \verbatim
00121 *>          LDT is INTEGER
00122 *>          The leading dimension of the array T. LDT >= K.
00123 *> \endverbatim
00124 *
00125 *  Authors:
00126 *  ========
00127 *
00128 *> \author Univ. of Tennessee 
00129 *> \author Univ. of California Berkeley 
00130 *> \author Univ. of Colorado Denver 
00131 *> \author NAG Ltd. 
00132 *
00133 *> \date November 2011
00134 *
00135 *> \ingroup realOTHERauxiliary
00136 *
00137 *> \par Further Details:
00138 *  =====================
00139 *>
00140 *> \verbatim
00141 *>
00142 *>  The shape of the matrix V and the storage of the vectors which define
00143 *>  the H(i) is best illustrated by the following example with n = 5 and
00144 *>  k = 3. The elements equal to 1 are not stored; the corresponding
00145 *>  array elements are modified but restored on exit. The rest of the
00146 *>  array is not used.
00147 *>
00148 *>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
00149 *>
00150 *>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
00151 *>                   ( v1  1    )                     (     1 v2 v2 v2 )
00152 *>                   ( v1 v2  1 )                     (        1 v3 v3 )
00153 *>                   ( v1 v2 v3 )
00154 *>                   ( v1 v2 v3 )
00155 *>
00156 *>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
00157 *>
00158 *>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
00159 *>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
00160 *>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
00161 *>                   (     1 v3 )
00162 *>                   (        1 )
00163 *> \endverbatim
00164 *>
00165 *  =====================================================================
00166       SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
00167 *
00168 *  -- LAPACK auxiliary routine (version 3.4.0) --
00169 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00170 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00171 *     November 2011
00172 *
00173 *     .. Scalar Arguments ..
00174       CHARACTER          DIRECT, STOREV
00175       INTEGER            K, LDT, LDV, N
00176 *     ..
00177 *     .. Array Arguments ..
00178       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
00179 *     ..
00180 *
00181 *  =====================================================================
00182 *
00183 *     .. Parameters ..
00184       REAL               ONE, ZERO
00185       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00186 *     ..
00187 *     .. Local Scalars ..
00188       INTEGER            I, J, PREVLASTV, LASTV
00189       REAL               VII
00190 *     ..
00191 *     .. External Subroutines ..
00192       EXTERNAL           SGEMV, STRMV
00193 *     ..
00194 *     .. External Functions ..
00195       LOGICAL            LSAME
00196       EXTERNAL           LSAME
00197 *     ..
00198 *     .. Executable Statements ..
00199 *
00200 *     Quick return if possible
00201 *
00202       IF( N.EQ.0 )
00203      $   RETURN
00204 *
00205       IF( LSAME( DIRECT, 'F' ) ) THEN
00206          PREVLASTV = N
00207          DO 20 I = 1, K
00208             PREVLASTV = MAX( I, PREVLASTV )
00209             IF( TAU( I ).EQ.ZERO ) THEN
00210 *
00211 *              H(i)  =  I
00212 *
00213                DO 10 J = 1, I
00214                   T( J, I ) = ZERO
00215    10          CONTINUE
00216             ELSE
00217 *
00218 *              general case
00219 *
00220                VII = V( I, I )
00221                V( I, I ) = ONE
00222                IF( LSAME( STOREV, 'C' ) ) THEN
00223 !                 Skip any trailing zeros.
00224                   DO LASTV = N, I+1, -1
00225                      IF( V( LASTV, I ).NE.ZERO ) EXIT
00226                   END DO
00227                   J = MIN( LASTV, PREVLASTV )
00228 *
00229 *                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
00230 *
00231                   CALL SGEMV( 'Transpose', J-I+1, I-1, -TAU( I ),
00232      $                        V( I, 1 ), LDV, V( I, I ), 1, ZERO,
00233      $                        T( 1, I ), 1 )
00234                ELSE
00235 !                 Skip any trailing zeros.
00236                   DO LASTV = N, I+1, -1
00237                      IF( V( I, LASTV ).NE.ZERO ) EXIT
00238                   END DO
00239                   J = MIN( LASTV, PREVLASTV )
00240 *
00241 *                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
00242 *
00243                   CALL SGEMV( 'No transpose', I-1, J-I+1, -TAU( I ),
00244      $                        V( 1, I ), LDV, V( I, I ), LDV, ZERO,
00245      $                        T( 1, I ), 1 )
00246                END IF
00247                V( I, I ) = VII
00248 *
00249 *              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
00250 *
00251                CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
00252      $                     LDT, T( 1, I ), 1 )
00253                T( I, I ) = TAU( I )
00254                IF( I.GT.1 ) THEN
00255                   PREVLASTV = MAX( PREVLASTV, LASTV )
00256                ELSE
00257                   PREVLASTV = LASTV
00258                END IF
00259             END IF
00260    20    CONTINUE
00261       ELSE
00262          PREVLASTV = 1
00263          DO 40 I = K, 1, -1
00264             IF( TAU( I ).EQ.ZERO ) THEN
00265 *
00266 *              H(i)  =  I
00267 *
00268                DO 30 J = I, K
00269                   T( J, I ) = ZERO
00270    30          CONTINUE
00271             ELSE
00272 *
00273 *              general case
00274 *
00275                IF( I.LT.K ) THEN
00276                   IF( LSAME( STOREV, 'C' ) ) THEN
00277                      VII = V( N-K+I, I )
00278                      V( N-K+I, I ) = ONE
00279 !                    Skip any leading zeros.
00280                      DO LASTV = 1, I-1
00281                         IF( V( LASTV, I ).NE.ZERO ) EXIT
00282                      END DO
00283                      J = MAX( LASTV, PREVLASTV )
00284 *
00285 *                    T(i+1:k,i) :=
00286 *                            - tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
00287 *
00288                      CALL SGEMV( 'Transpose', N-K+I-J+1, K-I, -TAU( I ),
00289      $                           V( J, I+1 ), LDV, V( J, I ), 1, ZERO,
00290      $                           T( I+1, I ), 1 )
00291                      V( N-K+I, I ) = VII
00292                   ELSE
00293                      VII = V( I, N-K+I )
00294                      V( I, N-K+I ) = ONE
00295 !                    Skip any leading zeros.
00296                      DO LASTV = 1, I-1
00297                         IF( V( I, LASTV ).NE.ZERO ) EXIT
00298                      END DO
00299                      J = MAX( LASTV, PREVLASTV )
00300 *
00301 *                    T(i+1:k,i) :=
00302 *                            - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
00303 *
00304                      CALL SGEMV( 'No transpose', K-I, N-K+I-J+1,
00305      $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
00306      $                    ZERO, T( I+1, I ), 1 )
00307                      V( I, N-K+I ) = VII
00308                   END IF
00309 *
00310 *                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
00311 *
00312                   CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
00313      $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
00314                   IF( I.GT.1 ) THEN
00315                      PREVLASTV = MIN( PREVLASTV, LASTV )
00316                   ELSE
00317                      PREVLASTV = LASTV
00318                   END IF
00319                END IF
00320                T( I, I ) = TAU( I )
00321             END IF
00322    40    CONTINUE
00323       END IF
00324       RETURN
00325 *
00326 *     End of SLARFT
00327 *
00328       END