LAPACK 3.3.1 Linear Algebra PACKage

slarft.f

Go to the documentation of this file.
```00001       SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
00002       IMPLICIT NONE
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.3.1) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *  -- April 2011                                                      --
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          DIRECT, STOREV
00011       INTEGER            K, LDT, LDV, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SLARFT forms the triangular factor T of a real block reflector H
00021 *  of order n, which is defined as a product of k elementary reflectors.
00022 *
00023 *  If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
00024 *
00025 *  If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
00026 *
00027 *  If STOREV = 'C', the vector which defines the elementary reflector
00028 *  H(i) is stored in the i-th column of the array V, and
00029 *
00030 *     H  =  I - V * T * V**T
00031 *
00032 *  If STOREV = 'R', the vector which defines the elementary reflector
00033 *  H(i) is stored in the i-th row of the array V, and
00034 *
00035 *     H  =  I - V**T * T * V
00036 *
00037 *  Arguments
00038 *  =========
00039 *
00040 *  DIRECT  (input) CHARACTER*1
00041 *          Specifies the order in which the elementary reflectors are
00042 *          multiplied to form the block reflector:
00043 *          = 'F': H = H(1) H(2) . . . H(k) (Forward)
00044 *          = 'B': H = H(k) . . . H(2) H(1) (Backward)
00045 *
00046 *  STOREV  (input) CHARACTER*1
00047 *          Specifies how the vectors which define the elementary
00049 *          = 'C': columnwise
00050 *          = 'R': rowwise
00051 *
00052 *  N       (input) INTEGER
00053 *          The order of the block reflector H. N >= 0.
00054 *
00055 *  K       (input) INTEGER
00056 *          The order of the triangular factor T (= the number of
00057 *          elementary reflectors). K >= 1.
00058 *
00059 *  V       (input/output) REAL array, dimension
00060 *                               (LDV,K) if STOREV = 'C'
00061 *                               (LDV,N) if STOREV = 'R'
00062 *          The matrix V. See further details.
00063 *
00064 *  LDV     (input) INTEGER
00065 *          The leading dimension of the array V.
00066 *          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
00067 *
00068 *  TAU     (input) REAL array, dimension (K)
00069 *          TAU(i) must contain the scalar factor of the elementary
00070 *          reflector H(i).
00071 *
00072 *  T       (output) REAL array, dimension (LDT,K)
00073 *          The k by k triangular factor T of the block reflector.
00074 *          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
00075 *          lower triangular. The rest of the array is not used.
00076 *
00077 *  LDT     (input) INTEGER
00078 *          The leading dimension of the array T. LDT >= K.
00079 *
00080 *  Further Details
00081 *  ===============
00082 *
00083 *  The shape of the matrix V and the storage of the vectors which define
00084 *  the H(i) is best illustrated by the following example with n = 5 and
00085 *  k = 3. The elements equal to 1 are not stored; the corresponding
00086 *  array elements are modified but restored on exit. The rest of the
00087 *  array is not used.
00088 *
00089 *  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
00090 *
00091 *               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
00092 *                   ( v1  1    )                     (     1 v2 v2 v2 )
00093 *                   ( v1 v2  1 )                     (        1 v3 v3 )
00094 *                   ( v1 v2 v3 )
00095 *                   ( v1 v2 v3 )
00096 *
00097 *  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
00098 *
00099 *               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
00100 *                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
00101 *                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
00102 *                   (     1 v3 )
00103 *                   (        1 )
00104 *
00105 *  =====================================================================
00106 *
00107 *     .. Parameters ..
00108       REAL               ONE, ZERO
00109       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00110 *     ..
00111 *     .. Local Scalars ..
00112       INTEGER            I, J, PREVLASTV, LASTV
00113       REAL               VII
00114 *     ..
00115 *     .. External Subroutines ..
00116       EXTERNAL           SGEMV, STRMV
00117 *     ..
00118 *     .. External Functions ..
00119       LOGICAL            LSAME
00120       EXTERNAL           LSAME
00121 *     ..
00122 *     .. Executable Statements ..
00123 *
00124 *     Quick return if possible
00125 *
00126       IF( N.EQ.0 )
00127      \$   RETURN
00128 *
00129       IF( LSAME( DIRECT, 'F' ) ) THEN
00130          PREVLASTV = N
00131          DO 20 I = 1, K
00132             PREVLASTV = MAX( I, PREVLASTV )
00133             IF( TAU( I ).EQ.ZERO ) THEN
00134 *
00135 *              H(i)  =  I
00136 *
00137                DO 10 J = 1, I
00138                   T( J, I ) = ZERO
00139    10          CONTINUE
00140             ELSE
00141 *
00142 *              general case
00143 *
00144                VII = V( I, I )
00145                V( I, I ) = ONE
00146                IF( LSAME( STOREV, 'C' ) ) THEN
00147 !                 Skip any trailing zeros.
00148                   DO LASTV = N, I+1, -1
00149                      IF( V( LASTV, I ).NE.ZERO ) EXIT
00150                   END DO
00151                   J = MIN( LASTV, PREVLASTV )
00152 *
00153 *                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
00154 *
00155                   CALL SGEMV( 'Transpose', J-I+1, I-1, -TAU( I ),
00156      \$                        V( I, 1 ), LDV, V( I, I ), 1, ZERO,
00157      \$                        T( 1, I ), 1 )
00158                ELSE
00159 !                 Skip any trailing zeros.
00160                   DO LASTV = N, I+1, -1
00161                      IF( V( I, LASTV ).NE.ZERO ) EXIT
00162                   END DO
00163                   J = MIN( LASTV, PREVLASTV )
00164 *
00165 *                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
00166 *
00167                   CALL SGEMV( 'No transpose', I-1, J-I+1, -TAU( I ),
00168      \$                        V( 1, I ), LDV, V( I, I ), LDV, ZERO,
00169      \$                        T( 1, I ), 1 )
00170                END IF
00171                V( I, I ) = VII
00172 *
00173 *              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
00174 *
00175                CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
00176      \$                     LDT, T( 1, I ), 1 )
00177                T( I, I ) = TAU( I )
00178                IF( I.GT.1 ) THEN
00179                   PREVLASTV = MAX( PREVLASTV, LASTV )
00180                ELSE
00181                   PREVLASTV = LASTV
00182                END IF
00183             END IF
00184    20    CONTINUE
00185       ELSE
00186          PREVLASTV = 1
00187          DO 40 I = K, 1, -1
00188             IF( TAU( I ).EQ.ZERO ) THEN
00189 *
00190 *              H(i)  =  I
00191 *
00192                DO 30 J = I, K
00193                   T( J, I ) = ZERO
00194    30          CONTINUE
00195             ELSE
00196 *
00197 *              general case
00198 *
00199                IF( I.LT.K ) THEN
00200                   IF( LSAME( STOREV, 'C' ) ) THEN
00201                      VII = V( N-K+I, I )
00202                      V( N-K+I, I ) = ONE
00203 !                    Skip any leading zeros.
00204                      DO LASTV = 1, I-1
00205                         IF( V( LASTV, I ).NE.ZERO ) EXIT
00206                      END DO
00207                      J = MAX( LASTV, PREVLASTV )
00208 *
00209 *                    T(i+1:k,i) :=
00210 *                            - tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
00211 *
00212                      CALL SGEMV( 'Transpose', N-K+I-J+1, K-I, -TAU( I ),
00213      \$                           V( J, I+1 ), LDV, V( J, I ), 1, ZERO,
00214      \$                           T( I+1, I ), 1 )
00215                      V( N-K+I, I ) = VII
00216                   ELSE
00217                      VII = V( I, N-K+I )
00218                      V( I, N-K+I ) = ONE
00219 !                    Skip any leading zeros.
00220                      DO LASTV = 1, I-1
00221                         IF( V( I, LASTV ).NE.ZERO ) EXIT
00222                      END DO
00223                      J = MAX( LASTV, PREVLASTV )
00224 *
00225 *                    T(i+1:k,i) :=
00226 *                            - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
00227 *
00228                      CALL SGEMV( 'No transpose', K-I, N-K+I-J+1,
00229      \$                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
00230      \$                    ZERO, T( I+1, I ), 1 )
00231                      V( I, N-K+I ) = VII
00232                   END IF
00233 *
00234 *                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
00235 *
00236                   CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
00237      \$                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
00238                   IF( I.GT.1 ) THEN
00239                      PREVLASTV = MIN( PREVLASTV, LASTV )
00240                   ELSE
00241                      PREVLASTV = LASTV
00242                   END IF
00243                END IF
00244                T( I, I ) = TAU( I )
00245             END IF
00246    40    CONTINUE
00247       END IF
00248       RETURN
00249 *
00250 *     End of SLARFT
00251 *
00252       END
```