001:       SUBROUTINE ZLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            K, LDA, LDT, LDY, N, NB
010: *     ..
011: *     .. Array Arguments ..
012:       COMPLEX*16         A( LDA, * ), T( LDT, NB ), TAU( NB ),
013:      $                   Y( LDY, NB )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZLAHRD reduces the first NB columns of a complex general n-by-(n-k+1)
020: *  matrix A so that elements below the k-th subdiagonal are zero. The
021: *  reduction is performed by a unitary similarity transformation
022: *  Q' * A * Q. The routine returns the matrices V and T which determine
023: *  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
024: *
025: *  This is an OBSOLETE auxiliary routine. 
026: *  This routine will be 'deprecated' in a  future release.
027: *  Please use the new routine ZLAHR2 instead.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  N       (input) INTEGER
033: *          The order of the matrix A.
034: *
035: *  K       (input) INTEGER
036: *          The offset for the reduction. Elements below the k-th
037: *          subdiagonal in the first NB columns are reduced to zero.
038: *
039: *  NB      (input) INTEGER
040: *          The number of columns to be reduced.
041: *
042: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N-K+1)
043: *          On entry, the n-by-(n-k+1) general matrix A.
044: *          On exit, the elements on and above the k-th subdiagonal in
045: *          the first NB columns are overwritten with the corresponding
046: *          elements of the reduced matrix; the elements below the k-th
047: *          subdiagonal, with the array TAU, represent the matrix Q as a
048: *          product of elementary reflectors. The other columns of A are
049: *          unchanged. See Further Details.
050: *
051: *  LDA     (input) INTEGER
052: *          The leading dimension of the array A.  LDA >= max(1,N).
053: *
054: *  TAU     (output) COMPLEX*16 array, dimension (NB)
055: *          The scalar factors of the elementary reflectors. See Further
056: *          Details.
057: *
058: *  T       (output) COMPLEX*16 array, dimension (LDT,NB)
059: *          The upper triangular matrix T.
060: *
061: *  LDT     (input) INTEGER
062: *          The leading dimension of the array T.  LDT >= NB.
063: *
064: *  Y       (output) COMPLEX*16 array, dimension (LDY,NB)
065: *          The n-by-nb matrix Y.
066: *
067: *  LDY     (input) INTEGER
068: *          The leading dimension of the array Y. LDY >= max(1,N).
069: *
070: *  Further Details
071: *  ===============
072: *
073: *  The matrix Q is represented as a product of nb elementary reflectors
074: *
075: *     Q = H(1) H(2) . . . H(nb).
076: *
077: *  Each H(i) has the form
078: *
079: *     H(i) = I - tau * v * v'
080: *
081: *  where tau is a complex scalar, and v is a complex vector with
082: *  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
083: *  A(i+k+1:n,i), and tau in TAU(i).
084: *
085: *  The elements of the vectors v together form the (n-k+1)-by-nb matrix
086: *  V which is needed, with T and Y, to apply the transformation to the
087: *  unreduced part of the matrix, using an update of the form:
088: *  A := (I - V*T*V') * (A - Y*V').
089: *
090: *  The contents of A on exit are illustrated by the following example
091: *  with n = 7, k = 3 and nb = 2:
092: *
093: *     ( a   h   a   a   a )
094: *     ( a   h   a   a   a )
095: *     ( a   h   a   a   a )
096: *     ( h   h   a   a   a )
097: *     ( v1  h   a   a   a )
098: *     ( v1  v2  a   a   a )
099: *     ( v1  v2  a   a   a )
100: *
101: *  where a denotes an element of the original matrix A, h denotes a
102: *  modified element of the upper Hessenberg matrix H, and vi denotes an
103: *  element of the vector defining H(i).
104: *
105: *  =====================================================================
106: *
107: *     .. Parameters ..
108:       COMPLEX*16         ZERO, ONE
109:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
110:      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
111: *     ..
112: *     .. Local Scalars ..
113:       INTEGER            I
114:       COMPLEX*16         EI
115: *     ..
116: *     .. External Subroutines ..
117:       EXTERNAL           ZAXPY, ZCOPY, ZGEMV, ZLACGV, ZLARFG, ZSCAL,
118:      $                   ZTRMV
119: *     ..
120: *     .. Intrinsic Functions ..
121:       INTRINSIC          MIN
122: *     ..
123: *     .. Executable Statements ..
124: *
125: *     Quick return if possible
126: *
127:       IF( N.LE.1 )
128:      $   RETURN
129: *
130:       DO 10 I = 1, NB
131:          IF( I.GT.1 ) THEN
132: *
133: *           Update A(1:n,i)
134: *
135: *           Compute i-th column of A - Y * V'
136: *
137:             CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
138:             CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
139:      $                  A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
140:             CALL ZLACGV( I-1, A( K+I-1, 1 ), LDA )
141: *
142: *           Apply I - V * T' * V' to this column (call it b) from the
143: *           left, using the last column of T as workspace
144: *
145: *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
146: *                    ( V2 )             ( b2 )
147: *
148: *           where V1 is unit lower triangular
149: *
150: *           w := V1' * b1
151: *
152:             CALL ZCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
153:             CALL ZTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1,
154:      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
155: *
156: *           w := w + V2'*b2
157: *
158:             CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
159:      $                  A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE,
160:      $                  T( 1, NB ), 1 )
161: *
162: *           w := T'*w
163: *
164:             CALL ZTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1,
165:      $                  T, LDT, T( 1, NB ), 1 )
166: *
167: *           b2 := b2 - V2*w
168: *
169:             CALL ZGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
170:      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
171: *
172: *           b1 := b1 - V1*w
173: *
174:             CALL ZTRMV( 'Lower', 'No transpose', 'Unit', I-1,
175:      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
176:             CALL ZAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
177: *
178:             A( K+I-1, I-1 ) = EI
179:          END IF
180: *
181: *        Generate the elementary reflector H(i) to annihilate
182: *        A(k+i+1:n,i)
183: *
184:          EI = A( K+I, I )
185:          CALL ZLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1,
186:      $                TAU( I ) )
187:          A( K+I, I ) = ONE
188: *
189: *        Compute  Y(1:n,i)
190: *
191:          CALL ZGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
192:      $               A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
193:          CALL ZGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
194:      $               A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ),
195:      $               1 )
196:          CALL ZGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
197:      $               ONE, Y( 1, I ), 1 )
198:          CALL ZSCAL( N, TAU( I ), Y( 1, I ), 1 )
199: *
200: *        Compute T(1:i,i)
201: *
202:          CALL ZSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
203:          CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
204:      $               T( 1, I ), 1 )
205:          T( I, I ) = TAU( I )
206: *
207:    10 CONTINUE
208:       A( K+NB, NB ) = EI
209: *
210:       RETURN
211: *
212: *     End of ZLAHRD
213: *
214:       END
215: