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,out] 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 November 2011
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; 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 ZLARFT( 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       COMPLEX*16         T( LDT, * ), TAU( * ), V( LDV, * )
00179 *     ..
00180 *
00181 *  =====================================================================
00182 *
00183 *     .. Parameters ..
00184       COMPLEX*16         ONE, ZERO
00185       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
00186      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
00187 *     ..
00188 *     .. Local Scalars ..
00189       INTEGER            I, J, PREVLASTV, LASTV
00190       COMPLEX*16         VII
00191 *     ..
00192 *     .. External Subroutines ..
00193       EXTERNAL           ZGEMV, ZLACGV, ZTRMV
00194 *     ..
00195 *     .. External Functions ..
00196       LOGICAL            LSAME
00197       EXTERNAL           LSAME
00198 *     ..
00199 *     .. Executable Statements ..
00200 *
00201 *     Quick return if possible
00202 *
00203       IF( N.EQ.0 )
00204      $   RETURN
00205 *
00206       IF( LSAME( DIRECT, 'F' ) ) THEN
00207          PREVLASTV = N
00208          DO 20 I = 1, K
00209             PREVLASTV = MAX( PREVLASTV, I )
00210             IF( TAU( I ).EQ.ZERO ) THEN
00211 *
00212 *              H(i)  =  I
00213 *
00214                DO 10 J = 1, I
00215                   T( J, I ) = ZERO
00216    10          CONTINUE
00217             ELSE
00218 *
00219 *              general case
00220 *
00221                VII = V( I, I )
00222                V( I, I ) = ONE
00223                IF( LSAME( STOREV, 'C' ) ) THEN
00224 !                 Skip any trailing zeros.
00225                   DO LASTV = N, I+1, -1
00226                      IF( V( LASTV, I ).NE.ZERO ) EXIT
00227                   END DO
00228                   J = MIN( LASTV, PREVLASTV )
00229 *
00230 *                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
00231 *
00232                   CALL ZGEMV( 'Conjugate transpose', J-I+1, I-1,
00233      $                        -TAU( I ), V( I, 1 ), LDV, V( I, I ), 1,
00234      $                        ZERO, T( 1, I ), 1 )
00235                ELSE
00236 !                 Skip any trailing zeros.
00237                   DO LASTV = N, I+1, -1
00238                      IF( V( I, LASTV ).NE.ZERO ) EXIT
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)**H
00243 *
00244                   IF( I.LT.J )
00245      $               CALL ZLACGV( J-I, V( I, I+1 ), LDV )
00246                   CALL ZGEMV( 'No transpose', I-1, J-I+1, -TAU( I ),
00247      $                        V( 1, I ), LDV, V( I, I ), LDV, ZERO,
00248      $                        T( 1, I ), 1 )
00249                   IF( I.LT.J )
00250      $               CALL ZLACGV( J-I, V( I, I+1 ), LDV )
00251                END IF
00252                V( I, I ) = VII
00253 *
00254 *              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
00255 *
00256                CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
00257      $                     LDT, T( 1, I ), 1 )
00258                T( I, I ) = TAU( I )
00259                IF( I.GT.1 ) THEN
00260                   PREVLASTV = MAX( PREVLASTV, LASTV )
00261                ELSE
00262                   PREVLASTV = LASTV
00263                END IF
00264              END IF
00265    20    CONTINUE
00266       ELSE
00267          PREVLASTV = 1
00268          DO 40 I = K, 1, -1
00269             IF( TAU( I ).EQ.ZERO ) THEN
00270 *
00271 *              H(i)  =  I
00272 *
00273                DO 30 J = I, K
00274                   T( J, I ) = ZERO
00275    30          CONTINUE
00276             ELSE
00277 *
00278 *              general case
00279 *
00280                IF( I.LT.K ) THEN
00281                   IF( LSAME( STOREV, 'C' ) ) THEN
00282                      VII = V( N-K+I, I )
00283                      V( N-K+I, I ) = ONE
00284 !                    Skip any leading zeros.
00285                      DO LASTV = 1, I-1
00286                         IF( V( LASTV, I ).NE.ZERO ) EXIT
00287                      END DO
00288                      J = MAX( LASTV, PREVLASTV )
00289 *
00290 *                    T(i+1:k,i) :=
00291 *                            - tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
00292 *
00293                      CALL ZGEMV( 'Conjugate transpose', N-K+I-J+1, K-I,
00294      $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
00295      $                           1, ZERO, T( I+1, I ), 1 )
00296                      V( N-K+I, I ) = VII
00297                   ELSE
00298                      VII = V( I, N-K+I )
00299                      V( I, N-K+I ) = ONE
00300 !                    Skip any leading zeros.
00301                      DO LASTV = 1, I-1
00302                         IF( V( I, LASTV ).NE.ZERO ) EXIT
00303                      END DO
00304                      J = MAX( LASTV, PREVLASTV )
00305 *
00306 *                    T(i+1:k,i) :=
00307 *                            - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
00308 *
00309                      CALL ZLACGV( N-K+I-1-J+1, V( I, J ), LDV )
00310                      CALL ZGEMV( 'No transpose', K-I, N-K+I-J+1,
00311      $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
00312      $                    ZERO, T( I+1, I ), 1 )
00313                      CALL ZLACGV( N-K+I-1-J+1, V( I, J ), LDV )
00314                      V( I, N-K+I ) = VII
00315                   END IF
00316 *
00317 *                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
00318 *
00319                   CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
00320      $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
00321                   IF( I.GT.1 ) THEN
00322                      PREVLASTV = MIN( PREVLASTV, LASTV )
00323                   ELSE
00324                      PREVLASTV = LASTV
00325                   END IF
00326                END IF
00327                T( I, I ) = TAU( I )
00328             END IF
00329    40    CONTINUE
00330       END IF
00331       RETURN
00332 *
00333 *     End of ZLARFT
00334 *
00335       END