101 DOUBLE PRECISION rmax
110 DOUBLE PRECISION zero, half, one
111 parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
112 DOUBLE PRECISION two, three, four
113 parameter ( two = 2.0d0, three = 3.0d0, four = 4.0d0 )
114 DOUBLE PRECISION seven, ten
115 parameter ( seven = 7.0d0, ten = 10.0d0 )
116 DOUBLE PRECISION twnone
117 parameter ( twnone = 21.0d0 )
120 INTEGER ia, ib, ica, id1, id2, info, ismin, itrans,
122 DOUBLE PRECISION bignum, ca, d1, d2, den, eps, res, scale, smin,
123 $ smlnum, tmp, unfl, wi, wr, xnorm
126 LOGICAL ltrans( 0: 1 )
127 DOUBLE PRECISION a( 2, 2 ), b( 2, 2 ), vab( 3 ), vca( 5 ),
128 $ vdd( 4 ), vsmin( 4 ), vwi( 4 ), vwr( 4 ),
139 INTRINSIC abs, max, sqrt
142 DATA ltrans / .false., .true. /
150 smlnum =
dlamch(
'S' ) / eps
151 bignum = one / smlnum
152 CALL dlabad( smlnum, bignum )
158 vsmin( 3 ) = one / ( ten*ten )
159 vsmin( 4 ) = one / eps
160 vab( 1 ) = sqrt( smlnum )
162 vab( 3 ) = sqrt( bignum )
171 vdd( 1 ) = sqrt( smlnum )
174 vdd( 4 ) = sqrt( bignum )
176 vca( 2 ) = sqrt( smlnum )
197 smin = vsmin( ismin )
202 a( 1, 1 ) = vab( ia )
204 b( 1, 1 ) = vab( ib )
206 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
208 wr = vwr( iwr )*a( 1, 1 )
213 CALL dlaln2( ltrans( itrans ), na, nw,
214 $ smin, ca, a, 2, d1, d2, b, 2,
215 $ wr, wi, x, 2, scale, xnorm,
218 $ ninfo( 1 ) = ninfo( 1 ) + 1
220 $ ninfo( 2 ) = ninfo( 2 ) + 1
221 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
222 $ x( 1, 1 )-scale*b( 1, 1 ) )
224 den = max( eps*( abs( ( ca*a( 1,
225 $ 1 )-wr*d1 )*x( 1, 1 ) ) ),
228 den = max( smin*abs( x( 1, 1 ) ),
232 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
233 $ abs( b( 1, 1 ) ).LE.smlnum*
234 $ abs( ca*a( 1, 1 )-wr*d1 ) )res = zero
236 $ res = res + one / eps
237 res = res + abs( xnorm-abs( x( 1, 1 ) ) )
238 $ / max( smlnum, xnorm ) / eps
239 IF( info.NE.0 .AND. info.NE.1 )
240 $ res = res + one / eps
242 IF( res.GT.rmax )
THEN
253 a( 1, 1 ) = vab( ia )
255 b( 1, 1 ) = vab( ib )
256 b( 1, 2 ) = -half*vab( ib )
258 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
260 wr = vwr( iwr )*a( 1, 1 )
265 IF( d1.EQ.one .AND. d2.EQ.one .AND.
267 wi = vwi( iwi )*a( 1, 1 )
271 CALL dlaln2( ltrans( itrans ), na, nw,
272 $ smin, ca, a, 2, d1, d2, b,
273 $ 2, wr, wi, x, 2, scale,
276 $ ninfo( 1 ) = ninfo( 1 ) + 1
278 $ ninfo( 2 ) = ninfo( 2 ) + 1
279 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
280 $ x( 1, 1 )+( wi*d1 )*x( 1, 2 )-
282 res = res + abs( ( -wi*d1 )*x( 1, 1 )+
283 $ ( ca*a( 1, 1 )-wr*d1 )*x( 1, 2 )-
286 den = max( eps*( max( abs( ca*a( 1,
287 $ 1 )-wr*d1 ), abs( d1*wi ) )*
288 $ ( abs( x( 1, 1 ) )+abs( x( 1,
289 $ 2 ) ) ) ), smlnum )
291 den = max( smin*( abs( x( 1,
292 $ 1 ) )+abs( x( 1, 2 ) ) ),
296 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
297 $ abs( x( 1, 2 ) ).LT.unfl .AND.
298 $ abs( b( 1, 1 ) ).LE.smlnum*
299 $ abs( ca*a( 1, 1 )-wr*d1 ) )
302 $ res = res + one / eps
303 res = res + abs( xnorm-
305 $ abs( x( 1, 2 ) ) ) /
306 $ max( smlnum, xnorm ) / eps
307 IF( info.NE.0 .AND. info.NE.1 )
308 $ res = res + one / eps
310 IF( res.GT.rmax )
THEN
322 a( 1, 1 ) = vab( ia )
323 a( 1, 2 ) = -three*vab( ia )
324 a( 2, 1 ) = -seven*vab( ia )
325 a( 2, 2 ) = twnone*vab( ia )
327 b( 1, 1 ) = vab( ib )
328 b( 2, 1 ) = -two*vab( ib )
330 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
332 wr = vwr( iwr )*a( 1, 1 )
337 CALL dlaln2( ltrans( itrans ), na, nw,
338 $ smin, ca, a, 2, d1, d2, b, 2,
339 $ wr, wi, x, 2, scale, xnorm,
342 $ ninfo( 1 ) = ninfo( 1 ) + 1
344 $ ninfo( 2 ) = ninfo( 2 ) + 1
345 IF( itrans.EQ.1 )
THEN
347 a( 1, 2 ) = a( 2, 1 )
350 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
351 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
352 $ x( 2, 1 )-scale*b( 1, 1 ) )
353 res = res + abs( ( ca*a( 2, 1 ) )*
354 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
355 $ x( 2, 1 )-scale*b( 2, 1 ) )
357 den = max( eps*( max( abs( ca*a( 1,
358 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
359 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
360 $ 2 )-wr*d2 ) )*max( abs( x( 1,
361 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
364 den = max( eps*( max( smin / eps,
366 $ 1 )-wr*d1 )+abs( ca*a( 1, 2 ) ),
367 $ abs( ca*a( 2, 1 ) )+abs( ca*a( 2,
368 $ 2 )-wr*d2 ) ) )*max( abs( x( 1,
369 $ 1 ) ), abs( x( 2, 1 ) ) ) ),
373 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
374 $ abs( x( 2, 1 ) ).LT.unfl .AND.
375 $ abs( b( 1, 1 ) )+abs( b( 2, 1 ) ).LE.
376 $ smlnum*( abs( ca*a( 1,
377 $ 1 )-wr*d1 )+abs( ca*a( 1,
378 $ 2 ) )+abs( ca*a( 2,
379 $ 1 ) )+abs( ca*a( 2, 2 )-wr*d2 ) ) )
382 $ res = res + one / eps
383 res = res + abs( xnorm-
384 $ max( abs( x( 1, 1 ) ), abs( x( 2,
385 $ 1 ) ) ) ) / max( smlnum, xnorm ) /
387 IF( info.NE.0 .AND. info.NE.1 )
388 $ res = res + one / eps
390 IF( res.GT.rmax )
THEN
401 a( 1, 1 ) = vab( ia )*two
402 a( 1, 2 ) = -three*vab( ia )
403 a( 2, 1 ) = -seven*vab( ia )
404 a( 2, 2 ) = twnone*vab( ia )
406 b( 1, 1 ) = vab( ib )
407 b( 2, 1 ) = -two*vab( ib )
408 b( 1, 2 ) = four*vab( ib )
409 b( 2, 2 ) = -seven*vab( ib )
411 IF( d1.EQ.one .AND. d2.EQ.one .AND. ca.EQ.
413 wr = vwr( iwr )*a( 1, 1 )
418 IF( d1.EQ.one .AND. d2.EQ.one .AND.
420 wi = vwi( iwi )*a( 1, 1 )
424 CALL dlaln2( ltrans( itrans ), na, nw,
425 $ smin, ca, a, 2, d1, d2, b,
426 $ 2, wr, wi, x, 2, scale,
429 $ ninfo( 1 ) = ninfo( 1 ) + 1
431 $ ninfo( 2 ) = ninfo( 2 ) + 1
432 IF( itrans.EQ.1 )
THEN
434 a( 1, 2 ) = a( 2, 1 )
437 res = abs( ( ca*a( 1, 1 )-wr*d1 )*
438 $ x( 1, 1 )+( ca*a( 1, 2 ) )*
439 $ x( 2, 1 )+( wi*d1 )*x( 1, 2 )-
441 res = res + abs( ( ca*a( 1,
442 $ 1 )-wr*d1 )*x( 1, 2 )+
443 $ ( ca*a( 1, 2 ) )*x( 2, 2 )-
444 $ ( wi*d1 )*x( 1, 1 )-scale*
446 res = res + abs( ( ca*a( 2, 1 ) )*
447 $ x( 1, 1 )+( ca*a( 2, 2 )-wr*d2 )*
448 $ x( 2, 1 )+( wi*d2 )*x( 2, 2 )-
450 res = res + abs( ( ca*a( 2, 1 ) )*
451 $ x( 1, 2 )+( ca*a( 2, 2 )-wr*d2 )*
452 $ x( 2, 2 )-( wi*d2 )*x( 2, 1 )-
455 den = max( eps*( max( abs( ca*a( 1,
456 $ 1 )-wr*d1 )+abs( ca*a( 1,
457 $ 2 ) )+abs( wi*d1 ),
459 $ 1 ) )+abs( ca*a( 2,
460 $ 2 )-wr*d2 )+abs( wi*d2 ) )*
462 $ 1 ) )+abs( x( 2, 1 ) ),
463 $ abs( x( 1, 2 ) )+abs( x( 2,
464 $ 2 ) ) ) ), smlnum )
466 den = max( eps*( max( smin / eps,
468 $ 1 )-wr*d1 )+abs( ca*a( 1,
469 $ 2 ) )+abs( wi*d1 ),
471 $ 1 ) )+abs( ca*a( 2,
472 $ 2 )-wr*d2 )+abs( wi*d2 ) ) )*
474 $ 1 ) )+abs( x( 2, 1 ) ),
475 $ abs( x( 1, 2 ) )+abs( x( 2,
476 $ 2 ) ) ) ), smlnum )
479 IF( abs( x( 1, 1 ) ).LT.unfl .AND.
480 $ abs( x( 2, 1 ) ).LT.unfl .AND.
481 $ abs( x( 1, 2 ) ).LT.unfl .AND.
482 $ abs( x( 2, 2 ) ).LT.unfl .AND.
484 $ abs( b( 2, 1 ) ).LE.smlnum*
485 $ ( abs( ca*a( 1, 1 )-wr*d1 )+
486 $ abs( ca*a( 1, 2 ) )+abs( ca*a( 2,
487 $ 1 ) )+abs( ca*a( 2,
488 $ 2 )-wr*d2 )+abs( wi*d2 )+abs( wi*
491 $ res = res + one / eps
492 res = res + abs( xnorm-
493 $ max( abs( x( 1, 1 ) )+abs( x( 1,
495 $ 1 ) )+abs( x( 2, 2 ) ) ) ) /
496 $ max( smlnum, xnorm ) / eps
497 IF( info.NE.0 .AND. info.NE.1 )
498 $ res = res + one / eps
500 IF( res.GT.rmax )
THEN
double precision function dlamch(CMACH)
DLAMCH
subroutine dlaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
subroutine dlabad(SMALL, LARGE)
DLABAD