 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dlagsy()

 subroutine dlagsy ( integer N, integer K, double precision, dimension( * ) D, double precision, dimension( lda, * ) A, integer LDA, integer, dimension( 4 ) ISEED, double precision, dimension( * ) WORK, integer INFO )

DLAGSY

Purpose:
``` DLAGSY generates a real symmetric matrix A, by pre- and post-
multiplying a real diagonal matrix D with a random orthogonal matrix:
A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (2*N)` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```

Definition at line 100 of file dlagsy.f.

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 *
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:97
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
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 dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:147
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: