216 SUBROUTINE dlaln2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
217 $ LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
225 INTEGER INFO, LDA, LDB, LDX, NA, NW
226 DOUBLE PRECISION CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
229 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), X( LDX, * )
235 DOUBLE PRECISION ZERO, ONE
236 parameter( zero = 0.0d0, one = 1.0d0 )
238 parameter( two = 2.0d0 )
242 DOUBLE PRECISION BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
243 $ ci22, cmax, cnorm, cr21, cr22, csi, csr, li21,
244 $ lr21, smini, smlnum, temp, u22abs, ui11, ui11r,
245 $ ui12, ui12s, ui22, ur11, ur11r, ur12, ur12s,
246 $ ur22, xi1, xi2, xr1, xr2
249 LOGICAL RSWAP( 4 ), ZSWAP( 4 )
250 INTEGER IPIVOT( 4, 4 )
251 DOUBLE PRECISION CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
254 DOUBLE PRECISION DLAMCH
264 equivalence( ci( 1, 1 ), civ( 1 ) ),
265 $ ( cr( 1, 1 ), crv( 1 ) )
268 DATA zswap / .false., .false., .true., .true. /
269 DATA rswap / .false., .true., .false., .true. /
270 DATA ipivot / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
277 smlnum = two*dlamch(
'Safe minimum' )
278 bignum = one / smlnum
279 smini = max( smin, smlnum )
299 csr = ca*a( 1, 1 ) - wr*d1
304 IF( cnorm.LT.smini )
THEN
312 bnorm = abs( b( 1, 1 ) )
313 IF( cnorm.LT.one .AND. bnorm.GT.one )
THEN
314 IF( bnorm.GT.bignum*cnorm )
315 $ scale = one / bnorm
320 x( 1, 1 ) = ( b( 1, 1 )*scale ) / csr
321 xnorm = abs( x( 1, 1 ) )
328 csr = ca*a( 1, 1 ) - wr*d1
330 cnorm = abs( csr ) + abs( csi )
334 IF( cnorm.LT.smini )
THEN
343 bnorm = abs( b( 1, 1 ) ) + abs( b( 1, 2 ) )
344 IF( cnorm.LT.one .AND. bnorm.GT.one )
THEN
345 IF( bnorm.GT.bignum*cnorm )
346 $ scale = one / bnorm
351 CALL dladiv( scale*b( 1, 1 ), scale*b( 1, 2 ), csr, csi,
352 $ x( 1, 1 ), x( 1, 2 ) )
353 xnorm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
362 cr( 1, 1 ) = ca*a( 1, 1 ) - wr*d1
363 cr( 2, 2 ) = ca*a( 2, 2 ) - wr*d2
365 cr( 1, 2 ) = ca*a( 2, 1 )
366 cr( 2, 1 ) = ca*a( 1, 2 )
368 cr( 2, 1 ) = ca*a( 2, 1 )
369 cr( 1, 2 ) = ca*a( 1, 2 )
382 IF( abs( crv( j ) ).GT.cmax )
THEN
383 cmax = abs( crv( j ) )
390 IF( cmax.LT.smini )
THEN
391 bnorm = max( abs( b( 1, 1 ) ), abs( b( 2, 1 ) ) )
392 IF( smini.LT.one .AND. bnorm.GT.one )
THEN
393 IF( bnorm.GT.bignum*smini )
394 $ scale = one / bnorm
397 x( 1, 1 ) = temp*b( 1, 1 )
398 x( 2, 1 ) = temp*b( 2, 1 )
407 cr21 = crv( ipivot( 2, icmax ) )
408 ur12 = crv( ipivot( 3, icmax ) )
409 cr22 = crv( ipivot( 4, icmax ) )
412 ur22 = cr22 - ur12*lr21
416 IF( abs( ur22 ).LT.smini )
THEN
420 IF( rswap( icmax ) )
THEN
428 bbnd = max( abs( br1*( ur22*ur11r ) ), abs( br2 ) )
429 IF( bbnd.GT.one .AND. abs( ur22 ).LT.one )
THEN
430 IF( bbnd.GE.bignum*abs( ur22 ) )
434 xr2 = ( br2*scale ) / ur22
435 xr1 = ( scale*br1 )*ur11r - xr2*( ur11r*ur12 )
436 IF( zswap( icmax ) )
THEN
443 xnorm = max( abs( xr1 ), abs( xr2 ) )
447 IF( xnorm.GT.one .AND. cmax.GT.one )
THEN
448 IF( xnorm.GT.bignum / cmax )
THEN
450 x( 1, 1 ) = temp*x( 1, 1 )
451 x( 2, 1 ) = temp*x( 2, 1 )
470 IF( abs( crv( j ) )+abs( civ( j ) ).GT.cmax )
THEN
471 cmax = abs( crv( j ) ) + abs( civ( j ) )
478 IF( cmax.LT.smini )
THEN
479 bnorm = max( abs( b( 1, 1 ) )+abs( b( 1, 2 ) ),
480 $ abs( b( 2, 1 ) )+abs( b( 2, 2 ) ) )
481 IF( smini.LT.one .AND. bnorm.GT.one )
THEN
482 IF( bnorm.GT.bignum*smini )
483 $ scale = one / bnorm
486 x( 1, 1 ) = temp*b( 1, 1 )
487 x( 2, 1 ) = temp*b( 2, 1 )
488 x( 1, 2 ) = temp*b( 1, 2 )
489 x( 2, 2 ) = temp*b( 2, 2 )
499 cr21 = crv( ipivot( 2, icmax ) )
500 ci21 = civ( ipivot( 2, icmax ) )
501 ur12 = crv( ipivot( 3, icmax ) )
502 ui12 = civ( ipivot( 3, icmax ) )
503 cr22 = crv( ipivot( 4, icmax ) )
504 ci22 = civ( ipivot( 4, icmax ) )
505 IF( icmax.EQ.1 .OR. icmax.EQ.4 )
THEN
509 IF( abs( ur11 ).GT.abs( ui11 ) )
THEN
511 ur11r = one / ( ur11*( one+temp**2 ) )
515 ui11r = -one / ( ui11*( one+temp**2 ) )
522 ur22 = cr22 - ur12*lr21
523 ui22 = ci22 - ur12*li21
534 ur22 = cr22 - ur12*lr21 + ui12*li21
535 ui22 = -ur12*li21 - ui12*lr21
537 u22abs = abs( ur22 ) + abs( ui22 )
541 IF( u22abs.LT.smini )
THEN
546 IF( rswap( icmax ) )
THEN
557 br2 = br2 - lr21*br1 + li21*bi1
558 bi2 = bi2 - li21*br1 - lr21*bi1
559 bbnd = max( ( abs( br1 )+abs( bi1 ) )*
560 $ ( u22abs*( abs( ur11r )+abs( ui11r ) ) ),
561 $ abs( br2 )+abs( bi2 ) )
562 IF( bbnd.GT.one .AND. u22abs.LT.one )
THEN
563 IF( bbnd.GE.bignum*u22abs )
THEN
572 CALL dladiv( br2, bi2, ur22, ui22, xr2, xi2 )
573 xr1 = ur11r*br1 - ui11r*bi1 - ur12s*xr2 + ui12s*xi2
574 xi1 = ui11r*br1 + ur11r*bi1 - ui12s*xr2 - ur12s*xi2
575 IF( zswap( icmax ) )
THEN
586 xnorm = max( abs( xr1 )+abs( xi1 ), abs( xr2 )+abs( xi2 ) )
590 IF( xnorm.GT.one .AND. cmax.GT.one )
THEN
591 IF( xnorm.GT.bignum / cmax )
THEN
593 x( 1, 1 ) = temp*x( 1, 1 )
594 x( 2, 1 ) = temp*x( 2, 1 )
595 x( 1, 2 ) = temp*x( 1, 2 )
596 x( 2, 2 ) = temp*x( 2, 2 )
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 dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.