SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlaghe()

subroutine zlaghe ( integer  n,
integer  k,
double precision, dimension( * )  d,
complex*16, dimension( lda, * )  a,
integer  lda,
integer, dimension( 4 )  iseed,
complex*16, dimension( * )  work,
integer  info 
)

Definition at line 1 of file zlaghe.f.

2*
3* -- LAPACK auxiliary test routine (version 3.1) --
4* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5* November 2006
6*
7* .. Scalar Arguments ..
8 INTEGER INFO, K, LDA, N
9* ..
10* .. Array Arguments ..
11 INTEGER ISEED( 4 )
12 DOUBLE PRECISION D( * )
13 COMPLEX*16 A( LDA, * ), WORK( * )
14* ..
15*
16* Purpose
17* =======
18*
19* ZLAGHE generates a complex hermitian matrix A, by pre- and post-
20* multiplying a real diagonal matrix D with a random unitary matrix:
21* A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
22* unitary transformations.
23*
24* Arguments
25* =========
26*
27* N (input) INTEGER
28* The order of the matrix A. N >= 0.
29*
30* K (input) INTEGER
31* The number of nonzero subdiagonals within the band of A.
32* 0 <= K <= N-1.
33*
34* D (input) DOUBLE PRECISION array, dimension (N)
35* The diagonal elements of the diagonal matrix D.
36*
37* A (output) COMPLEX*16 array, dimension (LDA,N)
38* The generated n by n hermitian matrix A (the full matrix is
39* stored).
40*
41* LDA (input) INTEGER
42* The leading dimension of the array A. LDA >= N.
43*
44* ISEED (input/output) INTEGER array, dimension (4)
45* On entry, the seed of the random number generator; the array
46* elements must be between 0 and 4095, and ISEED(4) must be
47* odd.
48* On exit, the seed is updated.
49*
50* WORK (workspace) COMPLEX*16 array, dimension (2*N)
51*
52* INFO (output) INTEGER
53* = 0: successful exit
54* < 0: if INFO = -i, the i-th argument had an illegal value
55*
56* =====================================================================
57*
58* .. Parameters ..
59 COMPLEX*16 ZERO, ONE, HALF
60 parameter( zero = ( 0.0d+0, 0.0d+0 ),
61 $ one = ( 1.0d+0, 0.0d+0 ),
62 $ half = ( 0.5d+0, 0.0d+0 ) )
63* ..
64* .. Local Scalars ..
65 INTEGER I, J
66 DOUBLE PRECISION WN
67 COMPLEX*16 ALPHA, TAU, WA, WB, DOTC
68* ..
69* .. External Subroutines ..
70 EXTERNAL xerbla, zaxpy, zgemv, zgerc, zhemv, zher2,
71 $ zlarnv, zscal, zzdotc
72* ..
73* .. External Functions ..
74 DOUBLE PRECISION DZNRM2
75 EXTERNAL dznrm2
76* ..
77* .. Intrinsic Functions ..
78 INTRINSIC abs, dble, dconjg, max
79* ..
80* .. Executable Statements ..
81*
82* Test the input arguments
83*
84 info = 0
85 IF( n.LT.0 ) THEN
86 info = -1
87 ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
88 info = -2
89 ELSE IF( lda.LT.max( 1, n ) ) THEN
90 info = -5
91 END IF
92 IF( info.LT.0 ) THEN
93 CALL xerbla( 'ZLAGHE', -info )
94 RETURN
95 END IF
96*
97* initialize lower triangle of A to diagonal matrix
98*
99 DO 20 j = 1, n
100 DO 10 i = j + 1, n
101 a( i, j ) = zero
102 10 CONTINUE
103 20 CONTINUE
104 DO 30 i = 1, n
105 a( i, i ) = d( i )
106 30 CONTINUE
107*
108* Generate lower triangle of hermitian matrix
109*
110 DO 40 i = n - 1, 1, -1
111*
112* generate random reflection
113*
114 CALL zlarnv( 3, iseed, n-i+1, work )
115 wn = dznrm2( n-i+1, work, 1 )
116 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
117 IF( wn.EQ.zero ) THEN
118 tau = zero
119 ELSE
120 wb = work( 1 ) + wa
121 CALL zscal( n-i, one / wb, work( 2 ), 1 )
122 work( 1 ) = one
123 tau = dble( wb / wa )
124 END IF
125*
126* apply random reflection to A(i:n,i:n) from the left
127* and the right
128*
129* compute y := tau * A * u
130*
131 CALL zhemv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
132 $ work( n+1 ), 1 )
133*
134* compute v := y - 1/2 * tau * ( y, u ) * u
135*
136 CALL zzdotc( n-i+1, dotc, work( n+1 ), 1, work, 1 )
137 alpha = -half*tau*dotc
138 CALL zaxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
139*
140* apply the transformation as a rank-2 update to A(i:n,i:n)
141*
142 CALL zher2( 'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
143 $ a( i, i ), lda )
144 40 CONTINUE
145*
146* Reduce number of subdiagonals to K
147*
148 DO 60 i = 1, n - 1 - k
149*
150* generate reflection to annihilate A(k+i+1:n,i)
151*
152 wn = dznrm2( n-k-i+1, a( k+i, i ), 1 )
153 wa = ( wn / abs( a( k+i, i ) ) )*a( k+i, i )
154 IF( wn.EQ.zero ) THEN
155 tau = zero
156 ELSE
157 wb = a( k+i, i ) + wa
158 CALL zscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
159 a( k+i, i ) = one
160 tau = dble( wb / wa )
161 END IF
162*
163* apply reflection to A(k+i:n,i+1:k+i-1) from the left
164*
165 CALL zgemv( 'Conjugate transpose', n-k-i+1, k-1, one,
166 $ a( k+i, i+1 ), lda, a( k+i, i ), 1, zero, work, 1 )
167 CALL zgerc( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
168 $ a( k+i, i+1 ), lda )
169*
170* apply reflection to A(k+i:n,k+i:n) from the left and the right
171*
172* compute y := tau * A * u
173*
174 CALL zhemv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
175 $ a( k+i, i ), 1, zero, work, 1 )
176*
177* compute v := y - 1/2 * tau * ( y, u ) * u
178*
179 CALL zzdotc( n-k-i+1, dotc, work, 1, a( k+i, i ), 1 )
180 alpha = -half*tau*dotc
181 CALL zaxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
182*
183* apply hermitian rank-2 update to A(k+i:n,k+i:n)
184*
185 CALL zher2( 'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
186 $ a( k+i, k+i ), lda )
187*
188 a( k+i, i ) = -wa
189 DO 50 j = k + i + 1, n
190 a( j, i ) = zero
191 50 CONTINUE
192 60 CONTINUE
193*
194* Store full hermitian matrix
195*
196 DO 80 j = 1, n
197 DO 70 i = j + 1, n
198 a( j, i ) = dconjg( a( i, j ) )
199 70 CONTINUE
200 80 CONTINUE
201 RETURN
202*
203* End of ZLAGHE
204*
max
#define max(A, B)
Definition pcgemr.c:180
zlarnv
subroutine zlarnv(idist, iseed, n, x)
Definition zlarnv.f:2
zzdotc
subroutine zzdotc(n, dotc, x, incx, y, incy)
Definition zzdotc.f:2
Here is the call graph for this function:
Here is the caller graph for this function: