LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlagsy.f
Go to the documentation of this file.
1*> \brief \b DLAGSY
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 DLAGSY( 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 A( LDA, * ), D( * ), WORK( * )
19* ..
20*
21*
22*> \par Purpose:
23* =============
24*>
25*> \verbatim
26*>
27*> DLAGSY generates a real symmetric matrix A, by pre- and post-
28*> multiplying a real diagonal matrix D with a random orthogonal matrix:
29*> A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
30*> orthogonal transformations.
31*> \endverbatim
32*
33* Arguments:
34* ==========
35*
36*> \param[in] N
37*> \verbatim
38*> N is INTEGER
39*> The order of the matrix A. N >= 0.
40*> \endverbatim
41*>
42*> \param[in] K
43*> \verbatim
44*> K is INTEGER
45*> The number of nonzero subdiagonals within the band of A.
46*> 0 <= K <= N-1.
47*> \endverbatim
48*>
49*> \param[in] D
50*> \verbatim
51*> D is DOUBLE PRECISION array, dimension (N)
52*> The diagonal elements of the diagonal matrix D.
53*> \endverbatim
54*>
55*> \param[out] A
56*> \verbatim
57*> A is DOUBLE PRECISION array, dimension (LDA,N)
58*> The generated n by n symmetric matrix A (the full matrix is
59*> stored).
60*> \endverbatim
61*>
62*> \param[in] LDA
63*> \verbatim
64*> LDA is INTEGER
65*> The leading dimension of the array A. LDA >= N.
66*> \endverbatim
67*>
68*> \param[in,out] ISEED
69*> \verbatim
70*> ISEED is INTEGER array, dimension (4)
71*> On entry, the seed of the random number generator; the array
72*> elements must be between 0 and 4095, and ISEED(4) must be
73*> odd.
74*> On exit, the seed is updated.
75*> \endverbatim
76*>
77*> \param[out] WORK
78*> \verbatim
79*> WORK is DOUBLE PRECISION array, dimension (2*N)
80*> \endverbatim
81*>
82*> \param[out] INFO
83*> \verbatim
84*> INFO is INTEGER
85*> = 0: successful exit
86*> < 0: if INFO = -i, the i-th argument had an illegal value
87*> \endverbatim
88*
89* Authors:
90* ========
91*
92*> \author Univ. of Tennessee
93*> \author Univ. of California Berkeley
94*> \author Univ. of Colorado Denver
95*> \author NAG Ltd.
96*
97*> \ingroup double_matgen
98*
99* =====================================================================
100 SUBROUTINE dlagsy( N, K, D, A, LDA, ISEED, WORK, INFO )
101*
102* -- LAPACK auxiliary routine --
103* -- LAPACK is a software package provided by Univ. of Tennessee, --
104* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
105*
106* .. Scalar Arguments ..
107 INTEGER INFO, K, LDA, N
108* ..
109* .. Array Arguments ..
110 INTEGER ISEED( 4 )
111 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
112* ..
113*
114* =====================================================================
115*
116* .. Parameters ..
117 DOUBLE PRECISION ZERO, ONE, HALF
118 parameter( zero = 0.0d+0, one = 1.0d+0, half = 0.5d+0 )
119* ..
120* .. Local Scalars ..
121 INTEGER I, J
122 DOUBLE PRECISION ALPHA, TAU, WA, WB, WN
123* ..
124* .. External Subroutines ..
125 EXTERNAL daxpy, dgemv, dger, dlarnv, dscal, dsymv,
126 $ dsyr2, xerbla
127* ..
128* .. External Functions ..
129 DOUBLE PRECISION DDOT, DNRM2
130 EXTERNAL ddot, dnrm2
131* ..
132* .. Intrinsic Functions ..
133 INTRINSIC max, sign
134* ..
135* .. Executable Statements ..
136*
137* Test the input arguments
138*
139 info = 0
140 IF( n.LT.0 ) THEN
141 info = -1
142 ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
143 info = -2
144 ELSE IF( lda.LT.max( 1, n ) ) THEN
145 info = -5
146 END IF
147 IF( info.LT.0 ) THEN
148 CALL xerbla( 'DLAGSY', -info )
149 RETURN
150 END IF
151*
152* initialize lower triangle of A to diagonal matrix
153*
154 DO 20 j = 1, n
155 DO 10 i = j + 1, n
156 a( i, j ) = zero
157 10 CONTINUE
158 20 CONTINUE
159 DO 30 i = 1, n
160 a( i, i ) = d( i )
161 30 CONTINUE
162*
163* Generate lower triangle of symmetric matrix
164*
165 DO 40 i = n - 1, 1, -1
166*
167* generate random reflection
168*
169 CALL dlarnv( 3, iseed, n-i+1, work )
170 wn = dnrm2( n-i+1, work, 1 )
171 wa = sign( wn, work( 1 ) )
172 IF( wn.EQ.zero ) THEN
173 tau = zero
174 ELSE
175 wb = work( 1 ) + wa
176 CALL dscal( n-i, one / wb, work( 2 ), 1 )
177 work( 1 ) = one
178 tau = wb / wa
179 END IF
180*
181* apply random reflection to A(i:n,i:n) from the left
182* and the right
183*
184* compute y := tau * A * u
185*
186 CALL dsymv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
187 $ work( n+1 ), 1 )
188*
189* compute v := y - 1/2 * tau * ( y, u ) * u
190*
191 alpha = -half*tau*ddot( n-i+1, work( n+1 ), 1, work, 1 )
192 CALL daxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
193*
194* apply the transformation as a rank-2 update to A(i:n,i:n)
195*
196 CALL dsyr2( 'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
197 $ a( i, i ), lda )
198 40 CONTINUE
199*
200* Reduce number of subdiagonals to K
201*
202 DO 60 i = 1, n - 1 - k
203*
204* generate reflection to annihilate A(k+i+1:n,i)
205*
206 wn = dnrm2( n-k-i+1, a( k+i, i ), 1 )
207 wa = sign( wn, a( k+i, i ) )
208 IF( wn.EQ.zero ) THEN
209 tau = zero
210 ELSE
211 wb = a( k+i, i ) + wa
212 CALL dscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
213 a( k+i, i ) = one
214 tau = wb / wa
215 END IF
216*
217* apply reflection to A(k+i:n,i+1:k+i-1) from the left
218*
219 CALL dgemv( 'Transpose', n-k-i+1, k-1, one, a( k+i, i+1 ), lda,
220 $ a( k+i, i ), 1, zero, work, 1 )
221 CALL dger( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
222 $ a( k+i, i+1 ), lda )
223*
224* apply reflection to A(k+i:n,k+i:n) from the left and the right
225*
226* compute y := tau * A * u
227*
228 CALL dsymv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
229 $ a( k+i, i ), 1, zero, work, 1 )
230*
231* compute v := y - 1/2 * tau * ( y, u ) * u
232*
233 alpha = -half*tau*ddot( n-k-i+1, work, 1, a( k+i, i ), 1 )
234 CALL daxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
235*
236* apply symmetric rank-2 update to A(k+i:n,k+i:n)
237*
238 CALL dsyr2( 'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
239 $ a( k+i, k+i ), lda )
240*
241 a( k+i, i ) = -wa
242 DO 50 j = k + i + 1, n
243 a( j, i ) = zero
244 50 CONTINUE
245 60 CONTINUE
246*
247* Store full symmetric matrix
248*
249 DO 80 j = 1, n
250 DO 70 i = j + 1, n
251 a( j, i ) = a( i, j )
252 70 CONTINUE
253 80 CONTINUE
254 RETURN
255*
256* End of DLAGSY
257*
258 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlagsy(n, k, d, a, lda, iseed, work, info)
DLAGSY
Definition dlagsy.f:101
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
subroutine dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DSYMV
Definition dsymv.f:152
subroutine dsyr2(uplo, n, alpha, x, incx, y, incy, a, lda)
DSYR2
Definition dsyr2.f:147
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79