1 SUBROUTINE SLAGSY( 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 REAL A( LDA, * ), D( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * SLAGSY 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) REAL array, dimension (N)
34 * The diagonal elements of the diagonal matrix D.
35 *
36 * A (output) REAL 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) REAL 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 REAL ZERO, ONE, HALF
59 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, HALF = 0.5E+0 )
60 * ..
61 * .. Local Scalars ..
62 INTEGER I, J
63 REAL ALPHA, TAU, WA, WB, WN
64 * ..
65 * .. External Subroutines ..
66 EXTERNAL SAXPY, SGEMV, SGER, SLARNV, SSCAL, SSYMV,
67 $ SSYR2, XERBLA
68 * ..
69 * .. External Functions ..
70 REAL SDOT, SNRM2
71 EXTERNAL SDOT, SNRM2
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( 'SLAGSY', -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 SLARNV( 3, ISEED, N-I+1, WORK )
111 WN = SNRM2( 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 SSCAL( 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 SSYMV( '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*SDOT( N-I+1, WORK( N+1 ), 1, WORK, 1 )
133 CALL SAXPY( 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 SSYR2( '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 = SNRM2( 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 SSCAL( 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 SGEMV( '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 SGER( 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 SSYMV( '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*SDOT( N-K-I+1, WORK, 1, A( K+I, I ), 1 )
175 CALL SAXPY( 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 SSYR2( '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 SLAGSY
198 *
199 END
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 REAL A( LDA, * ), D( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * SLAGSY 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) REAL array, dimension (N)
34 * The diagonal elements of the diagonal matrix D.
35 *
36 * A (output) REAL 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) REAL 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 REAL ZERO, ONE, HALF
59 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, HALF = 0.5E+0 )
60 * ..
61 * .. Local Scalars ..
62 INTEGER I, J
63 REAL ALPHA, TAU, WA, WB, WN
64 * ..
65 * .. External Subroutines ..
66 EXTERNAL SAXPY, SGEMV, SGER, SLARNV, SSCAL, SSYMV,
67 $ SSYR2, XERBLA
68 * ..
69 * .. External Functions ..
70 REAL SDOT, SNRM2
71 EXTERNAL SDOT, SNRM2
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( 'SLAGSY', -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 SLARNV( 3, ISEED, N-I+1, WORK )
111 WN = SNRM2( 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 SSCAL( 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 SSYMV( '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*SDOT( N-I+1, WORK( N+1 ), 1, WORK, 1 )
133 CALL SAXPY( 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 SSYR2( '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 = SNRM2( 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 SSCAL( 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 SGEMV( '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 SGER( 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 SSYMV( '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*SDOT( N-K-I+1, WORK, 1, A( K+I, I ), 1 )
175 CALL SAXPY( 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 SSYR2( '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 SLAGSY
198 *
199 END