ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
dlagsy.f
Go to the documentation of this file.
1  SUBROUTINE dlagsy( N, K, D, A, LDA, ISEED, WORK, INFO )
2 *
3 * -- LAPACK auxiliary test routine (version 3.1)
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8  INTEGER INFO, K, LDA, N
9 * ..
10 * .. Array Arguments ..
11  INTEGER ISEED( 4 )
12  DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DLAGSY generates a real symmetric matrix A, by pre- and post-
19 * multiplying a real diagonal matrix D with a random orthogonal matrix:
20 * A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
21 * orthogonal transformations.
22 *
23 * Arguments
24 * =========
25 *
26 * N (input) INTEGER
27 * The order of the matrix A. N >= 0.
28 *
29 * K (input) INTEGER
30 * The number of nonzero subdiagonals within the band of A.
31 * 0 <= K <= N-1.
32 *
33 * D (input) DOUBLE PRECISION array, dimension (N)
34 * The diagonal elements of the diagonal matrix D.
35 *
36 * A (output) DOUBLE PRECISION array, dimension (LDA,N)
37 * The generated n by n symmetric matrix A (the full matrix is
38 * stored).
39 *
40 * LDA (input) INTEGER
41 * The leading dimension of the array A. LDA >= N.
42 *
43 * ISEED (input/output) INTEGER array, dimension (4)
44 * On entry, the seed of the random number generator; the array
45 * elements must be between 0 and 4095, and ISEED(4) must be
46 * odd.
47 * On exit, the seed is updated.
48 *
49 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
50 *
51 * INFO (output) INTEGER
52 * = 0: successful exit
53 * < 0: if INFO = -i, the i-th argument had an illegal value
54 *
55 * =====================================================================
56 *
57 * .. Parameters ..
58  DOUBLE PRECISION ZERO, ONE, HALF
59  parameter( zero = 0.0d+0, one = 1.0d+0, half = 0.5d+0 )
60 * ..
61 * .. Local Scalars ..
62  INTEGER I, J
63  DOUBLE PRECISION ALPHA, TAU, WA, WB, WN
64 * ..
65 * .. External Subroutines ..
66  EXTERNAL daxpy, dgemv, dger, dlarnv, dscal, dsymv,
67  $ dsyr2, xerbla
68 * ..
69 * .. External Functions ..
70  DOUBLE PRECISION DDOT, DNRM2
71  EXTERNAL ddot, dnrm2
72 * ..
73 * .. Intrinsic Functions ..
74  INTRINSIC max, sign
75 * ..
76 * .. Executable Statements ..
77 *
78 * Test the input arguments
79 *
80  info = 0
81  IF( n.LT.0 ) THEN
82  info = -1
83  ELSE IF( k.LT.0 .OR. k.GT.n-1 ) THEN
84  info = -2
85  ELSE IF( lda.LT.max( 1, n ) ) THEN
86  info = -5
87  END IF
88  IF( info.LT.0 ) THEN
89  CALL xerbla( 'DLAGSY', -info )
90  RETURN
91  END IF
92 *
93 * initialize lower triangle of A to diagonal matrix
94 *
95  DO 20 j = 1, n
96  DO 10 i = j + 1, n
97  a( i, j ) = zero
98  10 CONTINUE
99  20 CONTINUE
100  DO 30 i = 1, n
101  a( i, i ) = d( i )
102  30 CONTINUE
103 *
104 * Generate lower triangle of symmetric matrix
105 *
106  DO 40 i = n - 1, 1, -1
107 *
108 * generate random reflection
109 *
110  CALL dlarnv( 3, iseed, n-i+1, work )
111  wn = dnrm2( n-i+1, work, 1 )
112  wa = sign( wn, work( 1 ) )
113  IF( wn.EQ.zero ) THEN
114  tau = zero
115  ELSE
116  wb = work( 1 ) + wa
117  CALL dscal( n-i, one / wb, work( 2 ), 1 )
118  work( 1 ) = one
119  tau = wb / wa
120  END IF
121 *
122 * apply random reflection to A(i:n,i:n) from the left
123 * and the right
124 *
125 * compute y := tau * A * u
126 *
127  CALL dsymv( 'Lower', n-i+1, tau, a( i, i ), lda, work, 1, zero,
128  $ work( n+1 ), 1 )
129 *
130 * compute v := y - 1/2 * tau * ( y, u ) * u
131 *
132  alpha = -half*tau*ddot( n-i+1, work( n+1 ), 1, work, 1 )
133  CALL daxpy( n-i+1, alpha, work, 1, work( n+1 ), 1 )
134 *
135 * apply the transformation as a rank-2 update to A(i:n,i:n)
136 *
137  CALL dsyr2( 'Lower', n-i+1, -one, work, 1, work( n+1 ), 1,
138  $ a( i, i ), lda )
139  40 CONTINUE
140 *
141 * Reduce number of subdiagonals to K
142 *
143  DO 60 i = 1, n - 1 - k
144 *
145 * generate reflection to annihilate A(k+i+1:n,i)
146 *
147  wn = dnrm2( n-k-i+1, a( k+i, i ), 1 )
148  wa = sign( wn, a( k+i, i ) )
149  IF( wn.EQ.zero ) THEN
150  tau = zero
151  ELSE
152  wb = a( k+i, i ) + wa
153  CALL dscal( n-k-i, one / wb, a( k+i+1, i ), 1 )
154  a( k+i, i ) = one
155  tau = wb / wa
156  END IF
157 *
158 * apply reflection to A(k+i:n,i+1:k+i-1) from the left
159 *
160  CALL dgemv( 'Transpose', n-k-i+1, k-1, one, a( k+i, i+1 ), lda,
161  $ a( k+i, i ), 1, zero, work, 1 )
162  CALL dger( n-k-i+1, k-1, -tau, a( k+i, i ), 1, work, 1,
163  $ a( k+i, i+1 ), lda )
164 *
165 * apply reflection to A(k+i:n,k+i:n) from the left and the right
166 *
167 * compute y := tau * A * u
168 *
169  CALL dsymv( 'Lower', n-k-i+1, tau, a( k+i, k+i ), lda,
170  $ a( k+i, i ), 1, zero, work, 1 )
171 *
172 * compute v := y - 1/2 * tau * ( y, u ) * u
173 *
174  alpha = -half*tau*ddot( n-k-i+1, work, 1, a( k+i, i ), 1 )
175  CALL daxpy( n-k-i+1, alpha, a( k+i, i ), 1, work, 1 )
176 *
177 * apply symmetric rank-2 update to A(k+i:n,k+i:n)
178 *
179  CALL dsyr2( 'Lower', n-k-i+1, -one, a( k+i, i ), 1, work, 1,
180  $ a( k+i, k+i ), lda )
181 *
182  a( k+i, i ) = -wa
183  DO 50 j = k + i + 1, n
184  a( j, i ) = zero
185  50 CONTINUE
186  60 CONTINUE
187 *
188 * Store full symmetric matrix
189 *
190  DO 80 j = 1, n
191  DO 70 i = j + 1, n
192  a( j, i ) = a( i, j )
193  70 CONTINUE
194  80 CONTINUE
195  RETURN
196 *
197 * End of DLAGSY
198 *
199  END
max
#define max(A, B)
Definition: pcgemr.c:180
dlagsy
subroutine dlagsy(N, K, D, A, LDA, ISEED, WORK, INFO)
Definition: dlagsy.f:2