LAPACK 3.3.1 Linear Algebra PACKage

# zhetrd.f

Go to the documentation of this file.
```00001       SUBROUTINE ZHETRD( UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, LDA, LWORK, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), E( * )
00014       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZHETRD reduces a complex Hermitian matrix A to real symmetric
00021 *  tridiagonal form T by a unitary similarity transformation:
00022 *  Q**H * A * Q = T.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  UPLO    (input) CHARACTER*1
00028 *          = 'U':  Upper triangle of A is stored;
00029 *          = 'L':  Lower triangle of A is stored.
00030 *
00031 *  N       (input) INTEGER
00032 *          The order of the matrix A.  N >= 0.
00033 *
00034 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00035 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00036 *          N-by-N upper triangular part of A contains the upper
00037 *          triangular part of the matrix A, and the strictly lower
00038 *          triangular part of A is not referenced.  If UPLO = 'L', the
00039 *          leading N-by-N lower triangular part of A contains the lower
00040 *          triangular part of the matrix A, and the strictly upper
00041 *          triangular part of A is not referenced.
00042 *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
00043 *          of A are overwritten by the corresponding elements of the
00044 *          tridiagonal matrix T, and the elements above the first
00045 *          superdiagonal, with the array TAU, represent the unitary
00046 *          matrix Q as a product of elementary reflectors; if UPLO
00047 *          = 'L', the diagonal and first subdiagonal of A are over-
00048 *          written by the corresponding elements of the tridiagonal
00049 *          matrix T, and the elements below the first subdiagonal, with
00050 *          the array TAU, represent the unitary matrix Q as a product
00051 *          of elementary reflectors. See Further Details.
00052 *
00053 *  LDA     (input) INTEGER
00054 *          The leading dimension of the array A.  LDA >= max(1,N).
00055 *
00056 *  D       (output) DOUBLE PRECISION array, dimension (N)
00057 *          The diagonal elements of the tridiagonal matrix T:
00058 *          D(i) = A(i,i).
00059 *
00060 *  E       (output) DOUBLE PRECISION array, dimension (N-1)
00061 *          The off-diagonal elements of the tridiagonal matrix T:
00062 *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
00063 *
00064 *  TAU     (output) COMPLEX*16 array, dimension (N-1)
00065 *          The scalar factors of the elementary reflectors (see Further
00066 *          Details).
00067 *
00068 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
00069 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00070 *
00071 *  LWORK   (input) INTEGER
00072 *          The dimension of the array WORK.  LWORK >= 1.
00073 *          For optimum performance LWORK >= N*NB, where NB is the
00074 *          optimal blocksize.
00075 *
00076 *          If LWORK = -1, then a workspace query is assumed; the routine
00077 *          only calculates the optimal size of the WORK array, returns
00078 *          this value as the first entry of the WORK array, and no error
00079 *          message related to LWORK is issued by XERBLA.
00080 *
00081 *  INFO    (output) INTEGER
00082 *          = 0:  successful exit
00083 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00084 *
00085 *  Further Details
00086 *  ===============
00087 *
00088 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
00089 *  reflectors
00090 *
00091 *     Q = H(n-1) . . . H(2) H(1).
00092 *
00093 *  Each H(i) has the form
00094 *
00095 *     H(i) = I - tau * v * v**H
00096 *
00097 *  where tau is a complex scalar, and v is a complex vector with
00098 *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
00099 *  A(1:i-1,i+1), and tau in TAU(i).
00100 *
00101 *  If UPLO = 'L', the matrix Q is represented as a product of elementary
00102 *  reflectors
00103 *
00104 *     Q = H(1) H(2) . . . H(n-1).
00105 *
00106 *  Each H(i) has the form
00107 *
00108 *     H(i) = I - tau * v * v**H
00109 *
00110 *  where tau is a complex scalar, and v is a complex vector with
00111 *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
00112 *  and tau in TAU(i).
00113 *
00114 *  The contents of A on exit are illustrated by the following examples
00115 *  with n = 5:
00116 *
00117 *  if UPLO = 'U':                       if UPLO = 'L':
00118 *
00119 *    (  d   e   v2  v3  v4 )              (  d                  )
00120 *    (      d   e   v3  v4 )              (  e   d              )
00121 *    (          d   e   v4 )              (  v1  e   d          )
00122 *    (              d   e  )              (  v1  v2  e   d      )
00123 *    (                  d  )              (  v1  v2  v3  e   d  )
00124 *
00125 *  where d and e denote diagonal and off-diagonal elements of T, and vi
00126 *  denotes an element of the vector defining H(i).
00127 *
00128 *  =====================================================================
00129 *
00130 *     .. Parameters ..
00131       DOUBLE PRECISION   ONE
00132       PARAMETER          ( ONE = 1.0D+0 )
00133       COMPLEX*16         CONE
00134       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
00135 *     ..
00136 *     .. Local Scalars ..
00137       LOGICAL            LQUERY, UPPER
00138       INTEGER            I, IINFO, IWS, J, KK, LDWORK, LWKOPT, NB,
00139      \$                   NBMIN, NX
00140 *     ..
00141 *     .. External Subroutines ..
00142       EXTERNAL           XERBLA, ZHER2K, ZHETD2, ZLATRD
00143 *     ..
00144 *     .. Intrinsic Functions ..
00145       INTRINSIC          MAX
00146 *     ..
00147 *     .. External Functions ..
00148       LOGICAL            LSAME
00149       INTEGER            ILAENV
00150       EXTERNAL           LSAME, ILAENV
00151 *     ..
00152 *     .. Executable Statements ..
00153 *
00154 *     Test the input parameters
00155 *
00156       INFO = 0
00157       UPPER = LSAME( UPLO, 'U' )
00158       LQUERY = ( LWORK.EQ.-1 )
00159       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00160          INFO = -1
00161       ELSE IF( N.LT.0 ) THEN
00162          INFO = -2
00163       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00164          INFO = -4
00165       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
00166          INFO = -9
00167       END IF
00168 *
00169       IF( INFO.EQ.0 ) THEN
00170 *
00171 *        Determine the block size.
00172 *
00173          NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
00174          LWKOPT = N*NB
00175          WORK( 1 ) = LWKOPT
00176       END IF
00177 *
00178       IF( INFO.NE.0 ) THEN
00179          CALL XERBLA( 'ZHETRD', -INFO )
00180          RETURN
00181       ELSE IF( LQUERY ) THEN
00182          RETURN
00183       END IF
00184 *
00185 *     Quick return if possible
00186 *
00187       IF( N.EQ.0 ) THEN
00188          WORK( 1 ) = 1
00189          RETURN
00190       END IF
00191 *
00192       NX = N
00193       IWS = 1
00194       IF( NB.GT.1 .AND. NB.LT.N ) THEN
00195 *
00196 *        Determine when to cross over from blocked to unblocked code
00197 *        (last block is always handled by unblocked code).
00198 *
00199          NX = MAX( NB, ILAENV( 3, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
00200          IF( NX.LT.N ) THEN
00201 *
00202 *           Determine if workspace is large enough for blocked code.
00203 *
00204             LDWORK = N
00205             IWS = LDWORK*NB
00206             IF( LWORK.LT.IWS ) THEN
00207 *
00208 *              Not enough workspace to use optimal NB:  determine the
00209 *              minimum value of NB, and reduce NB or force use of
00210 *              unblocked code by setting NX = N.
00211 *
00212                NB = MAX( LWORK / LDWORK, 1 )
00213                NBMIN = ILAENV( 2, 'ZHETRD', UPLO, N, -1, -1, -1 )
00214                IF( NB.LT.NBMIN )
00215      \$            NX = N
00216             END IF
00217          ELSE
00218             NX = N
00219          END IF
00220       ELSE
00221          NB = 1
00222       END IF
00223 *
00224       IF( UPPER ) THEN
00225 *
00226 *        Reduce the upper triangle of A.
00227 *        Columns 1:kk are handled by the unblocked method.
00228 *
00229          KK = N - ( ( N-NX+NB-1 ) / NB )*NB
00230          DO 20 I = N - NB + 1, KK + 1, -NB
00231 *
00232 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
00233 *           matrix W which is needed to update the unreduced part of
00234 *           the matrix
00235 *
00236             CALL ZLATRD( UPLO, I+NB-1, NB, A, LDA, E, TAU, WORK,
00237      \$                   LDWORK )
00238 *
00239 *           Update the unreduced submatrix A(1:i-1,1:i-1), using an
00240 *           update of the form:  A := A - V*W**H - W*V**H
00241 *
00242             CALL ZHER2K( UPLO, 'No transpose', I-1, NB, -CONE,
00243      \$                   A( 1, I ), LDA, WORK, LDWORK, ONE, A, LDA )
00244 *
00245 *           Copy superdiagonal elements back into A, and diagonal
00246 *           elements into D
00247 *
00248             DO 10 J = I, I + NB - 1
00249                A( J-1, J ) = E( J-1 )
00250                D( J ) = A( J, J )
00251    10       CONTINUE
00252    20    CONTINUE
00253 *
00254 *        Use unblocked code to reduce the last or only block
00255 *
00256          CALL ZHETD2( UPLO, KK, A, LDA, D, E, TAU, IINFO )
00257       ELSE
00258 *
00259 *        Reduce the lower triangle of A
00260 *
00261          DO 40 I = 1, N - NX, NB
00262 *
00263 *           Reduce columns i:i+nb-1 to tridiagonal form and form the
00264 *           matrix W which is needed to update the unreduced part of
00265 *           the matrix
00266 *
00267             CALL ZLATRD( UPLO, N-I+1, NB, A( I, I ), LDA, E( I ),
00268      \$                   TAU( I ), WORK, LDWORK )
00269 *
00270 *           Update the unreduced submatrix A(i+nb:n,i+nb:n), using
00271 *           an update of the form:  A := A - V*W**H - W*V**H
00272 *
00273             CALL ZHER2K( UPLO, 'No transpose', N-I-NB+1, NB, -CONE,
00274      \$                   A( I+NB, I ), LDA, WORK( NB+1 ), LDWORK, ONE,
00275      \$                   A( I+NB, I+NB ), LDA )
00276 *
00277 *           Copy subdiagonal elements back into A, and diagonal
00278 *           elements into D
00279 *
00280             DO 30 J = I, I + NB - 1
00281                A( J+1, J ) = E( J )
00282                D( J ) = A( J, J )
00283    30       CONTINUE
00284    40    CONTINUE
00285 *
00286 *        Use unblocked code to reduce the last or only block
00287 *
00288          CALL ZHETD2( UPLO, N-I+1, A( I, I ), LDA, D( I ), E( I ),
00289      \$                TAU( I ), IINFO )
00290       END IF
00291 *
00292       WORK( 1 ) = LWKOPT
00293       RETURN
00294 *
00295 *     End of ZHETRD
00296 *
00297       END
```