179 REAL b( * ), t( ldt, * ), work( * ), x( * )
186 parameter ( zero = 0.0e+0, one = 1.0e+0 )
190 INTEGER i, ierr, j, j1, j2, jnext, k, n1, n2
191 REAL bignum, eps, rec, scaloc, si, smin, sminw,
192 $ smlnum, sr, tjj, tmp, xj, xmax, xnorm, z
195 REAL d( 2, 2 ), v( 2, 2 )
223 smlnum =
slamch(
'S' ) / eps
224 bignum = one / smlnum
226 xnorm =
slange(
'M', n, n, t, ldt, d )
228 $ xnorm = max( xnorm, abs( w ),
slange(
'M', n, 1, b, n, d ) )
229 smin = max( smlnum, eps*xnorm )
236 work( j ) =
sasum( 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 sscal( 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 sscal( n, rec, x, 1 )
306 x( j1 ) = x( j1 ) / tmp
314 IF( work( j1 ).GT.( bignum-xmax )*rec )
THEN
315 CALL sscal( n, rec, x, 1 )
320 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
334 CALL slaln2( .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 sscal( 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 sscal( n, rec, x, 1 )
363 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
364 CALL saxpy( 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 sscal( n, rec, x, 1 )
408 x( j1 ) = x( j1 ) -
sdot( 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 sscal( 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 sscal( n, rec, x, 1 )
448 d( 1, 1 ) = x( j1 ) -
sdot( j1-1, t( 1, j1 ), 1, x,
450 d( 2, 1 ) = x( j2 ) -
sdot( j1-1, t( 1, j2 ), 1, x,
453 CALL slaln2( .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 sscal( 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 sscal( n2, rec, x, 1 )
521 CALL sladiv( 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 sscal( n2, rec, x, 1 )
538 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
539 CALL saxpy( 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 slaln2( .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 sscal( 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 sscal( n2, rec, x, 1 )
592 CALL saxpy( j1-1, -x( j1 ), t( 1, j1 ), 1, x, 1 )
593 CALL saxpy( j1-1, -x( j2 ), t( 1, j2 ), 1, x, 1 )
595 CALL saxpy( j1-1, -x( n+j1 ), t( 1, j1 ), 1,
597 CALL saxpy( 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 sscal( n2, rec, x, 1 )
650 x( j1 ) = x( j1 ) -
sdot( j1-1, t( 1, j1 ), 1, x, 1 )
651 x( n+j1 ) = x( n+j1 ) -
sdot( 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 sscal( n2, rec, x, 1 )
682 CALL sladiv( 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 sscal( n2, rec, x, 1 )
706 d( 1, 1 ) = x( j1 ) -
sdot( j1-1, t( 1, j1 ), 1, x,
708 d( 2, 1 ) = x( j2 ) -
sdot( j1-1, t( 1, j2 ), 1, x,
710 d( 1, 2 ) = x( n+j1 ) -
sdot( j1-1, t( 1, j1 ), 1,
712 d( 2, 2 ) = x( n+j2 ) -
sdot( 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 slaln2( .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 sscal( 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 )
real function sdot(N, SX, INCX, SY, INCY)
SDOT
integer function isamax(N, SX, INCX)
ISAMAX
subroutine sladiv(A, B, C, D, P, Q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
subroutine slaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function sasum(N, SX, INCX)
SASUM
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sscal(N, SA, SX, INCX)
SSCAL
real function slamch(CMACH)
SLAMCH