1 SUBROUTINE SSTT22( N, M, KBAND, AD, AE, SD, SE, U, LDU, WORK,
2 $ LDWORK, RESULT )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER KBAND, LDU, LDWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 REAL AD( * ), AE( * ), RESULT( 2 ), SD( * ),
13 $ SE( * ), U( LDU, * ), WORK( LDWORK, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SSTT22 checks a set of M eigenvalues and eigenvectors,
20 *
21 * A U = U S
22 *
23 * where A is symmetric tridiagonal, the columns of U are orthogonal,
24 * and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
25 * Two tests are performed:
26 *
27 * RESULT(1) = | U' A U - S | / ( |A| m ulp )
28 *
29 * RESULT(2) = | I - U'U | / ( m ulp )
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The size of the matrix. If it is zero, SSTT22 does nothing.
36 * It must be at least zero.
37 *
38 * M (input) INTEGER
39 * The number of eigenpairs to check. If it is zero, SSTT22
40 * does nothing. It must be at least zero.
41 *
42 * KBAND (input) INTEGER
43 * The bandwidth of the matrix S. It may only be zero or one.
44 * If zero, then S is diagonal, and SE is not referenced. If
45 * one, then S is symmetric tri-diagonal.
46 *
47 * AD (input) REAL array, dimension (N)
48 * The diagonal of the original (unfactored) matrix A. A is
49 * assumed to be symmetric tridiagonal.
50 *
51 * AE (input) REAL array, dimension (N)
52 * The off-diagonal of the original (unfactored) matrix A. A
53 * is assumed to be symmetric tridiagonal. AE(1) is ignored,
54 * AE(2) is the (1,2) and (2,1) element, etc.
55 *
56 * SD (input) REAL array, dimension (N)
57 * The diagonal of the (symmetric tri-) diagonal matrix S.
58 *
59 * SE (input) REAL array, dimension (N)
60 * The off-diagonal of the (symmetric tri-) diagonal matrix S.
61 * Not referenced if KBSND=0. If KBAND=1, then AE(1) is
62 * ignored, SE(2) is the (1,2) and (2,1) element, etc.
63 *
64 * U (input) REAL array, dimension (LDU, N)
65 * The orthogonal matrix in the decomposition.
66 *
67 * LDU (input) INTEGER
68 * The leading dimension of U. LDU must be at least N.
69 *
70 * WORK (workspace) REAL array, dimension (LDWORK, M+1)
71 *
72 * LDWORK (input) INTEGER
73 * The leading dimension of WORK. LDWORK must be at least
74 * max(1,M).
75 *
76 * RESULT (output) REAL array, dimension (2)
77 * The values computed by the two tests described above. The
78 * values are currently limited to 1/ulp, to avoid overflow.
79 *
80 * =====================================================================
81 *
82 * .. Parameters ..
83 REAL ZERO, ONE
84 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
85 * ..
86 * .. Local Scalars ..
87 INTEGER I, J, K
88 REAL ANORM, AUKJ, ULP, UNFL, WNORM
89 * ..
90 * .. External Functions ..
91 REAL SLAMCH, SLANGE, SLANSY
92 EXTERNAL SLAMCH, SLANGE, SLANSY
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL SGEMM
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC ABS, MAX, MIN, REAL
99 * ..
100 * .. Executable Statements ..
101 *
102 RESULT( 1 ) = ZERO
103 RESULT( 2 ) = ZERO
104 IF( N.LE.0 .OR. M.LE.0 )
105 $ RETURN
106 *
107 UNFL = SLAMCH( 'Safe minimum' )
108 ULP = SLAMCH( 'Epsilon' )
109 *
110 * Do Test 1
111 *
112 * Compute the 1-norm of A.
113 *
114 IF( N.GT.1 ) THEN
115 ANORM = ABS( AD( 1 ) ) + ABS( AE( 1 ) )
116 DO 10 J = 2, N - 1
117 ANORM = MAX( ANORM, ABS( AD( J ) )+ABS( AE( J ) )+
118 $ ABS( AE( J-1 ) ) )
119 10 CONTINUE
120 ANORM = MAX( ANORM, ABS( AD( N ) )+ABS( AE( N-1 ) ) )
121 ELSE
122 ANORM = ABS( AD( 1 ) )
123 END IF
124 ANORM = MAX( ANORM, UNFL )
125 *
126 * Norm of U'AU - S
127 *
128 DO 40 I = 1, M
129 DO 30 J = 1, M
130 WORK( I, J ) = ZERO
131 DO 20 K = 1, N
132 AUKJ = AD( K )*U( K, J )
133 IF( K.NE.N )
134 $ AUKJ = AUKJ + AE( K )*U( K+1, J )
135 IF( K.NE.1 )
136 $ AUKJ = AUKJ + AE( K-1 )*U( K-1, J )
137 WORK( I, J ) = WORK( I, J ) + U( K, I )*AUKJ
138 20 CONTINUE
139 30 CONTINUE
140 WORK( I, I ) = WORK( I, I ) - SD( I )
141 IF( KBAND.EQ.1 ) THEN
142 IF( I.NE.1 )
143 $ WORK( I, I-1 ) = WORK( I, I-1 ) - SE( I-1 )
144 IF( I.NE.N )
145 $ WORK( I, I+1 ) = WORK( I, I+1 ) - SE( I )
146 END IF
147 40 CONTINUE
148 *
149 WNORM = SLANSY( '1', 'L', M, WORK, M, WORK( 1, M+1 ) )
150 *
151 IF( ANORM.GT.WNORM ) THEN
152 RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
153 ELSE
154 IF( ANORM.LT.ONE ) THEN
155 RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
156 ELSE
157 RESULT( 1 ) = MIN( WNORM / ANORM, REAL( M ) ) / ( M*ULP )
158 END IF
159 END IF
160 *
161 * Do Test 2
162 *
163 * Compute U'U - I
164 *
165 CALL SGEMM( 'T', 'N', M, M, N, ONE, U, LDU, U, LDU, ZERO, WORK,
166 $ M )
167 *
168 DO 50 J = 1, M
169 WORK( J, J ) = WORK( J, J ) - ONE
170 50 CONTINUE
171 *
172 RESULT( 2 ) = MIN( REAL( M ), SLANGE( '1', M, M, WORK, M, WORK( 1,
173 $ M+1 ) ) ) / ( M*ULP )
174 *
175 RETURN
176 *
177 * End of SSTT22
178 *
179 END
2 $ LDWORK, RESULT )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER KBAND, LDU, LDWORK, M, N
10 * ..
11 * .. Array Arguments ..
12 REAL AD( * ), AE( * ), RESULT( 2 ), SD( * ),
13 $ SE( * ), U( LDU, * ), WORK( LDWORK, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * SSTT22 checks a set of M eigenvalues and eigenvectors,
20 *
21 * A U = U S
22 *
23 * where A is symmetric tridiagonal, the columns of U are orthogonal,
24 * and S is diagonal (if KBAND=0) or symmetric tridiagonal (if KBAND=1).
25 * Two tests are performed:
26 *
27 * RESULT(1) = | U' A U - S | / ( |A| m ulp )
28 *
29 * RESULT(2) = | I - U'U | / ( m ulp )
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The size of the matrix. If it is zero, SSTT22 does nothing.
36 * It must be at least zero.
37 *
38 * M (input) INTEGER
39 * The number of eigenpairs to check. If it is zero, SSTT22
40 * does nothing. It must be at least zero.
41 *
42 * KBAND (input) INTEGER
43 * The bandwidth of the matrix S. It may only be zero or one.
44 * If zero, then S is diagonal, and SE is not referenced. If
45 * one, then S is symmetric tri-diagonal.
46 *
47 * AD (input) REAL array, dimension (N)
48 * The diagonal of the original (unfactored) matrix A. A is
49 * assumed to be symmetric tridiagonal.
50 *
51 * AE (input) REAL array, dimension (N)
52 * The off-diagonal of the original (unfactored) matrix A. A
53 * is assumed to be symmetric tridiagonal. AE(1) is ignored,
54 * AE(2) is the (1,2) and (2,1) element, etc.
55 *
56 * SD (input) REAL array, dimension (N)
57 * The diagonal of the (symmetric tri-) diagonal matrix S.
58 *
59 * SE (input) REAL array, dimension (N)
60 * The off-diagonal of the (symmetric tri-) diagonal matrix S.
61 * Not referenced if KBSND=0. If KBAND=1, then AE(1) is
62 * ignored, SE(2) is the (1,2) and (2,1) element, etc.
63 *
64 * U (input) REAL array, dimension (LDU, N)
65 * The orthogonal matrix in the decomposition.
66 *
67 * LDU (input) INTEGER
68 * The leading dimension of U. LDU must be at least N.
69 *
70 * WORK (workspace) REAL array, dimension (LDWORK, M+1)
71 *
72 * LDWORK (input) INTEGER
73 * The leading dimension of WORK. LDWORK must be at least
74 * max(1,M).
75 *
76 * RESULT (output) REAL array, dimension (2)
77 * The values computed by the two tests described above. The
78 * values are currently limited to 1/ulp, to avoid overflow.
79 *
80 * =====================================================================
81 *
82 * .. Parameters ..
83 REAL ZERO, ONE
84 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
85 * ..
86 * .. Local Scalars ..
87 INTEGER I, J, K
88 REAL ANORM, AUKJ, ULP, UNFL, WNORM
89 * ..
90 * .. External Functions ..
91 REAL SLAMCH, SLANGE, SLANSY
92 EXTERNAL SLAMCH, SLANGE, SLANSY
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL SGEMM
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC ABS, MAX, MIN, REAL
99 * ..
100 * .. Executable Statements ..
101 *
102 RESULT( 1 ) = ZERO
103 RESULT( 2 ) = ZERO
104 IF( N.LE.0 .OR. M.LE.0 )
105 $ RETURN
106 *
107 UNFL = SLAMCH( 'Safe minimum' )
108 ULP = SLAMCH( 'Epsilon' )
109 *
110 * Do Test 1
111 *
112 * Compute the 1-norm of A.
113 *
114 IF( N.GT.1 ) THEN
115 ANORM = ABS( AD( 1 ) ) + ABS( AE( 1 ) )
116 DO 10 J = 2, N - 1
117 ANORM = MAX( ANORM, ABS( AD( J ) )+ABS( AE( J ) )+
118 $ ABS( AE( J-1 ) ) )
119 10 CONTINUE
120 ANORM = MAX( ANORM, ABS( AD( N ) )+ABS( AE( N-1 ) ) )
121 ELSE
122 ANORM = ABS( AD( 1 ) )
123 END IF
124 ANORM = MAX( ANORM, UNFL )
125 *
126 * Norm of U'AU - S
127 *
128 DO 40 I = 1, M
129 DO 30 J = 1, M
130 WORK( I, J ) = ZERO
131 DO 20 K = 1, N
132 AUKJ = AD( K )*U( K, J )
133 IF( K.NE.N )
134 $ AUKJ = AUKJ + AE( K )*U( K+1, J )
135 IF( K.NE.1 )
136 $ AUKJ = AUKJ + AE( K-1 )*U( K-1, J )
137 WORK( I, J ) = WORK( I, J ) + U( K, I )*AUKJ
138 20 CONTINUE
139 30 CONTINUE
140 WORK( I, I ) = WORK( I, I ) - SD( I )
141 IF( KBAND.EQ.1 ) THEN
142 IF( I.NE.1 )
143 $ WORK( I, I-1 ) = WORK( I, I-1 ) - SE( I-1 )
144 IF( I.NE.N )
145 $ WORK( I, I+1 ) = WORK( I, I+1 ) - SE( I )
146 END IF
147 40 CONTINUE
148 *
149 WNORM = SLANSY( '1', 'L', M, WORK, M, WORK( 1, M+1 ) )
150 *
151 IF( ANORM.GT.WNORM ) THEN
152 RESULT( 1 ) = ( WNORM / ANORM ) / ( M*ULP )
153 ELSE
154 IF( ANORM.LT.ONE ) THEN
155 RESULT( 1 ) = ( MIN( WNORM, M*ANORM ) / ANORM ) / ( M*ULP )
156 ELSE
157 RESULT( 1 ) = MIN( WNORM / ANORM, REAL( M ) ) / ( M*ULP )
158 END IF
159 END IF
160 *
161 * Do Test 2
162 *
163 * Compute U'U - I
164 *
165 CALL SGEMM( 'T', 'N', M, M, N, ONE, U, LDU, U, LDU, ZERO, WORK,
166 $ M )
167 *
168 DO 50 J = 1, M
169 WORK( J, J ) = WORK( J, J ) - ONE
170 50 CONTINUE
171 *
172 RESULT( 2 ) = MIN( REAL( M ), SLANGE( '1', M, M, WORK, M, WORK( 1,
173 $ M+1 ) ) ) / ( M*ULP )
174 *
175 RETURN
176 *
177 * End of SSTT22
178 *
179 END