LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clarft.f
Go to the documentation of this file.
00001 *> \brief \b CLARFT
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CLARFT + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CLARFT( 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 *       COMPLEX            T( LDT, * ), TAU( * ), V( LDV, * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> CLARFT forms the triangular factor T of a complex 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**H
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**H * 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 complexOTHERauxiliary
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 CLARFT( 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       COMPLEX            T( LDT, * ), TAU( * ), V( LDV, * )
00177 *     ..
00178 *
00179 *  =====================================================================
00180 *
00181 *     .. Parameters ..
00182       COMPLEX            ONE, ZERO
00183       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
00184      $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
00185 *     ..
00186 *     .. Local Scalars ..
00187       INTEGER            I, J, PREVLASTV, LASTV
00188 *     ..
00189 *     .. External Subroutines ..
00190       EXTERNAL           CGEMV, CLACGV, CTRMV
00191 *     ..
00192 *     .. External Functions ..
00193       LOGICAL            LSAME
00194       EXTERNAL           LSAME
00195 *     ..
00196 *     .. Executable Statements ..
00197 *
00198 *     Quick return if possible
00199 *
00200       IF( N.EQ.0 )
00201      $   RETURN
00202 *
00203       IF( LSAME( DIRECT, 'F' ) ) THEN
00204          PREVLASTV = N
00205          DO I = 1, K
00206             PREVLASTV = MAX( PREVLASTV, I )
00207             IF( TAU( I ).EQ.ZERO ) THEN
00208 *
00209 *              H(i)  =  I
00210 *
00211                DO J = 1, I
00212                   T( J, I ) = ZERO
00213                END DO
00214             ELSE
00215 *
00216 *              general case
00217 *
00218                IF( LSAME( STOREV, 'C' ) ) THEN
00219 *                 Skip any trailing zeros.
00220                   DO LASTV = N, I+1, -1
00221                      IF( V( LASTV, I ).NE.ZERO ) EXIT
00222                   END DO
00223                   DO J = 1, I-1
00224                      T( J, I ) = -TAU( I ) * CONJG( V( I , J ) )
00225                   END DO                     
00226                   J = MIN( LASTV, PREVLASTV )
00227 *
00228 *                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
00229 *
00230                   CALL CGEMV( 'Conjugate transpose', J-I, I-1,
00231      $                        -TAU( I ), V( I+1, 1 ), LDV, 
00232      $                        V( I+1, I ), 1,
00233      $                        ONE, 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                   DO J = 1, I-1
00240                      T( J, I ) = -TAU( I ) * V( J , I )
00241                   END DO                     
00242                   J = MIN( LASTV, PREVLASTV )
00243 *
00244 *                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
00245 *
00246                   CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
00247      $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
00248      $                        ONE, T( 1, I ), LDT )                  
00249                END IF
00250 *
00251 *              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
00252 *
00253                CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
00254      $                     LDT, T( 1, I ), 1 )
00255                T( I, I ) = TAU( I )
00256                IF( I.GT.1 ) THEN
00257                   PREVLASTV = MAX( PREVLASTV, LASTV )
00258                ELSE
00259                   PREVLASTV = LASTV
00260                END IF
00261             END IF
00262          END DO
00263       ELSE
00264          PREVLASTV = 1
00265          DO I = K, 1, -1
00266             IF( TAU( I ).EQ.ZERO ) THEN
00267 *
00268 *              H(i)  =  I
00269 *
00270                DO J = I, K
00271                   T( J, I ) = ZERO
00272                END DO
00273             ELSE
00274 *
00275 *              general case
00276 *
00277                IF( I.LT.K ) THEN
00278                   IF( LSAME( STOREV, 'C' ) ) THEN
00279 *                    Skip any leading zeros.
00280                      DO LASTV = 1, I-1
00281                         IF( V( LASTV, I ).NE.ZERO ) EXIT
00282                      END DO
00283                      DO J = I+1, K
00284                         T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
00285                      END DO                        
00286                      J = MAX( LASTV, PREVLASTV )
00287 *
00288 *                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
00289 *
00290                      CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I,
00291      $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
00292      $                           1, ONE, T( I+1, I ), 1 )
00293                   ELSE
00294 *                    Skip any leading zeros.
00295                      DO LASTV = 1, I-1
00296                         IF( V( I, LASTV ).NE.ZERO ) EXIT
00297                      END DO
00298                      DO J = I+1, K
00299                         T( J, I ) = -TAU( I ) * V( J, N-K+I )
00300                      END DO                      
00301                      J = MAX( LASTV, PREVLASTV )
00302 *
00303 *                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
00304 *
00305                      CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
00306      $                           V( I+1, J ), LDV, V( I, J ), LDV,
00307      $                           ONE, T( I+1, I ), LDT )                     
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 CTRMV( '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          END DO
00323       END IF
00324       RETURN
00325 *
00326 *     End of CLARFT
00327 *
00328       END
 All Files Functions