LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlarft.f
Go to the documentation of this file.
00001 *> \brief \b ZLARFT
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZLARFT + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarft.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarft.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarft.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZLARFT( 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*16         T( LDT, * ), TAU( * ), V( LDV, * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> ZLARFT 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*16 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*16 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*16 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 complex16OTHERauxiliary
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 ZLARFT( 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*16         T( LDT, * ), TAU( * ), V( LDV, * )
00177 *     ..
00178 *
00179 *  =====================================================================
00180 *
00181 *     .. Parameters ..
00182       COMPLEX*16         ONE, ZERO
00183       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
00184      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
00185 *     ..
00186 *     .. Local Scalars ..
00187       INTEGER            I, J, PREVLASTV, LASTV
00188 *     ..
00189 *     .. External Subroutines ..
00190       EXTERNAL           ZGEMV, ZLACGV, ZTRMV
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 ZGEMV( 'Conjugate transpose', J-I, I-1,
00231      $                        -TAU( I ), V( I+1, 1 ), LDV, 
00232      $                        V( I+1, I ), 1, ONE, T( 1, I ), 1 )
00233                ELSE
00234 *                 Skip any trailing zeros.
00235                   DO LASTV = N, I+1, -1
00236                      IF( V( I, LASTV ).NE.ZERO ) EXIT
00237                   END DO
00238                   DO J = 1, I-1
00239                      T( J, I ) = -TAU( I ) * V( J , I )
00240                   END DO                     
00241                   J = MIN( LASTV, PREVLASTV )
00242 *
00243 *                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
00244 *
00245                   CALL ZGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
00246      $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
00247      $                        ONE, T( 1, I ), LDT )                  
00248                END IF
00249 *
00250 *              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
00251 *
00252                CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
00253      $                     LDT, T( 1, I ), 1 )
00254                T( I, I ) = TAU( I )
00255                IF( I.GT.1 ) THEN
00256                   PREVLASTV = MAX( PREVLASTV, LASTV )
00257                ELSE
00258                   PREVLASTV = LASTV
00259                END IF
00260              END IF
00261          END DO
00262       ELSE
00263          PREVLASTV = 1
00264          DO I = K, 1, -1
00265             IF( TAU( I ).EQ.ZERO ) THEN
00266 *
00267 *              H(i)  =  I
00268 *
00269                DO J = I, K
00270                   T( J, I ) = ZERO
00271                END DO
00272             ELSE
00273 *
00274 *              general case
00275 *
00276                IF( I.LT.K ) THEN
00277                   IF( LSAME( STOREV, 'C' ) ) THEN
00278 *                    Skip any leading zeros.
00279                      DO LASTV = 1, I-1
00280                         IF( V( LASTV, I ).NE.ZERO ) EXIT
00281                      END DO
00282                      DO J = I+1, K
00283                         T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
00284                      END DO                        
00285                      J = MAX( LASTV, PREVLASTV )
00286 *
00287 *                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
00288 *
00289                      CALL ZGEMV( 'Conjugate transpose', N-K+I-J, K-I,
00290      $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
00291      $                           1, ONE, T( I+1, I ), 1 )
00292                   ELSE
00293 *                    Skip any leading zeros.
00294                      DO LASTV = 1, I-1
00295                         IF( V( I, LASTV ).NE.ZERO ) EXIT
00296                      END DO
00297                      DO J = I+1, K
00298                         T( J, I ) = -TAU( I ) * V( J, N-K+I )
00299                      END DO                                           
00300                      J = MAX( LASTV, PREVLASTV )
00301 *
00302 *                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
00303 *
00304                      CALL ZGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
00305      $                           V( I+1, J ), LDV, V( I, J ), LDV,
00306      $                           ONE, T( I+1, I ), LDT )                     
00307                   END IF
00308 *
00309 *                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
00310 *
00311                   CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
00312      $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
00313                   IF( I.GT.1 ) THEN
00314                      PREVLASTV = MIN( PREVLASTV, LASTV )
00315                   ELSE
00316                      PREVLASTV = LASTV
00317                   END IF
00318                END IF
00319                T( I, I ) = TAU( I )
00320             END IF
00321          END DO
00322       END IF
00323       RETURN
00324 *
00325 *     End of ZLARFT
00326 *
00327       END
 All Files Functions