LAPACK 3.3.0

zhetd2.f

Go to the documentation of this file.
00001       SUBROUTINE ZHETD2( UPLO, N, A, LDA, D, E, TAU, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, LDA, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   D( * ), E( * )
00014       COMPLEX*16         A( LDA, * ), TAU( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZHETD2 reduces a complex Hermitian matrix A to real symmetric
00021 *  tridiagonal form T by a unitary similarity transformation:
00022 *  Q' * A * Q = T.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  UPLO    (input) CHARACTER*1
00028 *          Specifies whether the upper or lower triangular part of the
00029 *          Hermitian matrix A is stored:
00030 *          = 'U':  Upper triangular
00031 *          = 'L':  Lower triangular
00032 *
00033 *  N       (input) INTEGER
00034 *          The order of the matrix A.  N >= 0.
00035 *
00036 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00037 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00038 *          n-by-n upper triangular part of A contains the upper
00039 *          triangular part of the matrix A, and the strictly lower
00040 *          triangular part of A is not referenced.  If UPLO = 'L', the
00041 *          leading n-by-n lower triangular part of A contains the lower
00042 *          triangular part of the matrix A, and the strictly upper
00043 *          triangular part of A is not referenced.
00044 *          On exit, if UPLO = 'U', the diagonal and first superdiagonal
00045 *          of A are overwritten by the corresponding elements of the
00046 *          tridiagonal matrix T, and the elements above the first
00047 *          superdiagonal, with the array TAU, represent the unitary
00048 *          matrix Q as a product of elementary reflectors; if UPLO
00049 *          = 'L', the diagonal and first subdiagonal of A are over-
00050 *          written by the corresponding elements of the tridiagonal
00051 *          matrix T, and the elements below the first subdiagonal, with
00052 *          the array TAU, represent the unitary matrix Q as a product
00053 *          of elementary reflectors. See Further Details.
00054 *
00055 *  LDA     (input) INTEGER
00056 *          The leading dimension of the array A.  LDA >= max(1,N).
00057 *
00058 *  D       (output) DOUBLE PRECISION array, dimension (N)
00059 *          The diagonal elements of the tridiagonal matrix T:
00060 *          D(i) = A(i,i).
00061 *
00062 *  E       (output) DOUBLE PRECISION array, dimension (N-1)
00063 *          The off-diagonal elements of the tridiagonal matrix T:
00064 *          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
00065 *
00066 *  TAU     (output) COMPLEX*16 array, dimension (N-1)
00067 *          The scalar factors of the elementary reflectors (see Further
00068 *          Details).
00069 *
00070 *  INFO    (output) INTEGER
00071 *          = 0:  successful exit
00072 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00073 *
00074 *  Further Details
00075 *  ===============
00076 *
00077 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
00078 *  reflectors
00079 *
00080 *     Q = H(n-1) . . . H(2) H(1).
00081 *
00082 *  Each H(i) has the form
00083 *
00084 *     H(i) = I - tau * v * v'
00085 *
00086 *  where tau is a complex scalar, and v is a complex vector with
00087 *  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
00088 *  A(1:i-1,i+1), and tau in TAU(i).
00089 *
00090 *  If UPLO = 'L', the matrix Q is represented as a product of elementary
00091 *  reflectors
00092 *
00093 *     Q = H(1) H(2) . . . H(n-1).
00094 *
00095 *  Each H(i) has the form
00096 *
00097 *     H(i) = I - tau * v * v'
00098 *
00099 *  where tau is a complex scalar, and v is a complex vector with
00100 *  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
00101 *  and tau in TAU(i).
00102 *
00103 *  The contents of A on exit are illustrated by the following examples
00104 *  with n = 5:
00105 *
00106 *  if UPLO = 'U':                       if UPLO = 'L':
00107 *
00108 *    (  d   e   v2  v3  v4 )              (  d                  )
00109 *    (      d   e   v3  v4 )              (  e   d              )
00110 *    (          d   e   v4 )              (  v1  e   d          )
00111 *    (              d   e  )              (  v1  v2  e   d      )
00112 *    (                  d  )              (  v1  v2  v3  e   d  )
00113 *
00114 *  where d and e denote diagonal and off-diagonal elements of T, and vi
00115 *  denotes an element of the vector defining H(i).
00116 *
00117 *  =====================================================================
00118 *
00119 *     .. Parameters ..
00120       COMPLEX*16         ONE, ZERO, HALF
00121       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
00122      $                   ZERO = ( 0.0D+0, 0.0D+0 ),
00123      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
00124 *     ..
00125 *     .. Local Scalars ..
00126       LOGICAL            UPPER
00127       INTEGER            I
00128       COMPLEX*16         ALPHA, TAUI
00129 *     ..
00130 *     .. External Subroutines ..
00131       EXTERNAL           XERBLA, ZAXPY, ZHEMV, ZHER2, ZLARFG
00132 *     ..
00133 *     .. External Functions ..
00134       LOGICAL            LSAME
00135       COMPLEX*16         ZDOTC
00136       EXTERNAL           LSAME, ZDOTC
00137 *     ..
00138 *     .. Intrinsic Functions ..
00139       INTRINSIC          DBLE, MAX, MIN
00140 *     ..
00141 *     .. Executable Statements ..
00142 *
00143 *     Test the input parameters
00144 *
00145       INFO = 0
00146       UPPER = LSAME( UPLO, 'U' )
00147       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00148          INFO = -1
00149       ELSE IF( N.LT.0 ) THEN
00150          INFO = -2
00151       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00152          INFO = -4
00153       END IF
00154       IF( INFO.NE.0 ) THEN
00155          CALL XERBLA( 'ZHETD2', -INFO )
00156          RETURN
00157       END IF
00158 *
00159 *     Quick return if possible
00160 *
00161       IF( N.LE.0 )
00162      $   RETURN
00163 *
00164       IF( UPPER ) THEN
00165 *
00166 *        Reduce the upper triangle of A
00167 *
00168          A( N, N ) = DBLE( A( N, N ) )
00169          DO 10 I = N - 1, 1, -1
00170 *
00171 *           Generate elementary reflector H(i) = I - tau * v * v'
00172 *           to annihilate A(1:i-1,i+1)
00173 *
00174             ALPHA = A( I, I+1 )
00175             CALL ZLARFG( I, ALPHA, A( 1, I+1 ), 1, TAUI )
00176             E( I ) = ALPHA
00177 *
00178             IF( TAUI.NE.ZERO ) THEN
00179 *
00180 *              Apply H(i) from both sides to A(1:i,1:i)
00181 *
00182                A( I, I+1 ) = ONE
00183 *
00184 *              Compute  x := tau * A * v  storing x in TAU(1:i)
00185 *
00186                CALL ZHEMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO,
00187      $                     TAU, 1 )
00188 *
00189 *              Compute  w := x - 1/2 * tau * (x'*v) * v
00190 *
00191                ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, A( 1, I+1 ), 1 )
00192                CALL ZAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 )
00193 *
00194 *              Apply the transformation as a rank-2 update:
00195 *                 A := A - v * w' - w * v'
00196 *
00197                CALL ZHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A,
00198      $                     LDA )
00199 *
00200             ELSE
00201                A( I, I ) = DBLE( A( I, I ) )
00202             END IF
00203             A( I, I+1 ) = E( I )
00204             D( I+1 ) = A( I+1, I+1 )
00205             TAU( I ) = TAUI
00206    10    CONTINUE
00207          D( 1 ) = A( 1, 1 )
00208       ELSE
00209 *
00210 *        Reduce the lower triangle of A
00211 *
00212          A( 1, 1 ) = DBLE( A( 1, 1 ) )
00213          DO 20 I = 1, N - 1
00214 *
00215 *           Generate elementary reflector H(i) = I - tau * v * v'
00216 *           to annihilate A(i+2:n,i)
00217 *
00218             ALPHA = A( I+1, I )
00219             CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAUI )
00220             E( I ) = ALPHA
00221 *
00222             IF( TAUI.NE.ZERO ) THEN
00223 *
00224 *              Apply H(i) from both sides to A(i+1:n,i+1:n)
00225 *
00226                A( I+1, I ) = ONE
00227 *
00228 *              Compute  x := tau * A * v  storing y in TAU(i:n-1)
00229 *
00230                CALL ZHEMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA,
00231      $                     A( I+1, I ), 1, ZERO, TAU( I ), 1 )
00232 *
00233 *              Compute  w := x - 1/2 * tau * (x'*v) * v
00234 *
00235                ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, A( I+1, I ),
00236      $                 1 )
00237                CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 )
00238 *
00239 *              Apply the transformation as a rank-2 update:
00240 *                 A := A - v * w' - w * v'
00241 *
00242                CALL ZHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1,
00243      $                     A( I+1, I+1 ), LDA )
00244 *
00245             ELSE
00246                A( I+1, I+1 ) = DBLE( A( I+1, I+1 ) )
00247             END IF
00248             A( I+1, I ) = E( I )
00249             D( I ) = A( I, I )
00250             TAU( I ) = TAUI
00251    20    CONTINUE
00252          D( N ) = A( N, N )
00253       END IF
00254 *
00255       RETURN
00256 *
00257 *     End of ZHETD2
00258 *
00259       END
 All Files Functions