LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 104 of file zlagsy.f.

104 *
105 * -- LAPACK auxiliary routine (version 3.4.0) --
106 * -- LAPACK is a software package provided by Univ. of Tennessee, --
107 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108 * November 2011
109 *
110 * .. Scalar Arguments ..
111  INTEGER info, k, lda, n
112 * ..
113 * .. Array Arguments ..
114  INTEGER iseed( 4 )
115  DOUBLE PRECISION d( * )
116  COMPLEX*16 a( lda, * ), work( * )
117 * ..
118 *
119 * =====================================================================
120 *
121 * .. Parameters ..
122  COMPLEX*16 zero, one, half
123  parameter ( zero = ( 0.0d+0, 0.0d+0 ),
124  $ one = ( 1.0d+0, 0.0d+0 ),
125  $ half = ( 0.5d+0, 0.0d+0 ) )
126 * ..
127 * .. Local Scalars ..
128  INTEGER i, ii, j, jj
129  DOUBLE PRECISION wn
130  COMPLEX*16 alpha, tau, wa, wb
131 * ..
132 * .. External Subroutines ..
133  EXTERNAL xerbla, zaxpy, zgemv, zgerc, zlacgv, zlarnv,
134  $ zscal, zsymv
135 * ..
136 * .. External Functions ..
137  DOUBLE PRECISION dznrm2
138  COMPLEX*16 zdotc
139  EXTERNAL dznrm2, zdotc
140 * ..
141 * .. Intrinsic Functions ..
142  INTRINSIC abs, dble, max
143 * ..
144 * .. Executable Statements ..
145 *
146 * Test the input arguments
147 *
148  info = 0
149  IF( n.LT.0 ) THEN
150  info = -1
151  ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
152  info = -2
153  ELSE IF( lda.LT.max( 1, n ) ) THEN
154  info = -5
155  END IF
156  IF( info.LT.0 ) THEN
157  CALL xerbla( 'ZLAGSY', -info )
158  RETURN
159  END IF
160 *
161 * initialize lower triangle of A to diagonal matrix
162 *
163  DO 20 j = 1, n
164  DO 10 i = j + 1, n
165  a( i, j ) = zero
166  10 CONTINUE
167  20 CONTINUE
168  DO 30 i = 1, n
169  a( i, i ) = d( i )
170  30 CONTINUE
171 *
172 * Generate lower triangle of symmetric matrix
173 *
174  DO 60 i = n - 1, 1, -1
175 *
176 * generate random reflection
177 *
178  CALL zlarnv( 3, iseed, n-i+1, work )
179  wn = dznrm2( n-i+1, work, 1 )
180  wa = ( wn / abs( work( 1 ) ) )*work( 1 )
181  IF( wn.EQ.zero ) THEN
182  tau = zero
183  ELSE
184  wb = work( 1 ) + wa
185  CALL zscal( n-i, one / wb, work( 2 ), 1 )
186  work( 1 ) = one
187  tau = dble( wb / wa )
188  END IF
189 *
190 * apply random reflection to A(i:n,i:n) from the left
191 * and the right
192 *
193 * compute y := tau * A * conjg(u)
194 *
195  CALL zlacgv( n-i+1, work, 1 )
196  CALL zsymv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
197  $ work( n+1 ), 1 )
198  CALL zlacgv( n-i+1, work, 1 )
199 *
200 * compute v := y - 1/2 * tau * ( u, y ) * u
201 *
202  alpha = -half*tau*zdotc( n-i+1, work, 1, work( n+1 ), 1 )
203  CALL zaxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
204 *
205 * apply the transformation as a rank-2 update to A(i:n,i:n)
206 *
207 * CALL ZSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1,
208 * $ A( I, I ), LDA )
209 *
210  DO 50 jj = i, n
211  DO 40 ii = jj, n
212  a( ii, jj ) = a( ii, jj ) -
213  $ work( ii-i+1 )*work( n+jj-i+1 ) -
214  $ work( n+ii-i+1 )*work( jj-i+1 )
215  40 CONTINUE
216  50 CONTINUE
217  60 CONTINUE
218 *
219 * Reduce number of subdiagonals to K
220 *
221  DO 100 i = 1, n - 1 - k
222 *
223 * generate reflection to annihilate A(k+i+1:n,i)
224 *
225  wn = dznrm2( n-k-i+1, a( k+i, i ), 1 )
226  wa = ( wn / abs( a( k+i, i ) ) )*a( k+i, i )
227  IF( wn.EQ.zero ) THEN
228  tau = zero
229  ELSE
230  wb = a( k+i, i ) + wa
231  CALL zscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
232  a( k+i, i ) = one
233  tau = dble( wb / wa )
234  END IF
235 *
236 * apply reflection to A(k+i:n,i+1:k+i-1) from the left
237 *
238  CALL zgemv( 'Conjugate transpose', n-k-i+1, k-1, one,
239  $ a( k+i, i+1 ), lda, a( k+i, i ), 1, zero, work, 1 )
240  CALL zgerc( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
241  $ a( k+i, i+1 ), lda )
242 *
243 * apply reflection to A(k+i:n,k+i:n) from the left and the right
244 *
245 * compute y := tau * A * conjg(u)
246 *
247  CALL zlacgv( n-k-i+1, a( k+i, i ), 1 )
248  CALL zsymv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
249  $ a( k+i, i ), 1, zero, work, 1 )
250  CALL zlacgv( n-k-i+1, a( k+i, i ), 1 )
251 *
252 * compute v := y - 1/2 * tau * ( u, y ) * u
253 *
254  alpha = -half*tau*zdotc( n-k-i+1, a( k+i, i ), 1, work, 1 )
255  CALL zaxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
256 *
257 * apply symmetric rank-2 update to A(k+i:n,k+i:n)
258 *
259 * CALL ZSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1,
260 * $ A( K+I, K+I ), LDA )
261 *
262  DO 80 jj = k + i, n
263  DO 70 ii = jj, n
264  a( ii, jj ) = a( ii, jj ) - a( ii, i )*work( jj-k-i+1 ) -
265  $ work( ii-k-i+1 )*a( jj, i )
266  70 CONTINUE
267  80 CONTINUE
268 *
269  a( k+i, i ) = -wa
270  DO 90 j = k + i + 1, n
271  a( j, i ) = zero
272  90 CONTINUE
273  100 CONTINUE
274 *
275 * Store full symmetric matrix
276 *
277  DO 120 j = 1, n
278  DO 110 i = j + 1, n
279  a( j, i ) = a( i, j )
280  110 CONTINUE
281  120 CONTINUE
282  RETURN
283 *
284 * End of ZLAGSY
285 *
zlarnv
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: zlarnv.f:101
zgemv
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
zgerc
subroutine zgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERC
Definition: zgerc.f:132
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zsymv
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:159
zdotc
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
Definition: zdotc.f:54
dznrm2
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:56
zaxpy
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:53
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54
zlacgv
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76

Here is the call graph for this function:

Here is the caller graph for this function: