LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlagsy()

subroutine zlagsy ( 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 
)

ZLAGSY

Purpose:
 ZLAGSY generates a complex symmetric matrix A, by pre- and post-
 multiplying a real diagonal matrix D with a random unitary matrix:
 A = U*D*U**T. The semi-bandwidth may then be reduced to k by
 additional unitary transformations.
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]K
          K is INTEGER
          The number of nonzero subdiagonals within the band of A.
          0 <= K <= N-1.
[in]D
          D is DOUBLE PRECISION array, dimension (N)
          The diagonal elements of the diagonal matrix D.
[out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          The generated n by n symmetric matrix A (the full matrix is
          stored).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= N.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry, the seed of the random number generator; the array
          elements must be between 0 and 4095, and ISEED(4) must be
          odd.
          On exit, the seed is updated.
[out]WORK
          WORK is COMPLEX*16 array, dimension (2*N)
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 101 of file zlagsy.f.

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, II, J, JJ
126 DOUBLE PRECISION WN
127 COMPLEX*16 ALPHA, TAU, WA, WB
128* ..
129* .. External Subroutines ..
130 EXTERNAL xerbla, zaxpy, zgemv, zgerc, zlacgv, zlarnv,
131 $ zscal, zsymv
132* ..
133* .. External Functions ..
134 DOUBLE PRECISION DZNRM2
135 COMPLEX*16 ZDOTC
136 EXTERNAL dznrm2, zdotc
137* ..
138* .. Intrinsic Functions ..
139 INTRINSIC abs, dble, 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( 'ZLAGSY', -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 symmetric matrix
170*
171 DO 60 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 * conjg(u)
191*
192 CALL zlacgv( n-i+1, work, 1 )
193 CALL zsymv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
194 $ work( n+1 ), 1 )
195 CALL zlacgv( n-i+1, work, 1 )
196*
197* compute v := y - 1/2 * tau * ( u, y ) * u
198*
199 alpha = -half*tau*zdotc( n-i+1, work, 1, work( n+1 ), 1 )
200 CALL zaxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
201*
202* apply the transformation as a rank-2 update to A(i:n,i:n)
203*
204* CALL ZSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1,
205* $ A( I, I ), LDA )
206*
207 DO 50 jj = i, n
208 DO 40 ii = jj, n
209 a( ii, jj ) = a( ii, jj ) -
210 $ work( ii-i+1 )*work( n+jj-i+1 ) -
211 $ work( n+ii-i+1 )*work( jj-i+1 )
212 40 CONTINUE
213 50 CONTINUE
214 60 CONTINUE
215*
216* Reduce number of subdiagonals to K
217*
218 DO 100 i = 1, n - 1 - k
219*
220* generate reflection to annihilate A(k+i+1:n,i)
221*
222 wn = dznrm2( n-k-i+1, a( k+i, i ), 1 )
223 wa = ( wn / abs( a( k+i, i ) ) )*a( k+i, i )
224 IF( wn.EQ.zero ) THEN
225 tau = zero
226 ELSE
227 wb = a( k+i, i ) + wa
228 CALL zscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
229 a( k+i, i ) = one
230 tau = dble( wb / wa )
231 END IF
232*
233* apply reflection to A(k+i:n,i+1:k+i-1) from the left
234*
235 CALL zgemv( 'Conjugate transpose', n-k-i+1, k-1, one,
236 $ a( k+i, i+1 ), lda, a( k+i, i ), 1, zero, work, 1 )
237 CALL zgerc( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
238 $ a( k+i, i+1 ), lda )
239*
240* apply reflection to A(k+i:n,k+i:n) from the left and the right
241*
242* compute y := tau * A * conjg(u)
243*
244 CALL zlacgv( n-k-i+1, a( k+i, i ), 1 )
245 CALL zsymv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
246 $ a( k+i, i ), 1, zero, work, 1 )
247 CALL zlacgv( n-k-i+1, a( k+i, i ), 1 )
248*
249* compute v := y - 1/2 * tau * ( u, y ) * u
250*
251 alpha = -half*tau*zdotc( n-k-i+1, a( k+i, i ), 1, work, 1 )
252 CALL zaxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
253*
254* apply symmetric rank-2 update to A(k+i:n,k+i:n)
255*
256* CALL ZSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1,
257* $ A( K+I, K+I ), LDA )
258*
259 DO 80 jj = k + i, n
260 DO 70 ii = jj, n
261 a( ii, jj ) = a( ii, jj ) - a( ii, i )*work( jj-k-i+1 ) -
262 $ work( ii-k-i+1 )*a( jj, i )
263 70 CONTINUE
264 80 CONTINUE
265*
266 a( k+i, i ) = -wa
267 DO 90 j = k + i + 1, n
268 a( j, i ) = zero
269 90 CONTINUE
270 100 CONTINUE
271*
272* Store full symmetric matrix
273*
274 DO 120 j = 1, n
275 DO 110 i = j + 1, n
276 a( j, i ) = a( i, j )
277 110 CONTINUE
278 120 CONTINUE
279 RETURN
280*
281* End of ZLAGSY
282*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
complex *16 function zdotc(n, zx, incx, zy, incy)
ZDOTC
Definition zdotc.f:83
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:160
subroutine zgerc(m, n, alpha, x, incx, y, incy, a, lda)
ZGERC
Definition zgerc.f:130
subroutine zsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
ZSYMV computes a matrix-vector product for a complex symmetric matrix.
Definition zsymv.f:157
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:74
subroutine zlarnv(idist, iseed, n, x)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition zlarnv.f:99
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: