1 SUBROUTINE DSPT21( ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP,
2 $ TAU, WORK, 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 ITYPE, KBAND, LDU, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), D( * ), E( * ), RESULT( 2 ), TAU( * ),
14 $ U( LDU, * ), VP( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPT21 generally checks a decomposition of the form
21 *
22 * A = U S U'
23 *
24 * where ' means transpose, A is symmetric (stored in packed format), U
25 * is orthogonal, and S is diagonal (if KBAND=0) or symmetric
26 * tridiagonal (if KBAND=1). If ITYPE=1, then U is represented as a
27 * dense matrix, otherwise the U is expressed as a product of
28 * Householder transformations, whose vectors are stored in the array
29 * "V" and whose scaling constants are in "TAU"; we shall use the
30 * letter "V" to refer to the product of Householder transformations
31 * (which should be equal to U).
32 *
33 * Specifically, if ITYPE=1, then:
34 *
35 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
36 * RESULT(2) = | I - UU' | / ( n ulp )
37 *
38 * If ITYPE=2, then:
39 *
40 * RESULT(1) = | A - V S V' | / ( |A| n ulp )
41 *
42 * If ITYPE=3, then:
43 *
44 * RESULT(1) = | I - VU' | / ( n ulp )
45 *
46 * Packed storage means that, for example, if UPLO='U', then the columns
47 * of the upper triangle of A are stored one after another, so that
48 * A(1,j+1) immediately follows A(j,j) in the array AP. Similarly, if
49 * UPLO='L', then the columns of the lower triangle of A are stored one
50 * after another in AP, so that A(j+1,j+1) immediately follows A(n,j)
51 * in the array AP. This means that A(i,j) is stored in:
52 *
53 * AP( i + j*(j-1)/2 ) if UPLO='U'
54 *
55 * AP( i + (2*n-j)*(j-1)/2 ) if UPLO='L'
56 *
57 * The array VP bears the same relation to the matrix V that A does to
58 * AP.
59 *
60 * For ITYPE > 1, the transformation U is expressed as a product
61 * of Householder transformations:
62 *
63 * If UPLO='U', then V = H(n-1)...H(1), where
64 *
65 * H(j) = I - tau(j) v(j) v(j)'
66 *
67 * and the first j-1 elements of v(j) are stored in V(1:j-1,j+1),
68 * (i.e., VP( j*(j+1)/2 + 1 : j*(j+1)/2 + j-1 ) ),
69 * the j-th element is 1, and the last n-j elements are 0.
70 *
71 * If UPLO='L', then V = H(1)...H(n-1), where
72 *
73 * H(j) = I - tau(j) v(j) v(j)'
74 *
75 * and the first j elements of v(j) are 0, the (j+1)-st is 1, and the
76 * (j+2)-nd through n-th elements are stored in V(j+2:n,j) (i.e.,
77 * in VP( (2*n-j)*(j-1)/2 + j+2 : (2*n-j)*(j-1)/2 + n ) .)
78 *
79 * Arguments
80 * =========
81 *
82 * ITYPE (input) INTEGER
83 * Specifies the type of tests to be performed.
84 * 1: U expressed as a dense orthogonal matrix:
85 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
86 * RESULT(2) = | I - UU' | / ( n ulp )
87 *
88 * 2: U expressed as a product V of Housholder transformations:
89 * RESULT(1) = | A - V S V' | / ( |A| n ulp )
90 *
91 * 3: U expressed both as a dense orthogonal matrix and
92 * as a product of Housholder transformations:
93 * RESULT(1) = | I - VU' | / ( n ulp )
94 *
95 * UPLO (input) CHARACTER
96 * If UPLO='U', AP and VP are considered to contain the upper
97 * triangle of A and V.
98 * If UPLO='L', AP and VP are considered to contain the lower
99 * triangle of A and V.
100 *
101 * N (input) INTEGER
102 * The size of the matrix. If it is zero, DSPT21 does nothing.
103 * It must be at least zero.
104 *
105 * KBAND (input) INTEGER
106 * The bandwidth of the matrix. It may only be zero or one.
107 * If zero, then S is diagonal, and E is not referenced. If
108 * one, then S is symmetric tri-diagonal.
109 *
110 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
111 * The original (unfactored) matrix. It is assumed to be
112 * symmetric, and contains the columns of just the upper
113 * triangle (UPLO='U') or only the lower triangle (UPLO='L'),
114 * packed one after another.
115 *
116 * D (input) DOUBLE PRECISION array, dimension (N)
117 * The diagonal of the (symmetric tri-) diagonal matrix.
118 *
119 * E (input) DOUBLE PRECISION array, dimension (N-1)
120 * The off-diagonal of the (symmetric tri-) diagonal matrix.
121 * E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
122 * (3,2) element, etc.
123 * Not referenced if KBAND=0.
124 *
125 * U (input) DOUBLE PRECISION array, dimension (LDU, N)
126 * If ITYPE=1 or 3, this contains the orthogonal matrix in
127 * the decomposition, expressed as a dense matrix. If ITYPE=2,
128 * then it is not referenced.
129 *
130 * LDU (input) INTEGER
131 * The leading dimension of U. LDU must be at least N and
132 * at least 1.
133 *
134 * VP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
135 * If ITYPE=2 or 3, the columns of this array contain the
136 * Householder vectors used to describe the orthogonal matrix
137 * in the decomposition, as described in purpose.
138 * *NOTE* If ITYPE=2 or 3, V is modified and restored. The
139 * subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
140 * is set to one, and later reset to its original value, during
141 * the course of the calculation.
142 * If ITYPE=1, then it is neither referenced nor modified.
143 *
144 * TAU (input) DOUBLE PRECISION array, dimension (N)
145 * If ITYPE >= 2, then TAU(j) is the scalar factor of
146 * v(j) v(j)' in the Householder transformation H(j) of
147 * the product U = H(1)...H(n-2)
148 * If ITYPE < 2, then TAU is not referenced.
149 *
150 * WORK (workspace) DOUBLE PRECISION array, dimension (N**2+N)
151 * Workspace.
152 *
153 * RESULT (output) DOUBLE PRECISION array, dimension (2)
154 * The values computed by the two tests described above. The
155 * values are currently limited to 1/ulp, to avoid overflow.
156 * RESULT(1) is always modified. RESULT(2) is modified only
157 * if ITYPE=1.
158 *
159 * =====================================================================
160 *
161 * .. Parameters ..
162 DOUBLE PRECISION ZERO, ONE, TEN
163 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
164 DOUBLE PRECISION HALF
165 PARAMETER ( HALF = 1.0D+0 / 2.0D+0 )
166 * ..
167 * .. Local Scalars ..
168 LOGICAL LOWER
169 CHARACTER CUPLO
170 INTEGER IINFO, J, JP, JP1, JR, LAP
171 DOUBLE PRECISION ANORM, TEMP, ULP, UNFL, VSAVE, WNORM
172 * ..
173 * .. External Functions ..
174 LOGICAL LSAME
175 DOUBLE PRECISION DDOT, DLAMCH, DLANGE, DLANSP
176 EXTERNAL LSAME, DDOT, DLAMCH, DLANGE, DLANSP
177 * ..
178 * .. External Subroutines ..
179 EXTERNAL DAXPY, DCOPY, DGEMM, DLACPY, DLASET, DOPMTR,
180 $ DSPMV, DSPR, DSPR2
181 * ..
182 * .. Intrinsic Functions ..
183 INTRINSIC DBLE, MAX, MIN
184 * ..
185 * .. Executable Statements ..
186 *
187 * 1) Constants
188 *
189 RESULT( 1 ) = ZERO
190 IF( ITYPE.EQ.1 )
191 $ RESULT( 2 ) = ZERO
192 IF( N.LE.0 )
193 $ RETURN
194 *
195 LAP = ( N*( N+1 ) ) / 2
196 *
197 IF( LSAME( UPLO, 'U' ) ) THEN
198 LOWER = .FALSE.
199 CUPLO = 'U'
200 ELSE
201 LOWER = .TRUE.
202 CUPLO = 'L'
203 END IF
204 *
205 UNFL = DLAMCH( 'Safe minimum' )
206 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
207 *
208 * Some Error Checks
209 *
210 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
211 RESULT( 1 ) = TEN / ULP
212 RETURN
213 END IF
214 *
215 * Do Test 1
216 *
217 * Norm of A:
218 *
219 IF( ITYPE.EQ.3 ) THEN
220 ANORM = ONE
221 ELSE
222 ANORM = MAX( DLANSP( '1', CUPLO, N, AP, WORK ), UNFL )
223 END IF
224 *
225 * Compute error matrix:
226 *
227 IF( ITYPE.EQ.1 ) THEN
228 *
229 * ITYPE=1: error = A - U S U'
230 *
231 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
232 CALL DCOPY( LAP, AP, 1, WORK, 1 )
233 *
234 DO 10 J = 1, N
235 CALL DSPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK )
236 10 CONTINUE
237 *
238 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
239 DO 20 J = 1, N - 1
240 CALL DSPR2( CUPLO, N, -E( J ), U( 1, J ), 1, U( 1, J+1 ),
241 $ 1, WORK )
242 20 CONTINUE
243 END IF
244 WNORM = DLANSP( '1', CUPLO, N, WORK, WORK( N**2+1 ) )
245 *
246 ELSE IF( ITYPE.EQ.2 ) THEN
247 *
248 * ITYPE=2: error = V S V' - A
249 *
250 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
251 *
252 IF( LOWER ) THEN
253 WORK( LAP ) = D( N )
254 DO 40 J = N - 1, 1, -1
255 JP = ( ( 2*N-J )*( J-1 ) ) / 2
256 JP1 = JP + N - J
257 IF( KBAND.EQ.1 ) THEN
258 WORK( JP+J+1 ) = ( ONE-TAU( J ) )*E( J )
259 DO 30 JR = J + 2, N
260 WORK( JP+JR ) = -TAU( J )*E( J )*VP( JP+JR )
261 30 CONTINUE
262 END IF
263 *
264 IF( TAU( J ).NE.ZERO ) THEN
265 VSAVE = VP( JP+J+1 )
266 VP( JP+J+1 ) = ONE
267 CALL DSPMV( 'L', N-J, ONE, WORK( JP1+J+1 ),
268 $ VP( JP+J+1 ), 1, ZERO, WORK( LAP+1 ), 1 )
269 TEMP = -HALF*TAU( J )*DDOT( N-J, WORK( LAP+1 ), 1,
270 $ VP( JP+J+1 ), 1 )
271 CALL DAXPY( N-J, TEMP, VP( JP+J+1 ), 1, WORK( LAP+1 ),
272 $ 1 )
273 CALL DSPR2( 'L', N-J, -TAU( J ), VP( JP+J+1 ), 1,
274 $ WORK( LAP+1 ), 1, WORK( JP1+J+1 ) )
275 VP( JP+J+1 ) = VSAVE
276 END IF
277 WORK( JP+J ) = D( J )
278 40 CONTINUE
279 ELSE
280 WORK( 1 ) = D( 1 )
281 DO 60 J = 1, N - 1
282 JP = ( J*( J-1 ) ) / 2
283 JP1 = JP + J
284 IF( KBAND.EQ.1 ) THEN
285 WORK( JP1+J ) = ( ONE-TAU( J ) )*E( J )
286 DO 50 JR = 1, J - 1
287 WORK( JP1+JR ) = -TAU( J )*E( J )*VP( JP1+JR )
288 50 CONTINUE
289 END IF
290 *
291 IF( TAU( J ).NE.ZERO ) THEN
292 VSAVE = VP( JP1+J )
293 VP( JP1+J ) = ONE
294 CALL DSPMV( 'U', J, ONE, WORK, VP( JP1+1 ), 1, ZERO,
295 $ WORK( LAP+1 ), 1 )
296 TEMP = -HALF*TAU( J )*DDOT( J, WORK( LAP+1 ), 1,
297 $ VP( JP1+1 ), 1 )
298 CALL DAXPY( J, TEMP, VP( JP1+1 ), 1, WORK( LAP+1 ),
299 $ 1 )
300 CALL DSPR2( 'U', J, -TAU( J ), VP( JP1+1 ), 1,
301 $ WORK( LAP+1 ), 1, WORK )
302 VP( JP1+J ) = VSAVE
303 END IF
304 WORK( JP1+J+1 ) = D( J+1 )
305 60 CONTINUE
306 END IF
307 *
308 DO 70 J = 1, LAP
309 WORK( J ) = WORK( J ) - AP( J )
310 70 CONTINUE
311 WNORM = DLANSP( '1', CUPLO, N, WORK, WORK( LAP+1 ) )
312 *
313 ELSE IF( ITYPE.EQ.3 ) THEN
314 *
315 * ITYPE=3: error = U V' - I
316 *
317 IF( N.LT.2 )
318 $ RETURN
319 CALL DLACPY( ' ', N, N, U, LDU, WORK, N )
320 CALL DOPMTR( 'R', CUPLO, 'T', N, N, VP, TAU, WORK, N,
321 $ WORK( N**2+1 ), IINFO )
322 IF( IINFO.NE.0 ) THEN
323 RESULT( 1 ) = TEN / ULP
324 RETURN
325 END IF
326 *
327 DO 80 J = 1, N
328 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE
329 80 CONTINUE
330 *
331 WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) )
332 END IF
333 *
334 IF( ANORM.GT.WNORM ) THEN
335 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
336 ELSE
337 IF( ANORM.LT.ONE ) THEN
338 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
339 ELSE
340 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
341 END IF
342 END IF
343 *
344 * Do Test 2
345 *
346 * Compute UU' - I
347 *
348 IF( ITYPE.EQ.1 ) THEN
349 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
350 $ N )
351 *
352 DO 90 J = 1, N
353 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE
354 90 CONTINUE
355 *
356 RESULT( 2 ) = MIN( DLANGE( '1', N, N, WORK, N,
357 $ WORK( N**2+1 ) ), DBLE( N ) ) / ( N*ULP )
358 END IF
359 *
360 RETURN
361 *
362 * End of DSPT21
363 *
364 END
2 $ TAU, WORK, 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 ITYPE, KBAND, LDU, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), D( * ), E( * ), RESULT( 2 ), TAU( * ),
14 $ U( LDU, * ), VP( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DSPT21 generally checks a decomposition of the form
21 *
22 * A = U S U'
23 *
24 * where ' means transpose, A is symmetric (stored in packed format), U
25 * is orthogonal, and S is diagonal (if KBAND=0) or symmetric
26 * tridiagonal (if KBAND=1). If ITYPE=1, then U is represented as a
27 * dense matrix, otherwise the U is expressed as a product of
28 * Householder transformations, whose vectors are stored in the array
29 * "V" and whose scaling constants are in "TAU"; we shall use the
30 * letter "V" to refer to the product of Householder transformations
31 * (which should be equal to U).
32 *
33 * Specifically, if ITYPE=1, then:
34 *
35 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
36 * RESULT(2) = | I - UU' | / ( n ulp )
37 *
38 * If ITYPE=2, then:
39 *
40 * RESULT(1) = | A - V S V' | / ( |A| n ulp )
41 *
42 * If ITYPE=3, then:
43 *
44 * RESULT(1) = | I - VU' | / ( n ulp )
45 *
46 * Packed storage means that, for example, if UPLO='U', then the columns
47 * of the upper triangle of A are stored one after another, so that
48 * A(1,j+1) immediately follows A(j,j) in the array AP. Similarly, if
49 * UPLO='L', then the columns of the lower triangle of A are stored one
50 * after another in AP, so that A(j+1,j+1) immediately follows A(n,j)
51 * in the array AP. This means that A(i,j) is stored in:
52 *
53 * AP( i + j*(j-1)/2 ) if UPLO='U'
54 *
55 * AP( i + (2*n-j)*(j-1)/2 ) if UPLO='L'
56 *
57 * The array VP bears the same relation to the matrix V that A does to
58 * AP.
59 *
60 * For ITYPE > 1, the transformation U is expressed as a product
61 * of Householder transformations:
62 *
63 * If UPLO='U', then V = H(n-1)...H(1), where
64 *
65 * H(j) = I - tau(j) v(j) v(j)'
66 *
67 * and the first j-1 elements of v(j) are stored in V(1:j-1,j+1),
68 * (i.e., VP( j*(j+1)/2 + 1 : j*(j+1)/2 + j-1 ) ),
69 * the j-th element is 1, and the last n-j elements are 0.
70 *
71 * If UPLO='L', then V = H(1)...H(n-1), where
72 *
73 * H(j) = I - tau(j) v(j) v(j)'
74 *
75 * and the first j elements of v(j) are 0, the (j+1)-st is 1, and the
76 * (j+2)-nd through n-th elements are stored in V(j+2:n,j) (i.e.,
77 * in VP( (2*n-j)*(j-1)/2 + j+2 : (2*n-j)*(j-1)/2 + n ) .)
78 *
79 * Arguments
80 * =========
81 *
82 * ITYPE (input) INTEGER
83 * Specifies the type of tests to be performed.
84 * 1: U expressed as a dense orthogonal matrix:
85 * RESULT(1) = | A - U S U' | / ( |A| n ulp ) *and*
86 * RESULT(2) = | I - UU' | / ( n ulp )
87 *
88 * 2: U expressed as a product V of Housholder transformations:
89 * RESULT(1) = | A - V S V' | / ( |A| n ulp )
90 *
91 * 3: U expressed both as a dense orthogonal matrix and
92 * as a product of Housholder transformations:
93 * RESULT(1) = | I - VU' | / ( n ulp )
94 *
95 * UPLO (input) CHARACTER
96 * If UPLO='U', AP and VP are considered to contain the upper
97 * triangle of A and V.
98 * If UPLO='L', AP and VP are considered to contain the lower
99 * triangle of A and V.
100 *
101 * N (input) INTEGER
102 * The size of the matrix. If it is zero, DSPT21 does nothing.
103 * It must be at least zero.
104 *
105 * KBAND (input) INTEGER
106 * The bandwidth of the matrix. It may only be zero or one.
107 * If zero, then S is diagonal, and E is not referenced. If
108 * one, then S is symmetric tri-diagonal.
109 *
110 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
111 * The original (unfactored) matrix. It is assumed to be
112 * symmetric, and contains the columns of just the upper
113 * triangle (UPLO='U') or only the lower triangle (UPLO='L'),
114 * packed one after another.
115 *
116 * D (input) DOUBLE PRECISION array, dimension (N)
117 * The diagonal of the (symmetric tri-) diagonal matrix.
118 *
119 * E (input) DOUBLE PRECISION array, dimension (N-1)
120 * The off-diagonal of the (symmetric tri-) diagonal matrix.
121 * E(1) is the (1,2) and (2,1) element, E(2) is the (2,3) and
122 * (3,2) element, etc.
123 * Not referenced if KBAND=0.
124 *
125 * U (input) DOUBLE PRECISION array, dimension (LDU, N)
126 * If ITYPE=1 or 3, this contains the orthogonal matrix in
127 * the decomposition, expressed as a dense matrix. If ITYPE=2,
128 * then it is not referenced.
129 *
130 * LDU (input) INTEGER
131 * The leading dimension of U. LDU must be at least N and
132 * at least 1.
133 *
134 * VP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
135 * If ITYPE=2 or 3, the columns of this array contain the
136 * Householder vectors used to describe the orthogonal matrix
137 * in the decomposition, as described in purpose.
138 * *NOTE* If ITYPE=2 or 3, V is modified and restored. The
139 * subdiagonal (if UPLO='L') or the superdiagonal (if UPLO='U')
140 * is set to one, and later reset to its original value, during
141 * the course of the calculation.
142 * If ITYPE=1, then it is neither referenced nor modified.
143 *
144 * TAU (input) DOUBLE PRECISION array, dimension (N)
145 * If ITYPE >= 2, then TAU(j) is the scalar factor of
146 * v(j) v(j)' in the Householder transformation H(j) of
147 * the product U = H(1)...H(n-2)
148 * If ITYPE < 2, then TAU is not referenced.
149 *
150 * WORK (workspace) DOUBLE PRECISION array, dimension (N**2+N)
151 * Workspace.
152 *
153 * RESULT (output) DOUBLE PRECISION array, dimension (2)
154 * The values computed by the two tests described above. The
155 * values are currently limited to 1/ulp, to avoid overflow.
156 * RESULT(1) is always modified. RESULT(2) is modified only
157 * if ITYPE=1.
158 *
159 * =====================================================================
160 *
161 * .. Parameters ..
162 DOUBLE PRECISION ZERO, ONE, TEN
163 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
164 DOUBLE PRECISION HALF
165 PARAMETER ( HALF = 1.0D+0 / 2.0D+0 )
166 * ..
167 * .. Local Scalars ..
168 LOGICAL LOWER
169 CHARACTER CUPLO
170 INTEGER IINFO, J, JP, JP1, JR, LAP
171 DOUBLE PRECISION ANORM, TEMP, ULP, UNFL, VSAVE, WNORM
172 * ..
173 * .. External Functions ..
174 LOGICAL LSAME
175 DOUBLE PRECISION DDOT, DLAMCH, DLANGE, DLANSP
176 EXTERNAL LSAME, DDOT, DLAMCH, DLANGE, DLANSP
177 * ..
178 * .. External Subroutines ..
179 EXTERNAL DAXPY, DCOPY, DGEMM, DLACPY, DLASET, DOPMTR,
180 $ DSPMV, DSPR, DSPR2
181 * ..
182 * .. Intrinsic Functions ..
183 INTRINSIC DBLE, MAX, MIN
184 * ..
185 * .. Executable Statements ..
186 *
187 * 1) Constants
188 *
189 RESULT( 1 ) = ZERO
190 IF( ITYPE.EQ.1 )
191 $ RESULT( 2 ) = ZERO
192 IF( N.LE.0 )
193 $ RETURN
194 *
195 LAP = ( N*( N+1 ) ) / 2
196 *
197 IF( LSAME( UPLO, 'U' ) ) THEN
198 LOWER = .FALSE.
199 CUPLO = 'U'
200 ELSE
201 LOWER = .TRUE.
202 CUPLO = 'L'
203 END IF
204 *
205 UNFL = DLAMCH( 'Safe minimum' )
206 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
207 *
208 * Some Error Checks
209 *
210 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
211 RESULT( 1 ) = TEN / ULP
212 RETURN
213 END IF
214 *
215 * Do Test 1
216 *
217 * Norm of A:
218 *
219 IF( ITYPE.EQ.3 ) THEN
220 ANORM = ONE
221 ELSE
222 ANORM = MAX( DLANSP( '1', CUPLO, N, AP, WORK ), UNFL )
223 END IF
224 *
225 * Compute error matrix:
226 *
227 IF( ITYPE.EQ.1 ) THEN
228 *
229 * ITYPE=1: error = A - U S U'
230 *
231 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
232 CALL DCOPY( LAP, AP, 1, WORK, 1 )
233 *
234 DO 10 J = 1, N
235 CALL DSPR( CUPLO, N, -D( J ), U( 1, J ), 1, WORK )
236 10 CONTINUE
237 *
238 IF( N.GT.1 .AND. KBAND.EQ.1 ) THEN
239 DO 20 J = 1, N - 1
240 CALL DSPR2( CUPLO, N, -E( J ), U( 1, J ), 1, U( 1, J+1 ),
241 $ 1, WORK )
242 20 CONTINUE
243 END IF
244 WNORM = DLANSP( '1', CUPLO, N, WORK, WORK( N**2+1 ) )
245 *
246 ELSE IF( ITYPE.EQ.2 ) THEN
247 *
248 * ITYPE=2: error = V S V' - A
249 *
250 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N )
251 *
252 IF( LOWER ) THEN
253 WORK( LAP ) = D( N )
254 DO 40 J = N - 1, 1, -1
255 JP = ( ( 2*N-J )*( J-1 ) ) / 2
256 JP1 = JP + N - J
257 IF( KBAND.EQ.1 ) THEN
258 WORK( JP+J+1 ) = ( ONE-TAU( J ) )*E( J )
259 DO 30 JR = J + 2, N
260 WORK( JP+JR ) = -TAU( J )*E( J )*VP( JP+JR )
261 30 CONTINUE
262 END IF
263 *
264 IF( TAU( J ).NE.ZERO ) THEN
265 VSAVE = VP( JP+J+1 )
266 VP( JP+J+1 ) = ONE
267 CALL DSPMV( 'L', N-J, ONE, WORK( JP1+J+1 ),
268 $ VP( JP+J+1 ), 1, ZERO, WORK( LAP+1 ), 1 )
269 TEMP = -HALF*TAU( J )*DDOT( N-J, WORK( LAP+1 ), 1,
270 $ VP( JP+J+1 ), 1 )
271 CALL DAXPY( N-J, TEMP, VP( JP+J+1 ), 1, WORK( LAP+1 ),
272 $ 1 )
273 CALL DSPR2( 'L', N-J, -TAU( J ), VP( JP+J+1 ), 1,
274 $ WORK( LAP+1 ), 1, WORK( JP1+J+1 ) )
275 VP( JP+J+1 ) = VSAVE
276 END IF
277 WORK( JP+J ) = D( J )
278 40 CONTINUE
279 ELSE
280 WORK( 1 ) = D( 1 )
281 DO 60 J = 1, N - 1
282 JP = ( J*( J-1 ) ) / 2
283 JP1 = JP + J
284 IF( KBAND.EQ.1 ) THEN
285 WORK( JP1+J ) = ( ONE-TAU( J ) )*E( J )
286 DO 50 JR = 1, J - 1
287 WORK( JP1+JR ) = -TAU( J )*E( J )*VP( JP1+JR )
288 50 CONTINUE
289 END IF
290 *
291 IF( TAU( J ).NE.ZERO ) THEN
292 VSAVE = VP( JP1+J )
293 VP( JP1+J ) = ONE
294 CALL DSPMV( 'U', J, ONE, WORK, VP( JP1+1 ), 1, ZERO,
295 $ WORK( LAP+1 ), 1 )
296 TEMP = -HALF*TAU( J )*DDOT( J, WORK( LAP+1 ), 1,
297 $ VP( JP1+1 ), 1 )
298 CALL DAXPY( J, TEMP, VP( JP1+1 ), 1, WORK( LAP+1 ),
299 $ 1 )
300 CALL DSPR2( 'U', J, -TAU( J ), VP( JP1+1 ), 1,
301 $ WORK( LAP+1 ), 1, WORK )
302 VP( JP1+J ) = VSAVE
303 END IF
304 WORK( JP1+J+1 ) = D( J+1 )
305 60 CONTINUE
306 END IF
307 *
308 DO 70 J = 1, LAP
309 WORK( J ) = WORK( J ) - AP( J )
310 70 CONTINUE
311 WNORM = DLANSP( '1', CUPLO, N, WORK, WORK( LAP+1 ) )
312 *
313 ELSE IF( ITYPE.EQ.3 ) THEN
314 *
315 * ITYPE=3: error = U V' - I
316 *
317 IF( N.LT.2 )
318 $ RETURN
319 CALL DLACPY( ' ', N, N, U, LDU, WORK, N )
320 CALL DOPMTR( 'R', CUPLO, 'T', N, N, VP, TAU, WORK, N,
321 $ WORK( N**2+1 ), IINFO )
322 IF( IINFO.NE.0 ) THEN
323 RESULT( 1 ) = TEN / ULP
324 RETURN
325 END IF
326 *
327 DO 80 J = 1, N
328 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE
329 80 CONTINUE
330 *
331 WNORM = DLANGE( '1', N, N, WORK, N, WORK( N**2+1 ) )
332 END IF
333 *
334 IF( ANORM.GT.WNORM ) THEN
335 RESULT( 1 ) = ( WNORM / ANORM ) / ( N*ULP )
336 ELSE
337 IF( ANORM.LT.ONE ) THEN
338 RESULT( 1 ) = ( MIN( WNORM, N*ANORM ) / ANORM ) / ( N*ULP )
339 ELSE
340 RESULT( 1 ) = MIN( WNORM / ANORM, DBLE( N ) ) / ( N*ULP )
341 END IF
342 END IF
343 *
344 * Do Test 2
345 *
346 * Compute UU' - I
347 *
348 IF( ITYPE.EQ.1 ) THEN
349 CALL DGEMM( 'N', 'C', N, N, N, ONE, U, LDU, U, LDU, ZERO, WORK,
350 $ N )
351 *
352 DO 90 J = 1, N
353 WORK( ( N+1 )*( J-1 )+1 ) = WORK( ( N+1 )*( J-1 )+1 ) - ONE
354 90 CONTINUE
355 *
356 RESULT( 2 ) = MIN( DLANGE( '1', N, N, WORK, N,
357 $ WORK( N**2+1 ) ), DBLE( N ) ) / ( N*ULP )
358 END IF
359 *
360 RETURN
361 *
362 * End of DSPT21
363 *
364 END