91 INTEGER knt, lmax, ninfo
99 parameter ( zero = 0.0e0, one = 1.0e0 )
100 REAL two, four, eight
101 parameter ( two = 2.0e0, four = 4.0e0, eight = 8.0e0 )
104 LOGICAL ltranl, ltranr
105 INTEGER ib, ib1, ib2, ib3, info, isgn, itl, itlscl,
106 $ itr, itranl, itranr, itrscl, n1, n2
107 REAL bignum, den, eps, res, scale, sgn, smlnum, tmp,
111 INTEGER itval( 2, 2, 8 )
112 REAL 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 =
slamch(
'S' ) / eps
136 bignum = one / smlnum
137 CALL slabad( 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 slasy2( 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 slasy2( 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 slasy2( 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 slasy2( 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
subroutine slasy2(LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO)
SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
subroutine slabad(SMALL, LARGE)
SLABAD
real function slamch(CMACH)
SLAMCH