LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slarft.f
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] 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 April 2012
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.
00145 *>
00146 *>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
00147 *>
00148 *>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
00149 *>                   ( v1  1    )                     (     1 v2 v2 v2 )
00150 *>                   ( v1 v2  1 )                     (        1 v3 v3 )
00151 *>                   ( v1 v2 v3 )
00152 *>                   ( v1 v2 v3 )
00153 *>
00154 *>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
00155 *>
00156 *>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
00157 *>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
00158 *>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
00159 *>                   (     1 v3 )
00160 *>                   (        1 )
00161 *> \endverbatim
00162 *>
00163 *  =====================================================================
00164       SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
00165 *
00166 *  -- LAPACK auxiliary routine (version 3.4.1) --
00167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00169 *     April 2012
00170 *
00171 *     .. Scalar Arguments ..
00172       CHARACTER          DIRECT, STOREV
00173       INTEGER            K, LDT, LDV, N
00174 *     ..
00175 *     .. Array Arguments ..
00176       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
00177 *     ..
00178 *
00179 *  =====================================================================
00180 *
00181 *     .. Parameters ..
00182       REAL               ONE, ZERO
00183       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00184 *     ..
00185 *     .. Local Scalars ..
00186       INTEGER            I, J, PREVLASTV, LASTV
00187 *     ..
00188 *     .. External Subroutines ..
00189       EXTERNAL           SGEMV, STRMV
00190 *     ..
00191 *     .. External Functions ..
00192       LOGICAL            LSAME
00193       EXTERNAL           LSAME
00194 *     ..
00195 *     .. Executable Statements ..
00196 *
00197 *     Quick return if possible
00198 *
00199       IF( N.EQ.0 )
00200      $   RETURN
00201 *
00202       IF( LSAME( DIRECT, 'F' ) ) THEN
00203          PREVLASTV = N
00204          DO I = 1, K
00205             PREVLASTV = MAX( I, PREVLASTV )
00206             IF( TAU( I ).EQ.ZERO ) THEN
00207 *
00208 *              H(i)  =  I
00209 *
00210                DO J = 1, I
00211                   T( J, I ) = ZERO
00212                END DO
00213             ELSE
00214 *
00215 *              general case
00216 *
00217                IF( LSAME( STOREV, 'C' ) ) THEN
00218 *                 Skip any trailing zeros.
00219                   DO LASTV = N, I+1, -1
00220                      IF( V( LASTV, I ).NE.ZERO ) EXIT
00221                   END DO
00222                   DO J = 1, I-1
00223                      T( J, I ) = -TAU( I ) * V( I , J )
00224                   END DO   
00225                   J = MIN( LASTV, PREVLASTV )
00226 *
00227 *                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
00228 *
00229                   CALL SGEMV( 'Transpose', J-I, I-1, -TAU( I ),
00230      $                        V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE,
00231      $                        T( 1, I ), 1 )
00232                ELSE
00233 *                 Skip any trailing zeros.
00234                   DO LASTV = N, I+1, -1
00235                      IF( V( I, LASTV ).NE.ZERO ) EXIT
00236                   END DO
00237                   DO J = 1, I-1
00238                      T( J, I ) = -TAU( I ) * V( J , I )
00239                   END DO   
00240                   J = MIN( LASTV, PREVLASTV )
00241 *
00242 *                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
00243 *
00244                   CALL SGEMV( 'No transpose', I-1, J-I, -TAU( I ),
00245      $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 
00246      $                        ONE, T( 1, I ), 1 )
00247                END IF
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          END DO
00261       ELSE
00262          PREVLASTV = 1
00263          DO I = K, 1, -1
00264             IF( TAU( I ).EQ.ZERO ) THEN
00265 *
00266 *              H(i)  =  I
00267 *
00268                DO J = I, K
00269                   T( J, I ) = ZERO
00270                END DO
00271             ELSE
00272 *
00273 *              general case
00274 *
00275                IF( I.LT.K ) THEN
00276                   IF( LSAME( STOREV, 'C' ) ) THEN
00277 *                    Skip any leading zeros.
00278                      DO LASTV = 1, I-1
00279                         IF( V( LASTV, I ).NE.ZERO ) EXIT
00280                      END DO
00281                      DO J = I+1, K
00282                         T( J, I ) = -TAU( I ) * V( N-K+I , J )
00283                      END DO   
00284                      J = MAX( LASTV, PREVLASTV )
00285 *
00286 *                    T(i+1:k,i) = -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, K-I, -TAU( I ),
00289      $                           V( J, I+1 ), LDV, V( J, I ), 1, ONE,
00290      $                           T( I+1, I ), 1 )
00291                   ELSE
00292 *                    Skip any leading zeros.
00293                      DO LASTV = 1, I-1
00294                         IF( V( I, LASTV ).NE.ZERO ) EXIT
00295                      END DO
00296                      DO J = I+1, K
00297                         T( J, I ) = -TAU( I ) * V( J, N-K+I )
00298                      END DO   
00299                      J = MAX( LASTV, PREVLASTV )
00300 *
00301 *                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
00302 *
00303                      CALL SGEMV( 'No transpose', K-I, N-K+I-J,
00304      $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
00305      $                    ONE, T( I+1, I ), 1 )
00306                   END IF
00307 *
00308 *                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
00309 *
00310                   CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
00311      $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
00312                   IF( I.GT.1 ) THEN
00313                      PREVLASTV = MIN( PREVLASTV, LASTV )
00314                   ELSE
00315                      PREVLASTV = LASTV
00316                   END IF
00317                END IF
00318                T( I, I ) = TAU( I )
00319             END IF
00320          END DO
00321       END IF
00322       RETURN
00323 *
00324 *     End of SLARFT
00325 *
00326       END
 All Files Functions