LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlatrd.f
Go to the documentation of this file.
00001 *> \brief \b ZLATRD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZLATRD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlatrd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlatrd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlatrd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            LDA, LDW, N, NB
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       DOUBLE PRECISION   E( * )
00029 *       COMPLEX*16         A( LDA, * ), TAU( * ), W( LDW, * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to
00039 *> Hermitian tridiagonal form by a unitary similarity
00040 *> transformation Q**H * A * Q, and returns the matrices V and W which are
00041 *> needed to apply the transformation to the unreduced part of A.
00042 *>
00043 *> If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a
00044 *> matrix, of which the upper triangle is supplied;
00045 *> if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a
00046 *> matrix, of which the lower triangle is supplied.
00047 *>
00048 *> This is an auxiliary routine called by ZHETRD.
00049 *> \endverbatim
00050 *
00051 *  Arguments:
00052 *  ==========
00053 *
00054 *> \param[in] UPLO
00055 *> \verbatim
00056 *>          UPLO is CHARACTER*1
00057 *>          Specifies whether the upper or lower triangular part of the
00058 *>          Hermitian matrix A is stored:
00059 *>          = 'U': Upper triangular
00060 *>          = 'L': Lower triangular
00061 *> \endverbatim
00062 *>
00063 *> \param[in] N
00064 *> \verbatim
00065 *>          N is INTEGER
00066 *>          The order of the matrix A.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] NB
00070 *> \verbatim
00071 *>          NB is INTEGER
00072 *>          The number of rows and columns to be reduced.
00073 *> \endverbatim
00074 *>
00075 *> \param[in,out] A
00076 *> \verbatim
00077 *>          A is COMPLEX*16 array, dimension (LDA,N)
00078 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00079 *>          n-by-n upper triangular part of A contains the upper
00080 *>          triangular part of the matrix A, and the strictly lower
00081 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00082 *>          leading n-by-n lower triangular part of A contains the lower
00083 *>          triangular part of the matrix A, and the strictly upper
00084 *>          triangular part of A is not referenced.
00085 *>          On exit:
00086 *>          if UPLO = 'U', the last NB columns have been reduced to
00087 *>            tridiagonal form, with the diagonal elements overwriting
00088 *>            the diagonal elements of A; the elements above the diagonal
00089 *>            with the array TAU, represent the unitary matrix Q as a
00090 *>            product of elementary reflectors;
00091 *>          if UPLO = 'L', the first NB columns have been reduced to
00092 *>            tridiagonal form, with the diagonal elements overwriting
00093 *>            the diagonal elements of A; the elements below the diagonal
00094 *>            with the array TAU, represent the  unitary matrix Q as a
00095 *>            product of elementary reflectors.
00096 *>          See Further Details.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] LDA
00100 *> \verbatim
00101 *>          LDA is INTEGER
00102 *>          The leading dimension of the array A.  LDA >= max(1,N).
00103 *> \endverbatim
00104 *>
00105 *> \param[out] E
00106 *> \verbatim
00107 *>          E is DOUBLE PRECISION array, dimension (N-1)
00108 *>          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
00109 *>          elements of the last NB columns of the reduced matrix;
00110 *>          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
00111 *>          the first NB columns of the reduced matrix.
00112 *> \endverbatim
00113 *>
00114 *> \param[out] TAU
00115 *> \verbatim
00116 *>          TAU is COMPLEX*16 array, dimension (N-1)
00117 *>          The scalar factors of the elementary reflectors, stored in
00118 *>          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
00119 *>          See Further Details.
00120 *> \endverbatim
00121 *>
00122 *> \param[out] W
00123 *> \verbatim
00124 *>          W is COMPLEX*16 array, dimension (LDW,NB)
00125 *>          The n-by-nb matrix W required to update the unreduced part
00126 *>          of A.
00127 *> \endverbatim
00128 *>
00129 *> \param[in] LDW
00130 *> \verbatim
00131 *>          LDW is INTEGER
00132 *>          The leading dimension of the array W. LDW >= max(1,N).
00133 *> \endverbatim
00134 *
00135 *  Authors:
00136 *  ========
00137 *
00138 *> \author Univ. of Tennessee 
00139 *> \author Univ. of California Berkeley 
00140 *> \author Univ. of Colorado Denver 
00141 *> \author NAG Ltd. 
00142 *
00143 *> \date November 2011
00144 *
00145 *> \ingroup complex16OTHERauxiliary
00146 *
00147 *> \par Further Details:
00148 *  =====================
00149 *>
00150 *> \verbatim
00151 *>
00152 *>  If UPLO = 'U', the matrix Q is represented as a product of elementary
00153 *>  reflectors
00154 *>
00155 *>     Q = H(n) H(n-1) . . . H(n-nb+1).
00156 *>
00157 *>  Each H(i) has the form
00158 *>
00159 *>     H(i) = I - tau * v * v**H
00160 *>
00161 *>  where tau is a complex scalar, and v is a complex vector with
00162 *>  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
00163 *>  and tau in TAU(i-1).
00164 *>
00165 *>  If UPLO = 'L', the matrix Q is represented as a product of elementary
00166 *>  reflectors
00167 *>
00168 *>     Q = H(1) H(2) . . . H(nb).
00169 *>
00170 *>  Each H(i) has the form
00171 *>
00172 *>     H(i) = I - tau * v * v**H
00173 *>
00174 *>  where tau is a complex scalar, and v is a complex vector with
00175 *>  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
00176 *>  and tau in TAU(i).
00177 *>
00178 *>  The elements of the vectors v together form the n-by-nb matrix V
00179 *>  which is needed, with W, to apply the transformation to the unreduced
00180 *>  part of the matrix, using a Hermitian rank-2k update of the form:
00181 *>  A := A - V*W**H - W*V**H.
00182 *>
00183 *>  The contents of A on exit are illustrated by the following examples
00184 *>  with n = 5 and nb = 2:
00185 *>
00186 *>  if UPLO = 'U':                       if UPLO = 'L':
00187 *>
00188 *>    (  a   a   a   v4  v5 )              (  d                  )
00189 *>    (      a   a   v4  v5 )              (  1   d              )
00190 *>    (          a   1   v5 )              (  v1  1   a          )
00191 *>    (              d   1  )              (  v1  v2  a   a      )
00192 *>    (                  d  )              (  v1  v2  a   a   a  )
00193 *>
00194 *>  where d denotes a diagonal element of the reduced matrix, a denotes
00195 *>  an element of the original matrix that is unchanged, and vi denotes
00196 *>  an element of the vector defining H(i).
00197 *> \endverbatim
00198 *>
00199 *  =====================================================================
00200       SUBROUTINE ZLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
00201 *
00202 *  -- LAPACK auxiliary routine (version 3.4.0) --
00203 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00204 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00205 *     November 2011
00206 *
00207 *     .. Scalar Arguments ..
00208       CHARACTER          UPLO
00209       INTEGER            LDA, LDW, N, NB
00210 *     ..
00211 *     .. Array Arguments ..
00212       DOUBLE PRECISION   E( * )
00213       COMPLEX*16         A( LDA, * ), TAU( * ), W( LDW, * )
00214 *     ..
00215 *
00216 *  =====================================================================
00217 *
00218 *     .. Parameters ..
00219       COMPLEX*16         ZERO, ONE, HALF
00220       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
00221      $                   ONE = ( 1.0D+0, 0.0D+0 ),
00222      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
00223 *     ..
00224 *     .. Local Scalars ..
00225       INTEGER            I, IW
00226       COMPLEX*16         ALPHA
00227 *     ..
00228 *     .. External Subroutines ..
00229       EXTERNAL           ZAXPY, ZGEMV, ZHEMV, ZLACGV, ZLARFG, ZSCAL
00230 *     ..
00231 *     .. External Functions ..
00232       LOGICAL            LSAME
00233       COMPLEX*16         ZDOTC
00234       EXTERNAL           LSAME, ZDOTC
00235 *     ..
00236 *     .. Intrinsic Functions ..
00237       INTRINSIC          DBLE, MIN
00238 *     ..
00239 *     .. Executable Statements ..
00240 *
00241 *     Quick return if possible
00242 *
00243       IF( N.LE.0 )
00244      $   RETURN
00245 *
00246       IF( LSAME( UPLO, 'U' ) ) THEN
00247 *
00248 *        Reduce last NB columns of upper triangle
00249 *
00250          DO 10 I = N, N - NB + 1, -1
00251             IW = I - N + NB
00252             IF( I.LT.N ) THEN
00253 *
00254 *              Update A(1:i,i)
00255 *
00256                A( I, I ) = DBLE( A( I, I ) )
00257                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
00258                CALL ZGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
00259      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
00260                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
00261                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
00262                CALL ZGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
00263      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
00264                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
00265                A( I, I ) = DBLE( A( I, I ) )
00266             END IF
00267             IF( I.GT.1 ) THEN
00268 *
00269 *              Generate elementary reflector H(i) to annihilate
00270 *              A(1:i-2,i)
00271 *
00272                ALPHA = A( I-1, I )
00273                CALL ZLARFG( I-1, ALPHA, A( 1, I ), 1, TAU( I-1 ) )
00274                E( I-1 ) = ALPHA
00275                A( I-1, I ) = ONE
00276 *
00277 *              Compute W(1:i-1,i)
00278 *
00279                CALL ZHEMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
00280      $                     ZERO, W( 1, IW ), 1 )
00281                IF( I.LT.N ) THEN
00282                   CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE,
00283      $                        W( 1, IW+1 ), LDW, A( 1, I ), 1, ZERO,
00284      $                        W( I+1, IW ), 1 )
00285                   CALL ZGEMV( 'No transpose', I-1, N-I, -ONE,
00286      $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
00287      $                        W( 1, IW ), 1 )
00288                   CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE,
00289      $                        A( 1, I+1 ), LDA, A( 1, I ), 1, ZERO,
00290      $                        W( I+1, IW ), 1 )
00291                   CALL ZGEMV( 'No transpose', I-1, N-I, -ONE,
00292      $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
00293      $                        W( 1, IW ), 1 )
00294                END IF
00295                CALL ZSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
00296                ALPHA = -HALF*TAU( I-1 )*ZDOTC( I-1, W( 1, IW ), 1,
00297      $                 A( 1, I ), 1 )
00298                CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
00299             END IF
00300 *
00301    10    CONTINUE
00302       ELSE
00303 *
00304 *        Reduce first NB columns of lower triangle
00305 *
00306          DO 20 I = 1, NB
00307 *
00308 *           Update A(i:n,i)
00309 *
00310             A( I, I ) = DBLE( A( I, I ) )
00311             CALL ZLACGV( I-1, W( I, 1 ), LDW )
00312             CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
00313      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
00314             CALL ZLACGV( I-1, W( I, 1 ), LDW )
00315             CALL ZLACGV( I-1, A( I, 1 ), LDA )
00316             CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
00317      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
00318             CALL ZLACGV( I-1, A( I, 1 ), LDA )
00319             A( I, I ) = DBLE( A( I, I ) )
00320             IF( I.LT.N ) THEN
00321 *
00322 *              Generate elementary reflector H(i) to annihilate
00323 *              A(i+2:n,i)
00324 *
00325                ALPHA = A( I+1, I )
00326                CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1,
00327      $                      TAU( I ) )
00328                E( I ) = ALPHA
00329                A( I+1, I ) = ONE
00330 *
00331 *              Compute W(i+1:n,i)
00332 *
00333                CALL ZHEMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
00334      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
00335                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
00336      $                     W( I+1, 1 ), LDW, A( I+1, I ), 1, ZERO,
00337      $                     W( 1, I ), 1 )
00338                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
00339      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
00340                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
00341      $                     A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO,
00342      $                     W( 1, I ), 1 )
00343                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
00344      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
00345                CALL ZSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
00346                ALPHA = -HALF*TAU( I )*ZDOTC( N-I, W( I+1, I ), 1,
00347      $                 A( I+1, I ), 1 )
00348                CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
00349             END IF
00350 *
00351    20    CONTINUE
00352       END IF
00353 *
00354       RETURN
00355 *
00356 *     End of ZLATRD
00357 *
00358       END
 All Files Functions