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 dlabad(SMALL, LARGE)
DLABAD
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.