1 SUBROUTINE DGET39( RMAX, LMAX, NINFO, KNT )
2 *
3 * -- LAPACK 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 KNT, LMAX, NINFO
9 DOUBLE PRECISION RMAX
10 * ..
11 *
12 * Purpose
13 * =======
14 *
15 * DGET39 tests DLAQTR, a routine for solving the real or
16 * special complex quasi upper triangular system
17 *
18 * op(T)*p = scale*c,
19 * or
20 * op(T + iB)*(p+iq) = scale*(c+id),
21 *
22 * in real arithmetic. T is upper quasi-triangular.
23 * If it is complex, then the first diagonal block of T must be
24 * 1 by 1, B has the special structure
25 *
26 * B = [ b(1) b(2) ... b(n) ]
27 * [ w ]
28 * [ w ]
29 * [ . ]
30 * [ w ]
31 *
32 * op(A) = A or A', where A' denotes the conjugate transpose of
33 * the matrix A.
34 *
35 * On input, X = [ c ]. On output, X = [ p ].
36 * [ d ] [ q ]
37 *
38 * Scale is an output less than or equal to 1, chosen to avoid
39 * overflow in X.
40 * This subroutine is specially designed for the condition number
41 * estimation in the eigenproblem routine DTRSNA.
42 *
43 * The test code verifies that the following residual is order 1:
44 *
45 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
46 * -----------------------------------------
47 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
48 * (||T||+||B||)*smlnum/ulp,
49 * smlnum)
50 *
51 * (The (||T||+||B||)*smlnum/ulp term accounts for possible
52 * (gradual or nongradual) underflow in x1 and x2.)
53 *
54 * Arguments
55 * ==========
56 *
57 * RMAX (output) DOUBLE PRECISION
58 * Value of the largest test ratio.
59 *
60 * LMAX (output) INTEGER
61 * Example number where largest test ratio achieved.
62 *
63 * NINFO (output) INTEGER
64 * Number of examples where INFO is nonzero.
65 *
66 * KNT (output) INTEGER
67 * Total number of examples tested.
68 *
69 * =====================================================================
70 *
71 * .. Parameters ..
72 INTEGER LDT, LDT2
73 PARAMETER ( LDT = 10, LDT2 = 2*LDT )
74 DOUBLE PRECISION ZERO, ONE
75 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
76 * ..
77 * .. Local Scalars ..
78 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
79 $ NDIM
80 DOUBLE PRECISION BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
81 $ SCALE, SMLNUM, W, XNORM
82 * ..
83 * .. External Functions ..
84 INTEGER IDAMAX
85 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
86 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
87 * ..
88 * .. External Subroutines ..
89 EXTERNAL DCOPY, DGEMV, DLABAD, DLAQTR
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC ABS, COS, DBLE, MAX, SIN, SQRT
93 * ..
94 * .. Local Arrays ..
95 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 )
96 DOUBLE PRECISION B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
97 $ VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
98 $ VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
99 * ..
100 * .. Data statements ..
101 DATA IDIM / 4, 5*5 /
102 DATA IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
103 $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
104 $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
105 $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
106 $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
107 $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
108 $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
109 $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
110 * ..
111 * .. Executable Statements ..
112 *
113 * Get machine parameters
114 *
115 EPS = DLAMCH( 'P' )
116 SMLNUM = DLAMCH( 'S' )
117 BIGNUM = ONE / SMLNUM
118 CALL DLABAD( SMLNUM, BIGNUM )
119 *
120 * Set up test case parameters
121 *
122 VM1( 1 ) = ONE
123 VM1( 2 ) = SQRT( SMLNUM )
124 VM1( 3 ) = SQRT( VM1( 2 ) )
125 VM1( 4 ) = SQRT( BIGNUM )
126 VM1( 5 ) = SQRT( VM1( 4 ) )
127 *
128 VM2( 1 ) = ONE
129 VM2( 2 ) = SQRT( SMLNUM )
130 VM2( 3 ) = SQRT( VM2( 2 ) )
131 VM2( 4 ) = SQRT( BIGNUM )
132 VM2( 5 ) = SQRT( VM2( 4 ) )
133 *
134 VM3( 1 ) = ONE
135 VM3( 2 ) = SQRT( SMLNUM )
136 VM3( 3 ) = SQRT( VM3( 2 ) )
137 VM3( 4 ) = SQRT( BIGNUM )
138 VM3( 5 ) = SQRT( VM3( 4 ) )
139 *
140 VM4( 1 ) = ONE
141 VM4( 2 ) = SQRT( SMLNUM )
142 VM4( 3 ) = SQRT( VM4( 2 ) )
143 VM4( 4 ) = SQRT( BIGNUM )
144 VM4( 5 ) = SQRT( VM4( 4 ) )
145 *
146 VM5( 1 ) = ONE
147 VM5( 2 ) = EPS
148 VM5( 3 ) = SQRT( SMLNUM )
149 *
150 * Initalization
151 *
152 KNT = 0
153 RMAX = ZERO
154 NINFO = 0
155 SMLNUM = SMLNUM / EPS
156 *
157 * Begin test loop
158 *
159 DO 140 IVM5 = 1, 3
160 DO 130 IVM4 = 1, 5
161 DO 120 IVM3 = 1, 5
162 DO 110 IVM2 = 1, 5
163 DO 100 IVM1 = 1, 5
164 DO 90 NDIM = 1, 6
165 *
166 N = IDIM( NDIM )
167 DO 20 I = 1, N
168 DO 10 J = 1, N
169 T( I, J ) = DBLE( IVAL( I, J, NDIM ) )*
170 $ VM1( IVM1 )
171 IF( I.GE.J )
172 $ T( I, J ) = T( I, J )*VM5( IVM5 )
173 10 CONTINUE
174 20 CONTINUE
175 *
176 W = ONE*VM2( IVM2 )
177 *
178 DO 30 I = 1, N
179 B( I ) = COS( DBLE( I ) )*VM3( IVM3 )
180 30 CONTINUE
181 *
182 DO 40 I = 1, 2*N
183 D( I ) = SIN( DBLE( I ) )*VM4( IVM4 )
184 40 CONTINUE
185 *
186 NORM = DLANGE( '1', N, N, T, LDT, WORK )
187 K = IDAMAX( N, B, 1 )
188 NORMTB = NORM + ABS( B( K ) ) + ABS( W )
189 *
190 CALL DCOPY( N, D, 1, X, 1 )
191 KNT = KNT + 1
192 CALL DLAQTR( .FALSE., .TRUE., N, T, LDT, DUM,
193 $ DUMM, SCALE, X, WORK, INFO )
194 IF( INFO.NE.0 )
195 $ NINFO = NINFO + 1
196 *
197 * || T*x - scale*d || /
198 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
199 *
200 CALL DCOPY( N, D, 1, Y, 1 )
201 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
202 $ X, 1, -SCALE, Y, 1 )
203 XNORM = DASUM( N, X, 1 )
204 RESID = DASUM( N, Y, 1 )
205 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
206 $ ( NORM*EPS )*XNORM )
207 RESID = RESID / DOMIN
208 IF( RESID.GT.RMAX ) THEN
209 RMAX = RESID
210 LMAX = KNT
211 END IF
212 *
213 CALL DCOPY( N, D, 1, X, 1 )
214 KNT = KNT + 1
215 CALL DLAQTR( .TRUE., .TRUE., N, T, LDT, DUM,
216 $ DUMM, SCALE, X, WORK, INFO )
217 IF( INFO.NE.0 )
218 $ NINFO = NINFO + 1
219 *
220 * || T*x - scale*d || /
221 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
222 *
223 CALL DCOPY( N, D, 1, Y, 1 )
224 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
225 $ 1, -SCALE, Y, 1 )
226 XNORM = DASUM( N, X, 1 )
227 RESID = DASUM( N, Y, 1 )
228 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
229 $ ( NORM*EPS )*XNORM )
230 RESID = RESID / DOMIN
231 IF( RESID.GT.RMAX ) THEN
232 RMAX = RESID
233 LMAX = KNT
234 END IF
235 *
236 CALL DCOPY( 2*N, D, 1, X, 1 )
237 KNT = KNT + 1
238 CALL DLAQTR( .FALSE., .FALSE., N, T, LDT, B, W,
239 $ SCALE, X, WORK, INFO )
240 IF( INFO.NE.0 )
241 $ NINFO = NINFO + 1
242 *
243 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
244 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
245 * smlnum/ulp * (||T||+||B||), smlnum )
246 *
247 *
248 CALL DCOPY( 2*N, D, 1, Y, 1 )
249 Y( 1 ) = DDOT( N, B, 1, X( 1+N ), 1 ) +
250 $ SCALE*Y( 1 )
251 DO 50 I = 2, N
252 Y( I ) = W*X( I+N ) + SCALE*Y( I )
253 50 CONTINUE
254 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
255 $ X, 1, -ONE, Y, 1 )
256 *
257 Y( 1+N ) = DDOT( N, B, 1, X, 1 ) -
258 $ SCALE*Y( 1+N )
259 DO 60 I = 2, N
260 Y( I+N ) = W*X( I ) - SCALE*Y( I+N )
261 60 CONTINUE
262 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
263 $ X( 1+N ), 1, ONE, Y( 1+N ), 1 )
264 *
265 RESID = DASUM( 2*N, Y, 1 )
266 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
267 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
268 RESID = RESID / DOMIN
269 IF( RESID.GT.RMAX ) THEN
270 RMAX = RESID
271 LMAX = KNT
272 END IF
273 *
274 CALL DCOPY( 2*N, D, 1, X, 1 )
275 KNT = KNT + 1
276 CALL DLAQTR( .TRUE., .FALSE., N, T, LDT, B, W,
277 $ SCALE, X, WORK, INFO )
278 IF( INFO.NE.0 )
279 $ NINFO = NINFO + 1
280 *
281 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
282 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
283 * smlnum/ulp * (||T||+||B||), smlnum )
284 *
285 CALL DCOPY( 2*N, D, 1, Y, 1 )
286 Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 )
287 DO 70 I = 2, N
288 Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) -
289 $ SCALE*Y( I )
290 70 CONTINUE
291 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
292 $ 1, ONE, Y, 1 )
293 *
294 Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N )
295 DO 80 I = 2, N
296 Y( I+N ) = B( I )*X( 1 ) + W*X( I ) +
297 $ SCALE*Y( I+N )
298 80 CONTINUE
299 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT,
300 $ X( 1+N ), 1, -ONE, Y( 1+N ), 1 )
301 *
302 RESID = DASUM( 2*N, Y, 1 )
303 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
304 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
305 RESID = RESID / DOMIN
306 IF( RESID.GT.RMAX ) THEN
307 RMAX = RESID
308 LMAX = KNT
309 END IF
310 *
311 90 CONTINUE
312 100 CONTINUE
313 110 CONTINUE
314 120 CONTINUE
315 130 CONTINUE
316 140 CONTINUE
317 *
318 RETURN
319 *
320 * End of DGET39
321 *
322 END
2 *
3 * -- LAPACK 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 KNT, LMAX, NINFO
9 DOUBLE PRECISION RMAX
10 * ..
11 *
12 * Purpose
13 * =======
14 *
15 * DGET39 tests DLAQTR, a routine for solving the real or
16 * special complex quasi upper triangular system
17 *
18 * op(T)*p = scale*c,
19 * or
20 * op(T + iB)*(p+iq) = scale*(c+id),
21 *
22 * in real arithmetic. T is upper quasi-triangular.
23 * If it is complex, then the first diagonal block of T must be
24 * 1 by 1, B has the special structure
25 *
26 * B = [ b(1) b(2) ... b(n) ]
27 * [ w ]
28 * [ w ]
29 * [ . ]
30 * [ w ]
31 *
32 * op(A) = A or A', where A' denotes the conjugate transpose of
33 * the matrix A.
34 *
35 * On input, X = [ c ]. On output, X = [ p ].
36 * [ d ] [ q ]
37 *
38 * Scale is an output less than or equal to 1, chosen to avoid
39 * overflow in X.
40 * This subroutine is specially designed for the condition number
41 * estimation in the eigenproblem routine DTRSNA.
42 *
43 * The test code verifies that the following residual is order 1:
44 *
45 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
46 * -----------------------------------------
47 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
48 * (||T||+||B||)*smlnum/ulp,
49 * smlnum)
50 *
51 * (The (||T||+||B||)*smlnum/ulp term accounts for possible
52 * (gradual or nongradual) underflow in x1 and x2.)
53 *
54 * Arguments
55 * ==========
56 *
57 * RMAX (output) DOUBLE PRECISION
58 * Value of the largest test ratio.
59 *
60 * LMAX (output) INTEGER
61 * Example number where largest test ratio achieved.
62 *
63 * NINFO (output) INTEGER
64 * Number of examples where INFO is nonzero.
65 *
66 * KNT (output) INTEGER
67 * Total number of examples tested.
68 *
69 * =====================================================================
70 *
71 * .. Parameters ..
72 INTEGER LDT, LDT2
73 PARAMETER ( LDT = 10, LDT2 = 2*LDT )
74 DOUBLE PRECISION ZERO, ONE
75 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
76 * ..
77 * .. Local Scalars ..
78 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
79 $ NDIM
80 DOUBLE PRECISION BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
81 $ SCALE, SMLNUM, W, XNORM
82 * ..
83 * .. External Functions ..
84 INTEGER IDAMAX
85 DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
86 EXTERNAL IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
87 * ..
88 * .. External Subroutines ..
89 EXTERNAL DCOPY, DGEMV, DLABAD, DLAQTR
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC ABS, COS, DBLE, MAX, SIN, SQRT
93 * ..
94 * .. Local Arrays ..
95 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 )
96 DOUBLE PRECISION B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
97 $ VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
98 $ VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
99 * ..
100 * .. Data statements ..
101 DATA IDIM / 4, 5*5 /
102 DATA IVAL / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
103 $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
104 $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
105 $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
106 $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
107 $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
108 $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
109 $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
110 * ..
111 * .. Executable Statements ..
112 *
113 * Get machine parameters
114 *
115 EPS = DLAMCH( 'P' )
116 SMLNUM = DLAMCH( 'S' )
117 BIGNUM = ONE / SMLNUM
118 CALL DLABAD( SMLNUM, BIGNUM )
119 *
120 * Set up test case parameters
121 *
122 VM1( 1 ) = ONE
123 VM1( 2 ) = SQRT( SMLNUM )
124 VM1( 3 ) = SQRT( VM1( 2 ) )
125 VM1( 4 ) = SQRT( BIGNUM )
126 VM1( 5 ) = SQRT( VM1( 4 ) )
127 *
128 VM2( 1 ) = ONE
129 VM2( 2 ) = SQRT( SMLNUM )
130 VM2( 3 ) = SQRT( VM2( 2 ) )
131 VM2( 4 ) = SQRT( BIGNUM )
132 VM2( 5 ) = SQRT( VM2( 4 ) )
133 *
134 VM3( 1 ) = ONE
135 VM3( 2 ) = SQRT( SMLNUM )
136 VM3( 3 ) = SQRT( VM3( 2 ) )
137 VM3( 4 ) = SQRT( BIGNUM )
138 VM3( 5 ) = SQRT( VM3( 4 ) )
139 *
140 VM4( 1 ) = ONE
141 VM4( 2 ) = SQRT( SMLNUM )
142 VM4( 3 ) = SQRT( VM4( 2 ) )
143 VM4( 4 ) = SQRT( BIGNUM )
144 VM4( 5 ) = SQRT( VM4( 4 ) )
145 *
146 VM5( 1 ) = ONE
147 VM5( 2 ) = EPS
148 VM5( 3 ) = SQRT( SMLNUM )
149 *
150 * Initalization
151 *
152 KNT = 0
153 RMAX = ZERO
154 NINFO = 0
155 SMLNUM = SMLNUM / EPS
156 *
157 * Begin test loop
158 *
159 DO 140 IVM5 = 1, 3
160 DO 130 IVM4 = 1, 5
161 DO 120 IVM3 = 1, 5
162 DO 110 IVM2 = 1, 5
163 DO 100 IVM1 = 1, 5
164 DO 90 NDIM = 1, 6
165 *
166 N = IDIM( NDIM )
167 DO 20 I = 1, N
168 DO 10 J = 1, N
169 T( I, J ) = DBLE( IVAL( I, J, NDIM ) )*
170 $ VM1( IVM1 )
171 IF( I.GE.J )
172 $ T( I, J ) = T( I, J )*VM5( IVM5 )
173 10 CONTINUE
174 20 CONTINUE
175 *
176 W = ONE*VM2( IVM2 )
177 *
178 DO 30 I = 1, N
179 B( I ) = COS( DBLE( I ) )*VM3( IVM3 )
180 30 CONTINUE
181 *
182 DO 40 I = 1, 2*N
183 D( I ) = SIN( DBLE( I ) )*VM4( IVM4 )
184 40 CONTINUE
185 *
186 NORM = DLANGE( '1', N, N, T, LDT, WORK )
187 K = IDAMAX( N, B, 1 )
188 NORMTB = NORM + ABS( B( K ) ) + ABS( W )
189 *
190 CALL DCOPY( N, D, 1, X, 1 )
191 KNT = KNT + 1
192 CALL DLAQTR( .FALSE., .TRUE., N, T, LDT, DUM,
193 $ DUMM, SCALE, X, WORK, INFO )
194 IF( INFO.NE.0 )
195 $ NINFO = NINFO + 1
196 *
197 * || T*x - scale*d || /
198 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
199 *
200 CALL DCOPY( N, D, 1, Y, 1 )
201 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
202 $ X, 1, -SCALE, Y, 1 )
203 XNORM = DASUM( N, X, 1 )
204 RESID = DASUM( N, Y, 1 )
205 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
206 $ ( NORM*EPS )*XNORM )
207 RESID = RESID / DOMIN
208 IF( RESID.GT.RMAX ) THEN
209 RMAX = RESID
210 LMAX = KNT
211 END IF
212 *
213 CALL DCOPY( N, D, 1, X, 1 )
214 KNT = KNT + 1
215 CALL DLAQTR( .TRUE., .TRUE., N, T, LDT, DUM,
216 $ DUMM, SCALE, X, WORK, INFO )
217 IF( INFO.NE.0 )
218 $ NINFO = NINFO + 1
219 *
220 * || T*x - scale*d || /
221 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
222 *
223 CALL DCOPY( N, D, 1, Y, 1 )
224 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
225 $ 1, -SCALE, Y, 1 )
226 XNORM = DASUM( N, X, 1 )
227 RESID = DASUM( N, Y, 1 )
228 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORM,
229 $ ( NORM*EPS )*XNORM )
230 RESID = RESID / DOMIN
231 IF( RESID.GT.RMAX ) THEN
232 RMAX = RESID
233 LMAX = KNT
234 END IF
235 *
236 CALL DCOPY( 2*N, D, 1, X, 1 )
237 KNT = KNT + 1
238 CALL DLAQTR( .FALSE., .FALSE., N, T, LDT, B, W,
239 $ SCALE, X, WORK, INFO )
240 IF( INFO.NE.0 )
241 $ NINFO = NINFO + 1
242 *
243 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
244 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
245 * smlnum/ulp * (||T||+||B||), smlnum )
246 *
247 *
248 CALL DCOPY( 2*N, D, 1, Y, 1 )
249 Y( 1 ) = DDOT( N, B, 1, X( 1+N ), 1 ) +
250 $ SCALE*Y( 1 )
251 DO 50 I = 2, N
252 Y( I ) = W*X( I+N ) + SCALE*Y( I )
253 50 CONTINUE
254 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
255 $ X, 1, -ONE, Y, 1 )
256 *
257 Y( 1+N ) = DDOT( N, B, 1, X, 1 ) -
258 $ SCALE*Y( 1+N )
259 DO 60 I = 2, N
260 Y( I+N ) = W*X( I ) - SCALE*Y( I+N )
261 60 CONTINUE
262 CALL DGEMV( 'No transpose', N, N, ONE, T, LDT,
263 $ X( 1+N ), 1, ONE, Y( 1+N ), 1 )
264 *
265 RESID = DASUM( 2*N, Y, 1 )
266 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
267 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
268 RESID = RESID / DOMIN
269 IF( RESID.GT.RMAX ) THEN
270 RMAX = RESID
271 LMAX = KNT
272 END IF
273 *
274 CALL DCOPY( 2*N, D, 1, X, 1 )
275 KNT = KNT + 1
276 CALL DLAQTR( .TRUE., .FALSE., N, T, LDT, B, W,
277 $ SCALE, X, WORK, INFO )
278 IF( INFO.NE.0 )
279 $ NINFO = NINFO + 1
280 *
281 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
282 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
283 * smlnum/ulp * (||T||+||B||), smlnum )
284 *
285 CALL DCOPY( 2*N, D, 1, Y, 1 )
286 Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 )
287 DO 70 I = 2, N
288 Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) -
289 $ SCALE*Y( I )
290 70 CONTINUE
291 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT, X,
292 $ 1, ONE, Y, 1 )
293 *
294 Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N )
295 DO 80 I = 2, N
296 Y( I+N ) = B( I )*X( 1 ) + W*X( I ) +
297 $ SCALE*Y( I+N )
298 80 CONTINUE
299 CALL DGEMV( 'Transpose', N, N, ONE, T, LDT,
300 $ X( 1+N ), 1, -ONE, Y( 1+N ), 1 )
301 *
302 RESID = DASUM( 2*N, Y, 1 )
303 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
304 $ EPS*( NORMTB*DASUM( 2*N, X, 1 ) ) )
305 RESID = RESID / DOMIN
306 IF( RESID.GT.RMAX ) THEN
307 RMAX = RESID
308 LMAX = KNT
309 END IF
310 *
311 90 CONTINUE
312 100 CONTINUE
313 110 CONTINUE
314 120 CONTINUE
315 130 CONTINUE
316 140 CONTINUE
317 *
318 RETURN
319 *
320 * End of DGET39
321 *
322 END