83 SUBROUTINE dget34( RMAX, LMAX, NINFO, KNT )
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
subroutine dhst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
DHST01
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dget34(RMAX, LMAX, NINFO, KNT)
DGET34
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
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...