LAPACK 3.3.0

zlatrd.f

Go to the documentation of this file.
00001       SUBROUTINE ZLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
00002 *
00003 *  -- LAPACK auxiliary 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            LDA, LDW, N, NB
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   E( * )
00014       COMPLEX*16         A( LDA, * ), TAU( * ), W( LDW, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to
00021 *  Hermitian tridiagonal form by a unitary similarity
00022 *  transformation Q' * A * Q, and returns the matrices V and W which are
00023 *  needed to apply the transformation to the unreduced part of A.
00024 *
00025 *  If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a
00026 *  matrix, of which the upper triangle is supplied;
00027 *  if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a
00028 *  matrix, of which the lower triangle is supplied.
00029 *
00030 *  This is an auxiliary routine called by ZHETRD.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  UPLO    (input) CHARACTER*1
00036 *          Specifies whether the upper or lower triangular part of the
00037 *          Hermitian matrix A is stored:
00038 *          = 'U': Upper triangular
00039 *          = 'L': Lower triangular
00040 *
00041 *  N       (input) INTEGER
00042 *          The order of the matrix A.
00043 *
00044 *  NB      (input) INTEGER
00045 *          The number of rows and columns to be reduced.
00046 *
00047 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00048 *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00049 *          n-by-n upper triangular part of A contains the upper
00050 *          triangular part of the matrix A, and the strictly lower
00051 *          triangular part of A is not referenced.  If UPLO = 'L', the
00052 *          leading n-by-n lower triangular part of A contains the lower
00053 *          triangular part of the matrix A, and the strictly upper
00054 *          triangular part of A is not referenced.
00055 *          On exit:
00056 *          if UPLO = 'U', the last NB columns have been reduced to
00057 *            tridiagonal form, with the diagonal elements overwriting
00058 *            the diagonal elements of A; the elements above the diagonal
00059 *            with the array TAU, represent the unitary matrix Q as a
00060 *            product of elementary reflectors;
00061 *          if UPLO = 'L', the first NB columns have been reduced to
00062 *            tridiagonal form, with the diagonal elements overwriting
00063 *            the diagonal elements of A; the elements below the diagonal
00064 *            with the array TAU, represent the  unitary matrix Q as a
00065 *            product of elementary reflectors.
00066 *          See Further Details.
00067 *
00068 *  LDA     (input) INTEGER
00069 *          The leading dimension of the array A.  LDA >= max(1,N).
00070 *
00071 *  E       (output) DOUBLE PRECISION array, dimension (N-1)
00072 *          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
00073 *          elements of the last NB columns of the reduced matrix;
00074 *          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
00075 *          the first NB columns of the reduced matrix.
00076 *
00077 *  TAU     (output) COMPLEX*16 array, dimension (N-1)
00078 *          The scalar factors of the elementary reflectors, stored in
00079 *          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
00080 *          See Further Details.
00081 *
00082 *  W       (output) COMPLEX*16 array, dimension (LDW,NB)
00083 *          The n-by-nb matrix W required to update the unreduced part
00084 *          of A.
00085 *
00086 *  LDW     (input) INTEGER
00087 *          The leading dimension of the array W. LDW >= max(1,N).
00088 *
00089 *  Further Details
00090 *  ===============
00091 *
00092 *  If UPLO = 'U', the matrix Q is represented as a product of elementary
00093 *  reflectors
00094 *
00095 *     Q = H(n) H(n-1) . . . H(n-nb+1).
00096 *
00097 *  Each H(i) has the form
00098 *
00099 *     H(i) = I - tau * v * v'
00100 *
00101 *  where tau is a complex scalar, and v is a complex vector with
00102 *  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
00103 *  and tau in TAU(i-1).
00104 *
00105 *  If UPLO = 'L', the matrix Q is represented as a product of elementary
00106 *  reflectors
00107 *
00108 *     Q = H(1) H(2) . . . H(nb).
00109 *
00110 *  Each H(i) has the form
00111 *
00112 *     H(i) = I - tau * v * v'
00113 *
00114 *  where tau is a complex scalar, and v is a complex vector with
00115 *  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
00116 *  and tau in TAU(i).
00117 *
00118 *  The elements of the vectors v together form the n-by-nb matrix V
00119 *  which is needed, with W, to apply the transformation to the unreduced
00120 *  part of the matrix, using a Hermitian rank-2k update of the form:
00121 *  A := A - V*W' - W*V'.
00122 *
00123 *  The contents of A on exit are illustrated by the following examples
00124 *  with n = 5 and nb = 2:
00125 *
00126 *  if UPLO = 'U':                       if UPLO = 'L':
00127 *
00128 *    (  a   a   a   v4  v5 )              (  d                  )
00129 *    (      a   a   v4  v5 )              (  1   d              )
00130 *    (          a   1   v5 )              (  v1  1   a          )
00131 *    (              d   1  )              (  v1  v2  a   a      )
00132 *    (                  d  )              (  v1  v2  a   a   a  )
00133 *
00134 *  where d denotes a diagonal element of the reduced matrix, a denotes
00135 *  an element of the original matrix that is unchanged, and vi denotes
00136 *  an element of the vector defining H(i).
00137 *
00138 *  =====================================================================
00139 *
00140 *     .. Parameters ..
00141       COMPLEX*16         ZERO, ONE, HALF
00142       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
00143      $                   ONE = ( 1.0D+0, 0.0D+0 ),
00144      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
00145 *     ..
00146 *     .. Local Scalars ..
00147       INTEGER            I, IW
00148       COMPLEX*16         ALPHA
00149 *     ..
00150 *     .. External Subroutines ..
00151       EXTERNAL           ZAXPY, ZGEMV, ZHEMV, ZLACGV, ZLARFG, ZSCAL
00152 *     ..
00153 *     .. External Functions ..
00154       LOGICAL            LSAME
00155       COMPLEX*16         ZDOTC
00156       EXTERNAL           LSAME, ZDOTC
00157 *     ..
00158 *     .. Intrinsic Functions ..
00159       INTRINSIC          DBLE, MIN
00160 *     ..
00161 *     .. Executable Statements ..
00162 *
00163 *     Quick return if possible
00164 *
00165       IF( N.LE.0 )
00166      $   RETURN
00167 *
00168       IF( LSAME( UPLO, 'U' ) ) THEN
00169 *
00170 *        Reduce last NB columns of upper triangle
00171 *
00172          DO 10 I = N, N - NB + 1, -1
00173             IW = I - N + NB
00174             IF( I.LT.N ) THEN
00175 *
00176 *              Update A(1:i,i)
00177 *
00178                A( I, I ) = DBLE( A( I, I ) )
00179                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
00180                CALL ZGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
00181      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
00182                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
00183                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
00184                CALL ZGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
00185      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
00186                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
00187                A( I, I ) = DBLE( A( I, I ) )
00188             END IF
00189             IF( I.GT.1 ) THEN
00190 *
00191 *              Generate elementary reflector H(i) to annihilate
00192 *              A(1:i-2,i)
00193 *
00194                ALPHA = A( I-1, I )
00195                CALL ZLARFG( I-1, ALPHA, A( 1, I ), 1, TAU( I-1 ) )
00196                E( I-1 ) = ALPHA
00197                A( I-1, I ) = ONE
00198 *
00199 *              Compute W(1:i-1,i)
00200 *
00201                CALL ZHEMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
00202      $                     ZERO, W( 1, IW ), 1 )
00203                IF( I.LT.N ) THEN
00204                   CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE,
00205      $                        W( 1, IW+1 ), LDW, A( 1, I ), 1, ZERO,
00206      $                        W( I+1, IW ), 1 )
00207                   CALL ZGEMV( 'No transpose', I-1, N-I, -ONE,
00208      $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
00209      $                        W( 1, IW ), 1 )
00210                   CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE,
00211      $                        A( 1, I+1 ), LDA, A( 1, I ), 1, ZERO,
00212      $                        W( I+1, IW ), 1 )
00213                   CALL ZGEMV( 'No transpose', I-1, N-I, -ONE,
00214      $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
00215      $                        W( 1, IW ), 1 )
00216                END IF
00217                CALL ZSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
00218                ALPHA = -HALF*TAU( I-1 )*ZDOTC( I-1, W( 1, IW ), 1,
00219      $                 A( 1, I ), 1 )
00220                CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
00221             END IF
00222 *
00223    10    CONTINUE
00224       ELSE
00225 *
00226 *        Reduce first NB columns of lower triangle
00227 *
00228          DO 20 I = 1, NB
00229 *
00230 *           Update A(i:n,i)
00231 *
00232             A( I, I ) = DBLE( A( I, I ) )
00233             CALL ZLACGV( I-1, W( I, 1 ), LDW )
00234             CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
00235      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
00236             CALL ZLACGV( I-1, W( I, 1 ), LDW )
00237             CALL ZLACGV( I-1, A( I, 1 ), LDA )
00238             CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
00239      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
00240             CALL ZLACGV( I-1, A( I, 1 ), LDA )
00241             A( I, I ) = DBLE( A( I, I ) )
00242             IF( I.LT.N ) THEN
00243 *
00244 *              Generate elementary reflector H(i) to annihilate
00245 *              A(i+2:n,i)
00246 *
00247                ALPHA = A( I+1, I )
00248                CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1,
00249      $                      TAU( I ) )
00250                E( I ) = ALPHA
00251                A( I+1, I ) = ONE
00252 *
00253 *              Compute W(i+1:n,i)
00254 *
00255                CALL ZHEMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
00256      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
00257                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
00258      $                     W( I+1, 1 ), LDW, A( I+1, I ), 1, ZERO,
00259      $                     W( 1, I ), 1 )
00260                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
00261      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
00262                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
00263      $                     A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO,
00264      $                     W( 1, I ), 1 )
00265                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
00266      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
00267                CALL ZSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
00268                ALPHA = -HALF*TAU( I )*ZDOTC( N-I, W( I+1, I ), 1,
00269      $                 A( I+1, I ), 1 )
00270                CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
00271             END IF
00272 *
00273    20    CONTINUE
00274       END IF
00275 *
00276       RETURN
00277 *
00278 *     End of ZLATRD
00279 *
00280       END
 All Files Functions