LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlaghe.f
Go to the documentation of this file.
1*> \brief \b ZLAGHE
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE ZLAGHE( N, K, D, A, LDA, ISEED, WORK, INFO )
12*
13* .. Scalar Arguments ..
14* INTEGER INFO, K, LDA, N
15* ..
16* .. Array Arguments ..
17* INTEGER ISEED( 4 )
18* DOUBLE PRECISION D( * )
19* COMPLEX*16 A( LDA, * ), WORK( * )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> ZLAGHE generates a complex hermitian matrix A, by pre- and post-
29*> multiplying a real diagonal matrix D with a random unitary matrix:
30*> A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
31*> unitary transformations.
32*> \endverbatim
33*
34* Arguments:
35* ==========
36*
37*> \param[in] N
38*> \verbatim
39*> N is INTEGER
40*> The order of the matrix A. N >= 0.
41*> \endverbatim
42*>
43*> \param[in] K
44*> \verbatim
45*> K is INTEGER
46*> The number of nonzero subdiagonals within the band of A.
47*> 0 <= K <= N-1.
48*> \endverbatim
49*>
50*> \param[in] D
51*> \verbatim
52*> D is DOUBLE PRECISION array, dimension (N)
53*> The diagonal elements of the diagonal matrix D.
54*> \endverbatim
55*>
56*> \param[out] A
57*> \verbatim
58*> A is COMPLEX*16 array, dimension (LDA,N)
59*> The generated n by n hermitian matrix A (the full matrix is
60*> stored).
61*> \endverbatim
62*>
63*> \param[in] LDA
64*> \verbatim
65*> LDA is INTEGER
66*> The leading dimension of the array A. LDA >= N.
67*> \endverbatim
68*>
69*> \param[in,out] ISEED
70*> \verbatim
71*> ISEED is INTEGER array, dimension (4)
72*> On entry, the seed of the random number generator; the array
73*> elements must be between 0 and 4095, and ISEED(4) must be
74*> odd.
75*> On exit, the seed is updated.
76*> \endverbatim
77*>
78*> \param[out] WORK
79*> \verbatim
80*> WORK is COMPLEX*16 array, dimension (2*N)
81*> \endverbatim
82*>
83*> \param[out] INFO
84*> \verbatim
85*> INFO is INTEGER
86*> = 0: successful exit
87*> < 0: if INFO = -i, the i-th argument had an illegal value
88*> \endverbatim
89*
90* Authors:
91* ========
92*
93*> \author Univ. of Tennessee
94*> \author Univ. of California Berkeley
95*> \author Univ. of Colorado Denver
96*> \author NAG Ltd.
97*
98*> \ingroup complex16_matgen
99*
100* =====================================================================
101 SUBROUTINE zlaghe( N, K, D, A, LDA, ISEED, WORK, INFO )
102*
103* -- LAPACK auxiliary routine --
104* -- LAPACK is a software package provided by Univ. of Tennessee, --
105* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
106*
107* .. Scalar Arguments ..
108 INTEGER INFO, K, LDA, N
109* ..
110* .. Array Arguments ..
111 INTEGER ISEED( 4 )
112 DOUBLE PRECISION D( * )
113 COMPLEX*16 A( LDA, * ), WORK( * )
114* ..
115*
116* =====================================================================
117*
118* .. Parameters ..
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* .. Local Scalars ..
125 INTEGER I, J
126 DOUBLE PRECISION WN
127 COMPLEX*16 ALPHA, TAU, WA, WB
128* ..
129* .. External Subroutines ..
130 EXTERNAL xerbla, zaxpy, zgemv, zgerc, zhemv, zher2,
131 $ zlarnv, zscal
132* ..
133* .. External Functions ..
134 DOUBLE PRECISION DZNRM2
135 COMPLEX*16 ZDOTC
136 EXTERNAL dznrm2, zdotc
137* ..
138* .. Intrinsic Functions ..
139 INTRINSIC abs, dble, dconjg, max
140* ..
141* .. Executable Statements ..
142*
143* Test the input arguments
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* initialize lower triangle of A to diagonal matrix
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* Generate lower triangle of hermitian matrix
170*
171 DO 40 i = n - 1, 1, -1
172*
173* generate random reflection
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* apply random reflection to A(i:n,i:n) from the left
188* and the right
189*
190* compute y := tau * A * u
191*
192 CALL zhemv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
193 $ work( n+1 ), 1 )
194*
195* compute v := y - 1/2 * tau * ( y, u ) * u
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* apply the transformation as a rank-2 update to A(i:n,i:n)
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* Reduce number of subdiagonals to K
207*
208 DO 60 i = 1, n - 1 - k
209*
210* generate reflection to annihilate A(k+i+1:n,i)
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* apply reflection to A(k+i:n,i+1:k+i-1) from the left
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* apply reflection to A(k+i:n,k+i:n) from the left and the right
231*
232* compute y := tau * A * u
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* compute v := y - 1/2 * tau * ( y, u ) * u
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* apply hermitian rank-2 update to A(k+i:n,k+i:n)
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* Store full hermitian matrix
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* End of ZLAGHE
263*
264 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:88
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:78
subroutine zhemv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZHEMV
Definition: zhemv.f:154
subroutine zgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERC
Definition: zgerc.f:130
subroutine zher2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZHER2
Definition: zher2.f:150
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
subroutine zlaghe(N, K, D, A, LDA, ISEED, WORK, INFO)
ZLAGHE
Definition: zlaghe.f:102
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: zlarnv.f:99