1 SUBROUTINE ZHBT21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK,
2 $ 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 CHARACTER UPLO
10 INTEGER KA, KS, LDA, LDU, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * )
14 COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHBT21 generally checks a decomposition of the form
21 *
22 * A = U S U*
23 *
24 * where * means conjugate transpose, A is hermitian banded, U is
25 * unitary, and S is diagonal (if KS=0) or symmetric
26 * tridiagonal (if KS=1).
27 *
28 * Specifically:
29 *
30 * RESULT(1) = | A - U S U* | / ( |A| n ulp ) *and*
31 * RESULT(2) = | I - UU* | / ( n ulp )
32 *
33 * Arguments
34 * =========
35 *
36 * UPLO (input) CHARACTER
37 * If UPLO='U', the upper triangle of A and V will be used and
38 * the (strictly) lower triangle will not be referenced.
39 * If UPLO='L', the lower triangle of A and V will be used and
40 * the (strictly) upper triangle will not be referenced.
41 *
42 * N (input) INTEGER
43 * The size of the matrix. If it is zero, ZHBT21 does nothing.
44 * It must be at least zero.
45 *
46 * KA (input) INTEGER
47 * The bandwidth of the matrix A. It must be at least zero. If
48 * it is larger than N-1, then max( 0, N-1 ) will be used.
49 *
50 * KS (input) INTEGER
51 * The bandwidth of the matrix S. It may only be zero or one.
52 * If zero, then S is diagonal, and E is not referenced. If
53 * one, then S is symmetric tri-diagonal.
54 *
55 * A (input) COMPLEX*16 array, dimension (LDA, N)
56 * The original (unfactored) matrix. It is assumed to be
57 * hermitian, and only the upper (UPLO='U') or only the lower
58 * (UPLO='L') will be referenced.
59 *
60 * LDA (input) INTEGER
61 * The leading dimension of A. It must be at least 1
62 * and at least min( KA, N-1 ).
63 *
64 * D (input) DOUBLE PRECISION array, dimension (N)
65 * The diagonal of the (symmetric tri-) diagonal matrix S.
66 *
67 * E (input) DOUBLE PRECISION array, dimension (N-1)
68 * The off-diagonal of the (symmetric tri-) diagonal matrix S.
69 * E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
70 * (3,2) element, etc.
71 * Not referenced if KS=0.
72 *
73 * U (input) COMPLEX*16 array, dimension (LDU, N)
74 * The unitary matrix in the decomposition, expressed as a
75 * dense matrix (i.e., not as a product of Householder
76 * transformations, Givens transformations, etc.)
77 *
78 * LDU (input) INTEGER
79 * The leading dimension of U. LDU must be at least N and
80 * at least 1.
81 *
82 * WORK (workspace) COMPLEX*16 array, dimension (N**2)
83 *
84 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
85 *
86 * RESULT (output) DOUBLE PRECISION array, dimension (2)
87 * The values computed by the two tests described above. The
88 * values are currently limited to 1/ulp, to avoid overflow.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 COMPLEX*16 CZERO, CONE
94 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
95 $ CONE = ( 1.0D+0, 0.0D+0 ) )
96 DOUBLE PRECISION ZERO, ONE
97 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
98 * ..
99 * .. Local Scalars ..
100 LOGICAL LOWER
101 CHARACTER CUPLO
102 INTEGER IKA, J, JC, JR
103 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
104 * ..
105 * .. External Functions ..
106 LOGICAL LSAME
107 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHB, ZLANHP
108 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHB, ZLANHP
109 * ..
110 * .. External Subroutines ..
111 EXTERNAL ZGEMM, ZHPR, ZHPR2
112 * ..
113 * .. Intrinsic Functions ..
114 INTRINSIC DBLE, DCMPLX, MAX, MIN
115 * ..
116 * .. Executable Statements ..
117 *
118 * Constants
119 *
120 RESULT( 1 ) = ZERO
121 RESULT( 2 ) = ZERO
122 IF( N.LE.0 )
123 $ RETURN
124 *
125 IKA = MAX( 0, MIN( N-1, KA ) )
126 *
127 IF( LSAME( UPLO, 'U' ) ) THEN
128 LOWER = .FALSE.
129 CUPLO = 'U'
130 ELSE
131 LOWER = .TRUE.
132 CUPLO = 'L'
133 END IF
134 *
135 UNFL = DLAMCH( 'Safe minimum' )
136 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
137 *
138 * Some Error Checks
139 *
140 * Do Test 1
141 *
142 * Norm of A:
143 *
144 ANORM = MAX( ZLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL )
145 *
146 * Compute error matrix: Error = A - U S U*
147 *
148 * Copy A from SB to SP storage format.
149 *
150 J = 0
151 DO 50 JC = 1, N
152 IF( LOWER ) THEN
153 DO 10 JR = 1, MIN( IKA+1, N+1-JC )
154 J = J + 1
155 WORK( J ) = A( JR, JC )
156 10 CONTINUE
157 DO 20 JR = IKA + 2, N + 1 - JC
158 J = J + 1
159 WORK( J ) = ZERO
160 20 CONTINUE
161 ELSE
162 DO 30 JR = IKA + 2, JC
163 J = J + 1
164 WORK( J ) = ZERO
165 30 CONTINUE
166 DO 40 JR = MIN( IKA, JC-1 ), 0, -1
167 J = J + 1
168 WORK( J ) = A( IKA+1-JR, JC )
169 40 CONTINUE
170 END IF
171 50 CONTINUE
172 *
173 DO 60 J = 1, N
174 CALL ZHPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK )
175 60 CONTINUE
176 *
177 IF( N.GT.1 .AND. KS.EQ.1 ) THEN
178 DO 70 J = 1, N - 1
179 CALL ZHPR2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1,
180 $ U( 1, J+1 ), 1, WORK )
181 70 CONTINUE
182 END IF
183 WNORM = ZLANHP( '1', CUPLO, N, WORK, RWORK )
184 *
185 IF( ANORM.GT.WNORM ) THEN
186 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
187 ELSE
188 IF( ANORM.LT.ONE ) THEN
189 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
190 ELSE
191 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
192 END IF
193 END IF
194 *
195 * Do Test 2
196 *
197 * Compute UU* - I
198 *
199 CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
200 $ N )
201 *
202 DO 80 J = 1, N
203 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
204 80 CONTINUE
205 *
206 RESULT( 2 ) = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ),
207 $ DBLE( N ) ) / ( N*ULP )
208 *
209 RETURN
210 *
211 * End of ZHBT21
212 *
213 END
2 $ 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 CHARACTER UPLO
10 INTEGER KA, KS, LDA, LDU, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), RESULT( 2 ), RWORK( * )
14 COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZHBT21 generally checks a decomposition of the form
21 *
22 * A = U S U*
23 *
24 * where * means conjugate transpose, A is hermitian banded, U is
25 * unitary, and S is diagonal (if KS=0) or symmetric
26 * tridiagonal (if KS=1).
27 *
28 * Specifically:
29 *
30 * RESULT(1) = | A - U S U* | / ( |A| n ulp ) *and*
31 * RESULT(2) = | I - UU* | / ( n ulp )
32 *
33 * Arguments
34 * =========
35 *
36 * UPLO (input) CHARACTER
37 * If UPLO='U', the upper triangle of A and V will be used and
38 * the (strictly) lower triangle will not be referenced.
39 * If UPLO='L', the lower triangle of A and V will be used and
40 * the (strictly) upper triangle will not be referenced.
41 *
42 * N (input) INTEGER
43 * The size of the matrix. If it is zero, ZHBT21 does nothing.
44 * It must be at least zero.
45 *
46 * KA (input) INTEGER
47 * The bandwidth of the matrix A. It must be at least zero. If
48 * it is larger than N-1, then max( 0, N-1 ) will be used.
49 *
50 * KS (input) INTEGER
51 * The bandwidth of the matrix S. It may only be zero or one.
52 * If zero, then S is diagonal, and E is not referenced. If
53 * one, then S is symmetric tri-diagonal.
54 *
55 * A (input) COMPLEX*16 array, dimension (LDA, N)
56 * The original (unfactored) matrix. It is assumed to be
57 * hermitian, and only the upper (UPLO='U') or only the lower
58 * (UPLO='L') will be referenced.
59 *
60 * LDA (input) INTEGER
61 * The leading dimension of A. It must be at least 1
62 * and at least min( KA, N-1 ).
63 *
64 * D (input) DOUBLE PRECISION array, dimension (N)
65 * The diagonal of the (symmetric tri-) diagonal matrix S.
66 *
67 * E (input) DOUBLE PRECISION array, dimension (N-1)
68 * The off-diagonal of the (symmetric tri-) diagonal matrix S.
69 * E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
70 * (3,2) element, etc.
71 * Not referenced if KS=0.
72 *
73 * U (input) COMPLEX*16 array, dimension (LDU, N)
74 * The unitary matrix in the decomposition, expressed as a
75 * dense matrix (i.e., not as a product of Householder
76 * transformations, Givens transformations, etc.)
77 *
78 * LDU (input) INTEGER
79 * The leading dimension of U. LDU must be at least N and
80 * at least 1.
81 *
82 * WORK (workspace) COMPLEX*16 array, dimension (N**2)
83 *
84 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
85 *
86 * RESULT (output) DOUBLE PRECISION array, dimension (2)
87 * The values computed by the two tests described above. The
88 * values are currently limited to 1/ulp, to avoid overflow.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 COMPLEX*16 CZERO, CONE
94 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
95 $ CONE = ( 1.0D+0, 0.0D+0 ) )
96 DOUBLE PRECISION ZERO, ONE
97 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
98 * ..
99 * .. Local Scalars ..
100 LOGICAL LOWER
101 CHARACTER CUPLO
102 INTEGER IKA, J, JC, JR
103 DOUBLE PRECISION ANORM, ULP, UNFL, WNORM
104 * ..
105 * .. External Functions ..
106 LOGICAL LSAME
107 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHB, ZLANHP
108 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHB, ZLANHP
109 * ..
110 * .. External Subroutines ..
111 EXTERNAL ZGEMM, ZHPR, ZHPR2
112 * ..
113 * .. Intrinsic Functions ..
114 INTRINSIC DBLE, DCMPLX, MAX, MIN
115 * ..
116 * .. Executable Statements ..
117 *
118 * Constants
119 *
120 RESULT( 1 ) = ZERO
121 RESULT( 2 ) = ZERO
122 IF( N.LE.0 )
123 $ RETURN
124 *
125 IKA = MAX( 0, MIN( N-1, KA ) )
126 *
127 IF( LSAME( UPLO, 'U' ) ) THEN
128 LOWER = .FALSE.
129 CUPLO = 'U'
130 ELSE
131 LOWER = .TRUE.
132 CUPLO = 'L'
133 END IF
134 *
135 UNFL = DLAMCH( 'Safe minimum' )
136 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
137 *
138 * Some Error Checks
139 *
140 * Do Test 1
141 *
142 * Norm of A:
143 *
144 ANORM = MAX( ZLANHB( '1', CUPLO, N, IKA, A, LDA, RWORK ), UNFL )
145 *
146 * Compute error matrix: Error = A - U S U*
147 *
148 * Copy A from SB to SP storage format.
149 *
150 J = 0
151 DO 50 JC = 1, N
152 IF( LOWER ) THEN
153 DO 10 JR = 1, MIN( IKA+1, N+1-JC )
154 J = J + 1
155 WORK( J ) = A( JR, JC )
156 10 CONTINUE
157 DO 20 JR = IKA + 2, N + 1 - JC
158 J = J + 1
159 WORK( J ) = ZERO
160 20 CONTINUE
161 ELSE
162 DO 30 JR = IKA + 2, JC
163 J = J + 1
164 WORK( J ) = ZERO
165 30 CONTINUE
166 DO 40 JR = MIN( IKA, JC-1 ), 0, -1
167 J = J + 1
168 WORK( J ) = A( IKA+1-JR, JC )
169 40 CONTINUE
170 END IF
171 50 CONTINUE
172 *
173 DO 60 J = 1, N
174 CALL ZHPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK )
175 60 CONTINUE
176 *
177 IF( N.GT.1 .AND. KS.EQ.1 ) THEN
178 DO 70 J = 1, N - 1
179 CALL ZHPR2( CUPLO, N, -DCMPLX( E( J ) ), U( 1, J ), 1,
180 $ U( 1, J+1 ), 1, WORK )
181 70 CONTINUE
182 END IF
183 WNORM = ZLANHP( '1', CUPLO, N, WORK, RWORK )
184 *
185 IF( ANORM.GT.WNORM ) THEN
186 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
187 ELSE
188 IF( ANORM.LT.ONE ) THEN
189 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
190 ELSE
191 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
192 END IF
193 END IF
194 *
195 * Do Test 2
196 *
197 * Compute UU* - I
198 *
199 CALL ZGEMM( 'N', 'C', N, N, N, CONE, U, LDU, U, LDU, CZERO, WORK,
200 $ N )
201 *
202 DO 80 J = 1, N
203 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - CONE
204 80 CONTINUE
205 *
206 RESULT( 2 ) = MIN( ZLANGE( '1', N, N, WORK, N, RWORK ),
207 $ DBLE( N ) ) / ( N*ULP )
208 *
209 RETURN
210 *
211 * End of ZHBT21
212 *
213 END