LAPACK  3.8.0
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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016

Definition at line 103 of file dlagsy.f.

103 *
104 * -- LAPACK auxiliary routine (version 3.7.0) --
105 * -- LAPACK is a software package provided by Univ. of Tennessee, --
106 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
107 * December 2016
108 *
109 * .. Scalar Arguments ..
110  INTEGER info, k, lda, n
111 * ..
112 * .. Array Arguments ..
113  INTEGER iseed( 4 )
114  DOUBLE PRECISION a( lda, * ), d( * ), work( * )
115 * ..
116 *
117 * =====================================================================
118 *
119 * .. Parameters ..
120  DOUBLE PRECISION zero, one, half
121  parameter( zero = 0.0d+0, one = 1.0d+0, half = 0.5d+0 )
122 * ..
123 * .. Local Scalars ..
124  INTEGER i, j
125  DOUBLE PRECISION alpha, tau, wa, wb, wn
126 * ..
127 * .. External Subroutines ..
128  EXTERNAL daxpy, dgemv, dger, dlarnv, dscal, dsymv,
129  $ dsyr2, xerbla
130 * ..
131 * .. External Functions ..
132  DOUBLE PRECISION ddot, dnrm2
133  EXTERNAL ddot, dnrm2
134 * ..
135 * .. Intrinsic Functions ..
136  INTRINSIC max, sign
137 * ..
138 * .. Executable Statements ..
139 *
140 * Test the input arguments
141 *
142  info = 0
143  IF( n.LT.0 ) THEN
144  info = -1
145  ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
146  info = -2
147  ELSE IF( lda.LT.max( 1, n ) ) THEN
148  info = -5
149  END IF
150  IF( info.LT.0 ) THEN
151  CALL xerbla( 'DLAGSY', -info )
152  RETURN
153  END IF
154 *
155 * initialize lower triangle of A to diagonal matrix
156 *
157  DO 20 j = 1, n
158  DO 10 i = j + 1, n
159  a( i, j ) = zero
160  10 CONTINUE
161  20 CONTINUE
162  DO 30 i = 1, n
163  a( i, i ) = d( i )
164  30 CONTINUE
165 *
166 * Generate lower triangle of symmetric matrix
167 *
168  DO 40 i = n - 1, 1, -1
169 *
170 * generate random reflection
171 *
172  CALL dlarnv( 3, iseed, n-i+1, work )
173  wn = dnrm2( n-i+1, work, 1 )
174  wa = sign( wn, work( 1 ) )
175  IF( wn.EQ.zero ) THEN
176  tau = zero
177  ELSE
178  wb = work( 1 ) + wa
179  CALL dscal( n-i, one / wb, work( 2 ), 1 )
180  work( 1 ) = one
181  tau = wb / wa
182  END IF
183 *
184 * apply random reflection to A(i:n,i:n) from the left
185 * and the right
186 *
187 * compute y := tau * A * u
188 *
189  CALL dsymv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
190  $ work( n+1 ), 1 )
191 *
192 * compute v := y - 1/2 * tau * ( y, u ) * u
193 *
194  alpha = -half*tau*ddot( n-i+1, work( n+1 ), 1, work, 1 )
195  CALL daxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
196 *
197 * apply the transformation as a rank-2 update to A(i:n,i:n)
198 *
199  CALL dsyr2( 'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
200  $ a( i, i ), lda )
201  40 CONTINUE
202 *
203 * Reduce number of subdiagonals to K
204 *
205  DO 60 i = 1, n - 1 - k
206 *
207 * generate reflection to annihilate A(k+i+1:n,i)
208 *
209  wn = dnrm2( n-k-i+1, a( k+i, i ), 1 )
210  wa = sign( wn, a( k+i, i ) )
211  IF( wn.EQ.zero ) THEN
212  tau = zero
213  ELSE
214  wb = a( k+i, i ) + wa
215  CALL dscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
216  a( k+i, i ) = one
217  tau = wb / wa
218  END IF
219 *
220 * apply reflection to A(k+i:n,i+1:k+i-1) from the left
221 *
222  CALL dgemv( 'Transpose', n-k-i+1, k-1, one, a( k+i, i+1 ), lda,
223  $ a( k+i, i ), 1, zero, work, 1 )
224  CALL dger( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
225  $ a( k+i, i+1 ), lda )
226 *
227 * apply reflection to A(k+i:n,k+i:n) from the left and the right
228 *
229 * compute y := tau * A * u
230 *
231  CALL dsymv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
232  $ a( k+i, i ), 1, zero, work, 1 )
233 *
234 * compute v := y - 1/2 * tau * ( y, u ) * u
235 *
236  alpha = -half*tau*ddot( n-k-i+1, work, 1, a( k+i, i ), 1 )
237  CALL daxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
238 *
239 * apply symmetric rank-2 update to A(k+i:n,k+i:n)
240 *
241  CALL dsyr2( 'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
242  $ a( k+i, k+i ), lda )
243 *
244  a( k+i, i ) = -wa
245  DO 50 j = k + i + 1, n
246  a( j, i ) = zero
247  50 CONTINUE
248  60 CONTINUE
249 *
250 * Store full symmetric matrix
251 *
252  DO 80 j = 1, n
253  DO 70 i = j + 1, n
254  a( j, i ) = a( i, j )
255  70 CONTINUE
256  80 CONTINUE
257  RETURN
258 *
259 * End of DLAGSY
260 *
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:84
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
subroutine dsyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DSYR2
Definition: dsyr2.f:149
subroutine dsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DSYMV
Definition: dsymv.f:154
Here is the call graph for this function:
Here is the caller graph for this function: