165 SUBROUTINE dlaqtr( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
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 )
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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dscal(N, DA, DX, INCX)
DSCAL
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 ...
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dlaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
integer function idamax(N, DX, INCX)
IDAMAX
double precision function dasum(N, DX, INCX)
DASUM