102 parameter ( zero = 0.0e0, half = 0.5e0, one = 1.0e0 )
104 parameter ( two = 2.0e0, three = 3.0e0 )
106 parameter ( lwork = 32 )
109 INTEGER i, ia, ia11, ia12, ia21, ia22, iam, ib, ic,
110 $ ic11, ic12, ic21, ic22, icm, info, j
111 REAL bignum, eps, res, smlnum, tnrm
114 REAL q( 4, 4 ), result( 2 ), t( 4, 4 ), t1( 4, 4 ),
115 $ val( 9 ), vm( 2 ), work( lwork )
125 INTRINSIC abs, max,
REAL, sign, sqrt
132 smlnum =
slamch(
'S' ) / eps
133 bignum = one / smlnum
134 CALL slabad( 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 scopy( 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 scopy( 16, t, 1, t1, 1 )
170 CALL scopy( 16, val( 1 ), 0, q, 1 )
171 CALL scopy( 4, val( 3 ), 0, q, 5 )
172 CALL slaexc( .true., 2, t, 4, q, 4, 1, 1, 1, work,
175 $ ninfo( info ) = ninfo( info ) + 1
176 CALL shst01( 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 )*
REAL( 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 scopy( 16, t, 1, t1, 1 )
218 CALL scopy( 16, val( 1 ), 0, q, 1 )
219 CALL scopy( 4, val( 3 ), 0, q, 5 )
220 CALL slaexc( .true., 3, t, 4, q, 4, 1, 1, 2,
223 $ ninfo( info ) = ninfo( info ) + 1
224 CALL shst01( 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 )*
REAL( 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 scopy( 16, t, 1, t1, 1 )
274 CALL scopy( 16, val( 1 ), 0, q, 1 )
275 CALL scopy( 4, val( 3 ), 0, q, 5 )
276 CALL slaexc( .true., 3, t, 4, q, 4, 1, 2, 1,
279 $ ninfo( info ) = ninfo( info ) + 1
280 CALL shst01( 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 $
REAL( 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 scopy( 16, t, 1, t1, 1 )
350 CALL scopy( 16, val( 1 ), 0, q, 1 )
351 CALL scopy( 4, val( 3 ), 0, q, 5 )
352 CALL slaexc( .true., 4, t, 4, q, 4,
353 $ 1, 2, 2, work, info )
355 $ ninfo( info ) = ninfo( info ) + 1
356 CALL shst01( 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 slabad(SMALL, LARGE)
SLABAD
subroutine shst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
SHST01
subroutine slaexc(WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK, INFO)
SLAEXC swaps adjacent diagonal blocks of a real upper quasi-triangular matrix in Schur canonical form...
real function slamch(CMACH)
SLAMCH
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY