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