101 DOUBLE PRECISION zero, half, one
102 parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
103 DOUBLE PRECISION two, three
104 parameter ( two = 2.0d0, three = 3.0d0 )
106 parameter ( lwork = 32 )
109 INTEGER i, ia, ia11, ia12, ia21, ia22, iam, ib, ic,
110 $ ic11, ic12, ic21, ic22, icm, info, j
111 DOUBLE PRECISION bignum, eps, res, smlnum, tnrm
114 DOUBLE PRECISION q( 4, 4 ), result( 2 ), t( 4, 4 ), t1( 4, 4 ),
115 $ val( 9 ), vm( 2 ), work( lwork )
125 INTRINSIC abs, dble, max, sign, sqrt
132 smlnum =
dlamch(
'S' ) / eps
133 bignum = one / smlnum
134 CALL dlabad( smlnum, bignum )
139 val( 2 ) = sqrt( smlnum )
142 val( 5 ) = sqrt( bignum )
143 val( 6 ) = -sqrt( smlnum )
146 val( 9 ) = -sqrt( bignum )
148 vm( 2 ) = one + two*eps
149 CALL dcopy( 16, val( 4 ), 0, t( 1, 1 ), 1 )
163 t( 1, 1 ) = val( ia )*vm( iam )
164 t( 2, 2 ) = val( ic )
165 t( 1, 2 ) = val( ib )
167 tnrm = max( abs( t( 1, 1 ) ), abs( t( 2, 2 ) ),
169 CALL dcopy( 16, t, 1, t1, 1 )
170 CALL dcopy( 16, val( 1 ), 0, q, 1 )
171 CALL dcopy( 4, val( 3 ), 0, q, 5 )
172 CALL dlaexc( .true., 2, t, 4, q, 4, 1, 1, 1, work,
175 $ ninfo( info ) = ninfo( info ) + 1
176 CALL dhst01( 2, 1, 2, t1, 4, t, 4, q, 4, work, lwork,
178 res = result( 1 ) + result( 2 )
180 $ res = res + one / eps
181 IF( t( 1, 1 ).NE.t1( 2, 2 ) )
182 $ res = res + one / eps
183 IF( t( 2, 2 ).NE.t1( 1, 1 ) )
184 $ res = res + one / eps
185 IF( t( 2, 1 ).NE.zero )
186 $ res = res + one / eps
188 IF( res.GT.rmax )
THEN
203 DO 50 ic22 = -1, 1, 2
204 t( 1, 1 ) = val( ia )*vm( iam )
205 t( 1, 2 ) = val( ib )
206 t( 1, 3 ) = -two*val( ib )
208 t( 2, 2 ) = val( ic11 )
209 t( 2, 3 ) = val( ic12 )
211 t( 3, 2 ) = -val( ic21 )
212 t( 3, 3 ) = val( ic11 )*dble( ic22 )
213 tnrm = max( abs( t( 1, 1 ) ),
214 $ abs( t( 1, 2 ) ), abs( t( 1, 3 ) ),
215 $ abs( t( 2, 2 ) ), abs( t( 2, 3 ) ),
216 $ abs( t( 3, 2 ) ), abs( t( 3, 3 ) ) )
217 CALL dcopy( 16, t, 1, t1, 1 )
218 CALL dcopy( 16, val( 1 ), 0, q, 1 )
219 CALL dcopy( 4, val( 3 ), 0, q, 5 )
220 CALL dlaexc( .true., 3, t, 4, q, 4, 1, 1, 2,
223 $ ninfo( info ) = ninfo( info ) + 1
224 CALL dhst01( 3, 1, 3, t1, 4, t, 4, q, 4,
225 $ work, lwork, result )
226 res = result( 1 ) + result( 2 )
228 IF( t1( 1, 1 ).NE.t( 3, 3 ) )
229 $ res = res + one / eps
230 IF( t( 3, 1 ).NE.zero )
231 $ res = res + one / eps
232 IF( t( 3, 2 ).NE.zero )
233 $ res = res + one / eps
234 IF( t( 2, 1 ).NE.0 .AND.
235 $ ( t( 1, 1 ).NE.t( 2,
236 $ 2 ) .OR. sign( one, t( 1,
237 $ 2 ) ).EQ.sign( one, t( 2, 1 ) ) ) )
238 $ res = res + one / eps
241 IF( res.GT.rmax )
THEN
256 DO 150 ia22 = -1, 1, 2
260 t( 1, 1 ) = val( ia11 )
261 t( 1, 2 ) = val( ia12 )
262 t( 1, 3 ) = -two*val( ib )
263 t( 2, 1 ) = -val( ia21 )
264 t( 2, 2 ) = val( ia11 )*dble( ia22 )
265 t( 2, 3 ) = val( ib )
268 t( 3, 3 ) = val( ic )*vm( icm )
269 tnrm = max( abs( t( 1, 1 ) ),
270 $ abs( t( 1, 2 ) ), abs( t( 1, 3 ) ),
271 $ abs( t( 2, 2 ) ), abs( t( 2, 3 ) ),
272 $ abs( t( 3, 2 ) ), abs( t( 3, 3 ) ) )
273 CALL dcopy( 16, t, 1, t1, 1 )
274 CALL dcopy( 16, val( 1 ), 0, q, 1 )
275 CALL dcopy( 4, val( 3 ), 0, q, 5 )
276 CALL dlaexc( .true., 3, t, 4, q, 4, 1, 2, 1,
279 $ ninfo( info ) = ninfo( info ) + 1
280 CALL dhst01( 3, 1, 3, t1, 4, t, 4, q, 4,
281 $ work, lwork, result )
282 res = result( 1 ) + result( 2 )
284 IF( t1( 3, 3 ).NE.t( 1, 1 ) )
285 $ res = res + one / eps
286 IF( t( 2, 1 ).NE.zero )
287 $ res = res + one / eps
288 IF( t( 3, 1 ).NE.zero )
289 $ res = res + one / eps
290 IF( t( 3, 2 ).NE.0 .AND.
291 $ ( t( 2, 2 ).NE.t( 3,
292 $ 3 ) .OR. sign( one, t( 2,
293 $ 3 ) ).EQ.sign( one, t( 3, 2 ) ) ) )
294 $ res = res + one / eps
297 IF( res.GT.rmax )
THEN
312 DO 270 ia22 = -1, 1, 2
317 DO 220 ic22 = -1, 1, 2
320 t( 1, 1 ) = val( ia11 )*vm( iam )
321 t( 1, 2 ) = val( ia12 )*vm( iam )
322 t( 1, 3 ) = -two*val( ib )
323 t( 1, 4 ) = half*val( ib )
324 t( 2, 1 ) = -t( 1, 2 )*val( ia21 )
325 t( 2, 2 ) = val( ia11 )*
326 $ dble( ia22 )*vm( iam )
327 t( 2, 3 ) = val( ib )
328 t( 2, 4 ) = three*val( ib )
331 t( 3, 3 ) = val( ic11 )*
333 t( 3, 4 ) = val( ic12 )*
337 t( 4, 3 ) = -t( 3, 4 )*val( ic21 )*
339 t( 4, 4 ) = val( ic11 )*
349 CALL dcopy( 16, t, 1, t1, 1 )
350 CALL dcopy( 16, val( 1 ), 0, q, 1 )
351 CALL dcopy( 4, val( 3 ), 0, q, 5 )
352 CALL dlaexc( .true., 4, t, 4, q, 4,
353 $ 1, 2, 2, work, info )
355 $ ninfo( info ) = ninfo( info ) + 1
356 CALL dhst01( 4, 1, 4, t1, 4, t, 4,
359 res = result( 1 ) + result( 2 )
361 IF( t( 3, 1 ).NE.zero )
362 $ res = res + one / eps
363 IF( t( 4, 1 ).NE.zero )
364 $ res = res + one / eps
365 IF( t( 3, 2 ).NE.zero )
366 $ res = res + one / eps
367 IF( t( 4, 2 ).NE.zero )
368 $ res = res + one / eps
369 IF( t( 2, 1 ).NE.0 .AND.
370 $ ( t( 1, 1 ).NE.t( 2,
371 $ 2 ) .OR. sign( one, t( 1,
372 $ 2 ) ).EQ.sign( one, t( 2,
373 $ 1 ) ) ) )res = res +
375 IF( t( 4, 3 ).NE.0 .AND.
376 $ ( t( 3, 3 ).NE.t( 4,
377 $ 4 ) .OR. sign( one, t( 3,
378 $ 4 ) ).EQ.sign( one, t( 4,
379 $ 3 ) ) ) )res = res +
383 IF( res.GT.rmax )
THEN
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
double precision function dlamch(CMACH)
DLAMCH
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
subroutine dlabad(SMALL, LARGE)
DLABAD
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...