82
83
84
85
86
87
88 INTEGER KNT, LMAX
89 DOUBLE PRECISION RMAX
90
91
92 INTEGER NINFO( 2 )
93
94
95
96
97
98 DOUBLE PRECISION ZERO, HALF, ONE
99 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
100 DOUBLE PRECISION TWO, THREE
101 parameter( two = 2.0d0, three = 3.0d0 )
102 INTEGER LWORK
103 parameter( lwork = 32 )
104
105
106 INTEGER I, IA, IA11, IA12, IA21, IA22, IAM, IB, IC,
107 $ IC11, IC12, IC21, IC22, ICM, INFO, J
108 DOUBLE PRECISION BIGNUM, EPS, RES, SMLNUM, TNRM
109
110
111 DOUBLE PRECISION Q( 4, 4 ), RESULT( 2 ), T( 4, 4 ), T1( 4, 4 ),
112 $ VAL( 9 ), VM( 2 ), WORK( LWORK )
113
114
115 DOUBLE PRECISION DLAMCH
117
118
120
121
122 INTRINSIC abs, dble, max, sign, sqrt
123
124
125
126
127
129 smlnum =
dlamch(
'S' ) / eps
130 bignum = one / smlnum
131 CALL dlabad( smlnum, bignum )
132
133
134
135 val( 1 ) = zero
136 val( 2 ) = sqrt( smlnum )
137 val( 3 ) = one
138 val( 4 ) = two
139 val( 5 ) = sqrt( bignum )
140 val( 6 ) = -sqrt( smlnum )
141 val( 7 ) = -one
142 val( 8 ) = -two
143 val( 9 ) = -sqrt( bignum )
144 vm( 1 ) = one
145 vm( 2 ) = one + two*eps
146 CALL dcopy( 16, val( 4 ), 0, t( 1, 1 ), 1 )
147
148 ninfo( 1 ) = 0
149 ninfo( 2 ) = 0
150 knt = 0
151 lmax = 0
152 rmax = zero
153
154
155
156 DO 40 ia = 1, 9
157 DO 30 iam = 1, 2
158 DO 20 ib = 1, 9
159 DO 10 ic = 1, 9
160 t( 1, 1 ) = val( ia )*vm( iam )
161 t( 2, 2 ) = val( ic )
162 t( 1, 2 ) = val( ib )
163 t( 2, 1 ) = zero
164 tnrm = max( abs( t( 1, 1 ) ), abs( t( 2, 2 ) ),
165 $ abs( t( 1, 2 ) ) )
166 CALL dcopy( 16, t, 1, t1, 1 )
167 CALL dcopy( 16, val( 1 ), 0, q, 1 )
168 CALL dcopy( 4, val( 3 ), 0, q, 5 )
169 CALL dlaexc( .true., 2, t, 4, q, 4, 1, 1, 1, work,
170 $ info )
171 IF( info.NE.0 )
172 $ ninfo( info ) = ninfo( info ) + 1
173 CALL dhst01( 2, 1, 2, t1, 4, t, 4, q, 4, work, lwork,
174 $ result )
175 res = result( 1 ) + result( 2 )
176 IF( info.NE.0 )
177 $ res = res + one / eps
178 IF( t( 1, 1 ).NE.t1( 2, 2 ) )
179 $ res = res + one / eps
180 IF( t( 2, 2 ).NE.t1( 1, 1 ) )
181 $ res = res + one / eps
182 IF( t( 2, 1 ).NE.zero )
183 $ res = res + one / eps
184 knt = knt + 1
185 IF( res.GT.rmax ) THEN
186 lmax = knt
187 rmax = res
188 END IF
189 10 CONTINUE
190 20 CONTINUE
191 30 CONTINUE
192 40 CONTINUE
193
194 DO 110 ia = 1, 5
195 DO 100 iam = 1, 2
196 DO 90 ib = 1, 5
197 DO 80 ic11 = 1, 5
198 DO 70 ic12 = 2, 5
199 DO 60 ic21 = 2, 4
200 DO 50 ic22 = -1, 1, 2
201 t( 1, 1 ) = val( ia )*vm( iam )
202 t( 1, 2 ) = val( ib )
203 t( 1, 3 ) = -two*val( ib )
204 t( 2, 1 ) = zero
205 t( 2, 2 ) = val( ic11 )
206 t( 2, 3 ) = val( ic12 )
207 t( 3, 1 ) = zero
208 t( 3, 2 ) = -val( ic21 )
209 t( 3, 3 ) = val( ic11 )*dble( ic22 )
210 tnrm = max( abs( t( 1, 1 ) ),
211 $ abs( t( 1, 2 ) ), abs( t( 1, 3 ) ),
212 $ abs( t( 2, 2 ) ), abs( t( 2, 3 ) ),
213 $ abs( t( 3, 2 ) ), abs( t( 3, 3 ) ) )
214 CALL dcopy( 16, t, 1, t1, 1 )
215 CALL dcopy( 16, val( 1 ), 0, q, 1 )
216 CALL dcopy( 4, val( 3 ), 0, q, 5 )
217 CALL dlaexc( .true., 3, t, 4, q, 4, 1, 1, 2,
218 $ work, info )
219 IF( info.NE.0 )
220 $ ninfo( info ) = ninfo( info ) + 1
221 CALL dhst01( 3, 1, 3, t1, 4, t, 4, q, 4,
222 $ work, lwork, result )
223 res = result( 1 ) + result( 2 )
224 IF( info.EQ.0 ) THEN
225 IF( t1( 1, 1 ).NE.t( 3, 3 ) )
226 $ res = res + one / eps
227 IF( t( 3, 1 ).NE.zero )
228 $ res = res + one / eps
229 IF( t( 3, 2 ).NE.zero )
230 $ res = res + one / eps
231 IF( t( 2, 1 ).NE.0 .AND.
232 $ ( t( 1, 1 ).NE.t( 2,
233 $ 2 ) .OR. sign( one, t( 1,
234 $ 2 ) ).EQ.sign( one, t( 2, 1 ) ) ) )
235 $ res = res + one / eps
236 END IF
237 knt = knt + 1
238 IF( res.GT.rmax ) THEN
239 lmax = knt
240 rmax = res
241 END IF
242 50 CONTINUE
243 60 CONTINUE
244 70 CONTINUE
245 80 CONTINUE
246 90 CONTINUE
247 100 CONTINUE
248 110 CONTINUE
249
250 DO 180 ia11 = 1, 5
251 DO 170 ia12 = 2, 5
252 DO 160 ia21 = 2, 4
253 DO 150 ia22 = -1, 1, 2
254 DO 140 icm = 1, 2
255 DO 130 ib = 1, 5
256 DO 120 ic = 1, 5
257 t( 1, 1 ) = val( ia11 )
258 t( 1, 2 ) = val( ia12 )
259 t( 1, 3 ) = -two*val( ib )
260 t( 2, 1 ) = -val( ia21 )
261 t( 2, 2 ) = val( ia11 )*dble( ia22 )
262 t( 2, 3 ) = val( ib )
263 t( 3, 1 ) = zero
264 t( 3, 2 ) = zero
265 t( 3, 3 ) = val( ic )*vm( icm )
266 tnrm = max( abs( t( 1, 1 ) ),
267 $ abs( t( 1, 2 ) ), abs( t( 1, 3 ) ),
268 $ abs( t( 2, 2 ) ), abs( t( 2, 3 ) ),
269 $ abs( t( 3, 2 ) ), abs( t( 3, 3 ) ) )
270 CALL dcopy( 16, t, 1, t1, 1 )
271 CALL dcopy( 16, val( 1 ), 0, q, 1 )
272 CALL dcopy( 4, val( 3 ), 0, q, 5 )
273 CALL dlaexc( .true., 3, t, 4, q, 4, 1, 2, 1,
274 $ work, info )
275 IF( info.NE.0 )
276 $ ninfo( info ) = ninfo( info ) + 1
277 CALL dhst01( 3, 1, 3, t1, 4, t, 4, q, 4,
278 $ work, lwork, result )
279 res = result( 1 ) + result( 2 )
280 IF( info.EQ.0 ) THEN
281 IF( t1( 3, 3 ).NE.t( 1, 1 ) )
282 $ res = res + one / eps
283 IF( t( 2, 1 ).NE.zero )
284 $ res = res + one / eps
285 IF( t( 3, 1 ).NE.zero )
286 $ res = res + one / eps
287 IF( t( 3, 2 ).NE.0 .AND.
288 $ ( t( 2, 2 ).NE.t( 3,
289 $ 3 ) .OR. sign( one, t( 2,
290 $ 3 ) ).EQ.sign( one, t( 3, 2 ) ) ) )
291 $ res = res + one / eps
292 END IF
293 knt = knt + 1
294 IF( res.GT.rmax ) THEN
295 lmax = knt
296 rmax = res
297 END IF
298 120 CONTINUE
299 130 CONTINUE
300 140 CONTINUE
301 150 CONTINUE
302 160 CONTINUE
303 170 CONTINUE
304 180 CONTINUE
305
306 DO 300 ia11 = 1, 5
307 DO 290 ia12 = 2, 5
308 DO 280 ia21 = 2, 4
309 DO 270 ia22 = -1, 1, 2
310 DO 260 ib = 1, 5
311 DO 250 ic11 = 3, 4
312 DO 240 ic12 = 3, 4
313 DO 230 ic21 = 3, 4
314 DO 220 ic22 = -1, 1, 2
315 DO 210 icm = 5, 7
316 iam = 1
317 t( 1, 1 ) = val( ia11 )*vm( iam )
318 t( 1, 2 ) = val( ia12 )*vm( iam )
319 t( 1, 3 ) = -two*val( ib )
320 t( 1, 4 ) = half*val( ib )
321 t( 2, 1 ) = -t( 1, 2 )*val( ia21 )
322 t( 2, 2 ) = val( ia11 )*
323 $ dble( ia22 )*vm( iam )
324 t( 2, 3 ) = val( ib )
325 t( 2, 4 ) = three*val( ib )
326 t( 3, 1 ) = zero
327 t( 3, 2 ) = zero
328 t( 3, 3 ) = val( ic11 )*
329 $ abs( val( icm ) )
330 t( 3, 4 ) = val( ic12 )*
331 $ abs( val( icm ) )
332 t( 4, 1 ) = zero
333 t( 4, 2 ) = zero
334 t( 4, 3 ) = -t( 3, 4 )*val( ic21 )*
335 $ abs( val( icm ) )
336 t( 4, 4 ) = val( ic11 )*
337 $ dble( ic22 )*
338 $ abs( val( icm ) )
339 tnrm = zero
340 DO 200 i = 1, 4
341 DO 190 j = 1, 4
342 tnrm = max( tnrm,
343 $ abs( t( i, j ) ) )
344 190 CONTINUE
345 200 CONTINUE
346 CALL dcopy( 16, t, 1, t1, 1 )
347 CALL dcopy( 16, val( 1 ), 0, q, 1 )
348 CALL dcopy( 4, val( 3 ), 0, q, 5 )
349 CALL dlaexc( .true., 4, t, 4, q, 4,
350 $ 1, 2, 2, work, info )
351 IF( info.NE.0 )
352 $ ninfo( info ) = ninfo( info ) + 1
353 CALL dhst01( 4, 1, 4, t1, 4, t, 4,
354 $ q, 4, work, lwork,
355 $ result )
356 res = result( 1 ) + result( 2 )
357 IF( info.EQ.0 ) THEN
358 IF( t( 3, 1 ).NE.zero )
359 $ res = res + one / eps
360 IF( t( 4, 1 ).NE.zero )
361 $ res = res + one / eps
362 IF( t( 3, 2 ).NE.zero )
363 $ res = res + one / eps
364 IF( t( 4, 2 ).NE.zero )
365 $ res = res + one / eps
366 IF( t( 2, 1 ).NE.0 .AND.
367 $ ( t( 1, 1 ).NE.t( 2,
368 $ 2 ) .OR. sign( one, t( 1,
369 $ 2 ) ).EQ.sign( one, t( 2,
370 $ 1 ) ) ) )res = res +
371 $ one / eps
372 IF( t( 4, 3 ).NE.0 .AND.
373 $ ( t( 3, 3 ).NE.t( 4,
374 $ 4 ) .OR. sign( one, t( 3,
375 $ 4 ) ).EQ.sign( one, t( 4,
376 $ 3 ) ) ) )res = res +
377 $ one / eps
378 END IF
379 knt = knt + 1
380 IF( res.GT.rmax ) THEN
381 lmax = knt
382 rmax = res
383 END IF
384 210 CONTINUE
385 220 CONTINUE
386 230 CONTINUE
387 240 CONTINUE
388 250 CONTINUE
389 260 CONTINUE
390 270 CONTINUE
391 280 CONTINUE
392 290 CONTINUE
393 300 CONTINUE
394
395 RETURN
396
397
398
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
subroutine dlaexc(WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, INFO)
DLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form...