176 DOUBLE PRECISION scale, w
179 DOUBLE PRECISION b( * ), t( ldt, * ), work( * ), x( * )
185 DOUBLE PRECISION zero, one
186 parameter ( zero = 0.0d+0, one = 1.0d+0 )
190 INTEGER i, ierr, j, j1, j2, jnext, k, n1, n2
191 DOUBLE PRECISION bignum, eps, rec, scaloc, si, smin, sminw,
192 $ smlnum, sr, tjj, tmp, xj, xmax, xnorm, z
195 DOUBLE PRECISION d( 2, 2 ), v( 2, 2 )
223 smlnum =
dlamch(
'S' ) / eps
224 bignum = one / smlnum
226 xnorm =
dlange(
'M', n, n, t, ldt, d )
228 $ xnorm = max( xnorm, abs( w ),
dlange(
'M', n, 1, b, n, d ) )
229 smin = max( smlnum, eps*xnorm )
236 work( j ) =
dasum( j-1, t( 1, j ), 1 )
239 IF( .NOT.lreal )
THEN
241 work( i ) = work( i ) + abs( b( i ) )
253 IF( xmax.GT.bignum )
THEN
254 scale = bignum / xmax
255 CALL dscal( n1, scale, x, 1 )
273 IF( t( j, j-1 ).NE.zero )
THEN
287 tjj = abs( t( j1, j1 ) )
289 IF( tjj.LT.smin )
THEN
298 IF( tjj.LT.one )
THEN
299 IF( xj.GT.bignum*tjj )
THEN
301 CALL dscal( n, rec, x, 1 )
306 x( j1 ) = x( j1 ) / tmp
314 IF( work( j1 ).GT.( bignum-xmax )*rec )
THEN
315 CALL dscal( n, rec, x, 1 )
320 CALL daxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
334 CALL dlaln2( .false., 2, 1, smin, one, t( j1, j1 ),
335 $ ldt, one, one, d, 2, zero, zero, v, 2,
336 $ scaloc, xnorm, ierr )
340 IF( scaloc.NE.one )
THEN
341 CALL dscal( n, scaloc, x, 1 )
350 xj = max( abs( v( 1, 1 ) ), abs( v( 2, 1 ) ) )
353 IF( max( work( j1 ), work( j2 ) ).GT.
354 $ ( bignum-xmax )*rec )
THEN
355 CALL dscal( n, rec, x, 1 )
363 CALL daxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
364 CALL daxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x, 1 )
385 IF( t( j+1, j ).NE.zero )
THEN
399 IF( xmax.GT.one )
THEN
401 IF( work( j1 ).GT.( bignum-xj )*rec )
THEN
402 CALL dscal( n, rec, x, 1 )
408 x( j1 ) = x( j1 ) -
ddot( j1-1, t( 1, j1 ), 1, x, 1 )
411 tjj = abs( t( j1, j1 ) )
413 IF( tjj.LT.smin )
THEN
419 IF( tjj.LT.one )
THEN
420 IF( xj.GT.bignum*tjj )
THEN
422 CALL dscal( n, rec, x, 1 )
427 x( j1 ) = x( j1 ) / tmp
428 xmax = max( xmax, abs( x( j1 ) ) )
437 xj = max( abs( x( j1 ) ), abs( x( j2 ) ) )
438 IF( xmax.GT.one )
THEN
440 IF( max( work( j2 ), work( j1 ) ).GT.( bignum-xj )*
442 CALL dscal( n, rec, x, 1 )
448 d( 1, 1 ) = x( j1 ) -
ddot( j1-1, t( 1, j1 ), 1, x,
450 d( 2, 1 ) = x( j2 ) -
ddot( j1-1, t( 1, j2 ), 1, x,
453 CALL dlaln2( .true., 2, 1, smin, one, t( j1, j1 ),
454 $ ldt, one, one, d, 2, zero, zero, v, 2,
455 $ scaloc, xnorm, ierr )
459 IF( scaloc.NE.one )
THEN
460 CALL dscal( n, scaloc, x, 1 )
465 xmax = max( abs( x( j1 ) ), abs( x( j2 ) ), xmax )
473 sminw = max( eps*abs( w ), smin )
486 IF( t( j, j-1 ).NE.zero )
THEN
501 xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
502 tjj = abs( t( j1, j1 ) ) + abs( z )
504 IF( tjj.LT.sminw )
THEN
513 IF( tjj.LT.one )
THEN
514 IF( xj.GT.bignum*tjj )
THEN
516 CALL dscal( n2, rec, x, 1 )
521 CALL dladiv( x( j1 ), x( n+j1 ), tmp, z, sr, si )
524 xj = abs( x( j1 ) ) + abs( x( n+j1 ) )
531 IF( work( j1 ).GT.( bignum-xmax )*rec )
THEN
532 CALL dscal( n2, rec, x, 1 )
538 CALL daxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
539 CALL daxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
542 x( 1 ) = x( 1 ) + b( j1 )*x( n+j1 )
543 x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 )
547 xmax = max( xmax, abs( x( k ) )+
558 d( 1, 2 ) = x( n+j1 )
559 d( 2, 2 ) = x( n+j2 )
560 CALL dlaln2( .false., 2, 2, sminw, one, t( j1, j1 ),
561 $ ldt, one, one, d, 2, zero, -w, v, 2,
562 $ scaloc, xnorm, ierr )
566 IF( scaloc.NE.one )
THEN
567 CALL dscal( 2*n, scaloc, x, 1 )
572 x( n+j1 ) = v( 1, 2 )
573 x( n+j2 ) = v( 2, 2 )
578 xj = max( abs( v( 1, 1 ) )+abs( v( 1, 2 ) ),
579 $ abs( v( 2, 1 ) )+abs( v( 2, 2 ) ) )
582 IF( max( work( j1 ), work( j2 ) ).GT.
583 $ ( bignum-xmax )*rec )
THEN
584 CALL dscal( n2, rec, x, 1 )
592 CALL daxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
593 CALL daxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x, 1 )
595 CALL daxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
597 CALL daxpy( j1-1, -x( n+j2 ), t( 1, j2 ), 1,
600 x( 1 ) = x( 1 ) + b( j1 )*x( n+j1 ) +
602 x( n+1 ) = x( n+1 ) - b( j1 )*x( j1 ) -
607 xmax = max( abs( x( k ) )+abs( x( k+n ) ),
627 IF( t( j+1, j ).NE.zero )
THEN
640 xj = abs( x( j1 ) ) + abs( x( j1+n ) )
641 IF( xmax.GT.one )
THEN
643 IF( work( j1 ).GT.( bignum-xj )*rec )
THEN
644 CALL dscal( n2, rec, x, 1 )
650 x( j1 ) = x( j1 ) -
ddot( j1-1, t( 1, j1 ), 1, x, 1 )
651 x( n+j1 ) = x( n+j1 ) -
ddot( j1-1, t( 1, j1 ), 1,
654 x( j1 ) = x( j1 ) - b( j1 )*x( n+1 )
655 x( n+j1 ) = x( n+j1 ) + b( j1 )*x( 1 )
657 xj = abs( x( j1 ) ) + abs( x( j1+n ) )
666 tjj = abs( t( j1, j1 ) ) + abs( z )
668 IF( tjj.LT.sminw )
THEN
674 IF( tjj.LT.one )
THEN
675 IF( xj.GT.bignum*tjj )
THEN
677 CALL dscal( n2, rec, x, 1 )
682 CALL dladiv( x( j1 ), x( n+j1 ), tmp, -z, sr, si )
685 xmax = max( abs( x( j1 ) )+abs( x( j1+n ) ), xmax )
694 xj = max( abs( x( j1 ) )+abs( x( n+j1 ) ),
695 $ abs( x( j2 ) )+abs( x( n+j2 ) ) )
696 IF( xmax.GT.one )
THEN
698 IF( max( work( j1 ), work( j2 ) ).GT.
699 $ ( bignum-xj ) / xmax )
THEN
700 CALL dscal( n2, rec, x, 1 )
706 d( 1, 1 ) = x( j1 ) -
ddot( j1-1, t( 1, j1 ), 1, x,
708 d( 2, 1 ) = x( j2 ) -
ddot( j1-1, t( 1, j2 ), 1, x,
710 d( 1, 2 ) = x( n+j1 ) -
ddot( j1-1, t( 1, j1 ), 1,
712 d( 2, 2 ) = x( n+j2 ) -
ddot( j1-1, t( 1, j2 ), 1,
714 d( 1, 1 ) = d( 1, 1 ) - b( j1 )*x( n+1 )
715 d( 2, 1 ) = d( 2, 1 ) - b( j2 )*x( n+1 )
716 d( 1, 2 ) = d( 1, 2 ) + b( j1 )*x( 1 )
717 d( 2, 2 ) = d( 2, 2 ) + b( j2 )*x( 1 )
719 CALL dlaln2( .true., 2, 2, sminw, one, t( j1, j1 ),
720 $ ldt, one, one, d, 2, zero, w, v, 2,
721 $ scaloc, xnorm, ierr )
725 IF( scaloc.NE.one )
THEN
726 CALL dscal( n2, scaloc, x, 1 )
731 x( n+j1 ) = v( 1, 2 )
732 x( n+j2 ) = v( 2, 2 )
733 xmax = max( abs( x( j1 ) )+abs( x( n+j1 ) ),
734 $ abs( x( j2 ) )+abs( x( n+j2 ) ), xmax )
integer function idamax(N, DX, INCX)
IDAMAX
double precision function dlamch(CMACH)
DLAMCH
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
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.
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dscal(N, DA, DX, INCX)
DSCAL
double precision function dasum(N, DX, INCX)
DASUM