91 INTEGER knt, lmax, ninfo
98 DOUBLE PRECISION zero, one
99 parameter ( zero = 0.0d0, one = 1.0d0 )
100 DOUBLE PRECISION two, four, eight
101 parameter ( two = 2.0d0, four = 4.0d0, eight = 8.0d0 )
104 LOGICAL ltranl, ltranr
105 INTEGER ib, ib1, ib2, ib3, info, isgn, itl, itlscl,
106 $ itr, itranl, itranr, itrscl, n1, n2
107 DOUBLE PRECISION bignum, den, eps, res, scale, sgn, smlnum, tmp,
111 INTEGER itval( 2, 2, 8 )
112 DOUBLE PRECISION b( 2, 2 ), tl( 2, 2 ), tr( 2, 2 ), val( 3 ),
123 INTRINSIC abs, max, min, sqrt
126 DATA itval / 8, 4, 2, 1, 4, 8, 1, 2, 2, 1, 8, 4, 1,
127 $ 2, 4, 8, 9, 4, 2, 1, 4, 9, 1, 2, 2, 1, 9, 4, 1,
135 smlnum =
dlamch(
'S' ) / eps
136 bignum = one / smlnum
137 CALL dlabad( smlnum, bignum )
141 val( 1 ) = sqrt( smlnum )
143 val( 3 ) = sqrt( bignum )
154 DO 210 isgn = -1, 1, 2
164 tl( 1, 1 ) = val( itl )
165 tr( 1, 1 ) = val( itr )
166 b( 1, 1 ) = val( ib )
168 CALL dlasy2( ltranl, ltranr, isgn, n1, n2, tl,
169 $ 2, tr, 2, b, 2, scale, x, 2, xnorm,
173 res = abs( ( tl( 1, 1 )+sgn*tr( 1, 1 ) )*
174 $ x( 1, 1 )-scale*b( 1, 1 ) )
176 den = max( eps*( ( abs( tr( 1,
177 $ 1 ) )+abs( tl( 1, 1 ) ) )*abs( x( 1,
180 den = smlnum*max( abs( x( 1, 1 ) ), one )
184 $ res = res + one / eps
185 res = res + abs( xnorm-abs( x( 1, 1 ) ) ) /
186 $ max( smlnum, xnorm ) / eps
187 IF( info.NE.0 .AND. info.NE.1 )
188 $ res = res + one / eps
189 IF( res.GT.rmax )
THEN
204 b( 1, 1 ) = val( ib1 )
205 b( 2, 1 ) = -four*val( ib2 )
206 tl( 1, 1 ) = itval( 1, 1, itl )*
208 tl( 2, 1 ) = itval( 2, 1, itl )*
210 tl( 1, 2 ) = itval( 1, 2, itl )*
212 tl( 2, 2 ) = itval( 2, 2, itl )*
214 tr( 1, 1 ) = val( itr )
216 CALL dlasy2( ltranl, ltranr, isgn, n1, n2,
217 $ tl, 2, tr, 2, b, 2, scale, x,
223 tl( 1, 2 ) = tl( 2, 1 )
226 res = abs( ( tl( 1, 1 )+sgn*tr( 1, 1 ) )*
227 $ x( 1, 1 )+tl( 1, 2 )*x( 2, 1 )-
229 res = res + abs( ( tl( 2, 2 )+sgn*tr( 1,
230 $ 1 ) )*x( 2, 1 )+tl( 2, 1 )*
231 $ x( 1, 1 )-scale*b( 2, 1 ) )
232 tnrm = abs( tr( 1, 1 ) ) +
233 $ abs( tl( 1, 1 ) ) +
234 $ abs( tl( 1, 2 ) ) +
235 $ abs( tl( 2, 1 ) ) +
237 xnrm = max( abs( x( 1, 1 ) ),
239 den = max( smlnum, smlnum*xnrm,
240 $ ( tnrm*eps )*xnrm )
243 $ res = res + one / eps
244 res = res + abs( xnorm-xnrm ) /
245 $ max( smlnum, xnorm ) / eps
246 IF( res.GT.rmax )
THEN
263 b( 1, 1 ) = val( ib1 )
264 b( 1, 2 ) = -two*val( ib2 )
265 tr( 1, 1 ) = itval( 1, 1, itr )*
267 tr( 2, 1 ) = itval( 2, 1, itr )*
269 tr( 1, 2 ) = itval( 1, 2, itr )*
271 tr( 2, 2 ) = itval( 2, 2, itr )*
273 tl( 1, 1 ) = val( itl )
275 CALL dlasy2( ltranl, ltranr, isgn, n1, n2,
276 $ tl, 2, tr, 2, b, 2, scale, x,
282 tr( 1, 2 ) = tr( 2, 1 )
285 tnrm = abs( tl( 1, 1 ) ) +
286 $ abs( tr( 1, 1 ) ) +
287 $ abs( tr( 1, 2 ) ) +
288 $ abs( tr( 2, 2 ) ) +
290 xnrm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
291 res = abs( ( ( tl( 1, 1 )+sgn*tr( 1,
292 $ 1 ) ) )*( x( 1, 1 ) )+
293 $ ( sgn*tr( 2, 1 ) )*( x( 1, 2 ) )-
294 $ ( scale*b( 1, 1 ) ) )
295 res = res + abs( ( ( tl( 1, 1 )+sgn*tr( 2,
296 $ 2 ) ) )*( x( 1, 2 ) )+
297 $ ( sgn*tr( 1, 2 ) )*( x( 1, 1 ) )-
298 $ ( scale*b( 1, 2 ) ) )
299 den = max( smlnum, smlnum*xnrm,
300 $ ( tnrm*eps )*xnrm )
303 $ res = res + one / eps
304 res = res + abs( xnorm-xnrm ) /
305 $ max( smlnum, xnorm ) / eps
306 IF( res.GT.rmax )
THEN
325 b( 1, 1 ) = val( ib1 )
326 b( 2, 1 ) = -four*val( ib2 )
327 b( 1, 2 ) = -two*val( ib3 )
329 $ min( val( ib1 ), val
330 $ ( ib2 ), val( ib3 ) )
331 tr( 1, 1 ) = itval( 1, 1, itr )*
333 tr( 2, 1 ) = itval( 2, 1, itr )*
335 tr( 1, 2 ) = itval( 1, 2, itr )*
337 tr( 2, 2 ) = itval( 2, 2, itr )*
339 tl( 1, 1 ) = itval( 1, 1, itl )*
341 tl( 2, 1 ) = itval( 2, 1, itl )*
343 tl( 1, 2 ) = itval( 1, 2, itl )*
345 tl( 2, 2 ) = itval( 2, 2, itl )*
348 CALL dlasy2( ltranl, ltranr, isgn,
349 $ n1, n2, tl, 2, tr, 2,
356 tr( 1, 2 ) = tr( 2, 1 )
361 tl( 1, 2 ) = tl( 2, 1 )
364 tnrm = abs( tr( 1, 1 ) ) +
365 $ abs( tr( 2, 1 ) ) +
366 $ abs( tr( 1, 2 ) ) +
367 $ abs( tr( 2, 2 ) ) +
368 $ abs( tl( 1, 1 ) ) +
369 $ abs( tl( 2, 1 ) ) +
370 $ abs( tl( 1, 2 ) ) +
372 xnrm = max( abs( x( 1, 1 ) )+
376 res = abs( ( ( tl( 1, 1 )+sgn*tr( 1,
377 $ 1 ) ) )*( x( 1, 1 ) )+
378 $ ( sgn*tr( 2, 1 ) )*
379 $ ( x( 1, 2 ) )+( tl( 1, 2 ) )*
381 $ ( scale*b( 1, 1 ) ) )
382 res = res + abs( ( tl( 1, 1 ) )*
384 $ ( sgn*tr( 1, 2 ) )*
386 $ ( sgn*tr( 2, 2 ) )*
387 $ ( x( 1, 2 ) )+( tl( 1, 2 ) )*
389 $ ( scale*b( 1, 2 ) ) )
390 res = res + abs( ( tl( 2, 1 ) )*
392 $ ( sgn*tr( 1, 1 ) )*
394 $ ( sgn*tr( 2, 1 ) )*
395 $ ( x( 2, 2 ) )+( tl( 2, 2 ) )*
397 $ ( scale*b( 2, 1 ) ) )
398 res = res + abs( ( ( tl( 2,
399 $ 2 )+sgn*tr( 2, 2 ) ) )*
401 $ ( sgn*tr( 1, 2 ) )*
402 $ ( x( 2, 1 ) )+( tl( 2, 1 ) )*
404 $ ( scale*b( 2, 2 ) ) )
405 den = max( smlnum, smlnum*xnrm,
406 $ ( tnrm*eps )*xnrm )
409 $ res = res + one / eps
410 res = res + abs( xnorm-xnrm ) /
411 $ max( smlnum, xnorm ) / eps
412 IF( res.GT.rmax )
THEN
double precision function dlamch(CMACH)
DLAMCH
subroutine dlasy2(LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO)
DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
subroutine dlabad(SMALL, LARGE)
DLABAD