102
103
104
105
106
107
108 INTEGER INFO, K, LDA, N
109
110
111 INTEGER ISEED( 4 )
112 DOUBLE PRECISION D( * )
113 COMPLEX*16 A( LDA, * ), WORK( * )
114
115
116
117
118
119 COMPLEX*16 ZERO, ONE, HALF
120 parameter( zero = ( 0.0d+0, 0.0d+0 ),
121 $ one = ( 1.0d+0, 0.0d+0 ),
122 $ half = ( 0.5d+0, 0.0d+0 ) )
123
124
125 INTEGER I, J
126 DOUBLE PRECISION WN
127 COMPLEX*16 ALPHA, TAU, WA, WB
128
129
132
133
134 DOUBLE PRECISION DZNRM2
135 COMPLEX*16 ZDOTC
137
138
139 INTRINSIC abs, dble, dconjg, max
140
141
142
143
144
145 info = 0
146 IF( n.LT.0 ) THEN
147 info = -1
148 ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
149 info = -2
150 ELSE IF( lda.LT.max( 1, n ) ) THEN
151 info = -5
152 END IF
153 IF( info.LT.0 ) THEN
154 CALL xerbla(
'ZLAGHE', -info )
155 RETURN
156 END IF
157
158
159
160 DO 20 j = 1, n
161 DO 10 i = j + 1, n
162 a( i, j ) = zero
163 10 CONTINUE
164 20 CONTINUE
165 DO 30 i = 1, n
166 a( i, i ) = d( i )
167 30 CONTINUE
168
169
170
171 DO 40 i = n - 1, 1, -1
172
173
174
175 CALL zlarnv( 3, iseed, n-i+1, work )
176 wn =
dznrm2( n-i+1, work, 1 )
177 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
178 IF( wn.EQ.zero ) THEN
179 tau = zero
180 ELSE
181 wb = work( 1 ) + wa
182 CALL zscal( n-i, one / wb, work( 2 ), 1 )
183 work( 1 ) = one
184 tau = dble( wb / wa )
185 END IF
186
187
188
189
190
191
192 CALL zhemv(
'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
193 $ work( n+1 ), 1 )
194
195
196
197 alpha = -half*tau*
zdotc( n-i+1, work( n+1 ), 1, work, 1 )
198 CALL zaxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
199
200
201
202 CALL zher2(
'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
203 $ a( i, i ), lda )
204 40 CONTINUE
205
206
207
208 DO 60 i = 1, n - 1 - k
209
210
211
212 wn =
dznrm2( n-k-i+1, a( k+i, i ), 1 )
213 wa = ( wn / abs( a( k+i, i ) ) )*a( k+i, i )
214 IF( wn.EQ.zero ) THEN
215 tau = zero
216 ELSE
217 wb = a( k+i, i ) + wa
218 CALL zscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
219 a( k+i, i ) = one
220 tau = dble( wb / wa )
221 END IF
222
223
224
225 CALL zgemv(
'Conjugate transpose', n-k-i+1, k-1, one,
226 $ a( k+i, i+1 ), lda, a( k+i, i ), 1, zero, work, 1 )
227 CALL zgerc( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
228 $ a( k+i, i+1 ), lda )
229
230
231
232
233
234 CALL zhemv(
'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
235 $ a( k+i, i ), 1, zero, work, 1 )
236
237
238
239 alpha = -half*tau*
zdotc( n-k-i+1, work, 1, a( k+i, i ), 1 )
240 CALL zaxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
241
242
243
244 CALL zher2(
'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
245 $ a( k+i, k+i ), lda )
246
247 a( k+i, i ) = -wa
248 DO 50 j = k + i + 1, n
249 a( j, i ) = zero
250 50 CONTINUE
251 60 CONTINUE
252
253
254
255 DO 80 j = 1, n
256 DO 70 i = j + 1, n
257 a( j, i ) = dconjg( a( i, j ) )
258 70 CONTINUE
259 80 CONTINUE
260 RETURN
261
262
263
subroutine xerbla(srname, info)
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
complex *16 function zdotc(n, zx, incx, zy, incy)
ZDOTC
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
subroutine zgerc(m, n, alpha, x, incx, y, incy, a, lda)
ZGERC
subroutine zhemv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
ZHEMV
subroutine zher2(uplo, n, alpha, x, incx, y, incy, a, lda)
ZHER2
subroutine zlarnv(idist, iseed, n, x)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
real(wp) function dznrm2(n, x, incx)
DZNRM2
subroutine zscal(n, za, zx, incx)
ZSCAL