Set scale early for robust triangular solvers (Reference-LAPACK PR712)
This commit is contained in:
Martin Kroeker
GitHub
parent
1d5a3aff0d
commit
31d2145988
|
|
@ -278,7 +278,7 @@
|
|||
$ CDOTU, CLADIV
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL CAXPY, CSSCAL, CTBSV, SLABAD, SSCAL, XERBLA
|
||||
EXTERNAL CAXPY, CSSCAL, CTBSV, SSCAL, XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
|
||||
|
|
@ -324,17 +324,14 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
* Determine machine dependent parameters to control overflow.
|
||||
*
|
||||
SMLNUM = SLAMCH( 'Safe minimum' )
|
||||
SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
CALL SLABAD( SMLNUM, BIGNUM )
|
||||
SMLNUM = SMLNUM / SLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
|
|||
|
|
@ -274,7 +274,7 @@
|
|||
$ CDOTU, CLADIV
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL CAXPY, CSSCAL, CTRSV, SLABAD, SSCAL, XERBLA
|
||||
EXTERNAL CAXPY, CSSCAL, CTRSV, SSCAL, XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS, AIMAG, CMPLX, CONJG, MAX, MIN, REAL
|
||||
|
|
@ -318,17 +318,14 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
* Determine machine dependent parameters to control overflow.
|
||||
*
|
||||
SMLNUM = SLAMCH( 'Safe minimum' )
|
||||
SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
CALL SLABAD( SMLNUM, BIGNUM )
|
||||
SMLNUM = SMLNUM / SLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
@ -360,8 +357,74 @@
|
|||
IF( TMAX.LE.BIGNUM*HALF ) THEN
|
||||
TSCAL = ONE
|
||||
ELSE
|
||||
TSCAL = HALF / ( SMLNUM*TMAX )
|
||||
CALL SSCAL( N, TSCAL, CNORM, 1 )
|
||||
*
|
||||
* Avoid NaN generation if entries in CNORM exceed the
|
||||
* overflow threshold
|
||||
*
|
||||
IF ( TMAX.LE.SLAMCH('Overflow') ) THEN
|
||||
* Case 1: All entries in CNORM are valid floating-point numbers
|
||||
TSCAL = HALF / ( SMLNUM*TMAX )
|
||||
CALL SSCAL( N, TSCAL, CNORM, 1 )
|
||||
ELSE
|
||||
* Case 2: At least one column norm of A cannot be
|
||||
* represented as a floating-point number. Find the
|
||||
* maximum offdiagonal absolute value
|
||||
* max( |Re(A(I,J))|, |Im(A(I,J)| ). If this entry is
|
||||
* not +/- Infinity, use this value as TSCAL.
|
||||
TMAX = ZERO
|
||||
IF( UPPER ) THEN
|
||||
*
|
||||
* A is upper triangular.
|
||||
*
|
||||
DO J = 2, N
|
||||
DO I = 1, J - 1
|
||||
TMAX = MAX( TMAX, ABS( REAL( A( I, J ) ) ),
|
||||
$ ABS( AIMAG(A ( I, J ) ) ) )
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* A is lower triangular.
|
||||
*
|
||||
DO J = 1, N - 1
|
||||
DO I = J + 1, N
|
||||
TMAX = MAX( TMAX, ABS( REAL( A( I, J ) ) ),
|
||||
$ ABS( AIMAG(A ( I, J ) ) ) )
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
IF( TMAX.LE.SLAMCH('Overflow') ) THEN
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
DO J = 1, N
|
||||
IF( CNORM( J ).LE.SLAMCH('Overflow') ) THEN
|
||||
CNORM( J ) = CNORM( J )*TSCAL
|
||||
ELSE
|
||||
* Recompute the 1-norm of each column without
|
||||
* introducing Infinity in the summation.
|
||||
TSCAL = TWO * TSCAL
|
||||
CNORM( J ) = ZERO
|
||||
IF( UPPER ) THEN
|
||||
DO I = 1, J - 1
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * CABS2( A( I, J ) )
|
||||
END DO
|
||||
ELSE
|
||||
DO I = J + 1, N
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * CABS2( A( I, J ) )
|
||||
END DO
|
||||
END IF
|
||||
TSCAL = TSCAL * HALF
|
||||
END IF
|
||||
END DO
|
||||
ELSE
|
||||
* At least one entry of A is not a valid floating-point
|
||||
* entry. Rely on TRSV to propagate Inf and NaN.
|
||||
CALL CTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
|
||||
RETURN
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
* Compute a bound on the computed solution vector to see if the
|
||||
|
|
|
|||
|
|
@ -310,6 +310,7 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
|
|
@ -317,7 +318,6 @@
|
|||
*
|
||||
SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
|
|||
|
|
@ -264,8 +264,8 @@
|
|||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
INTEGER IDAMAX
|
||||
DOUBLE PRECISION DASUM, DDOT, DLAMCH
|
||||
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH
|
||||
DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
|
||||
EXTERNAL LSAME, IDAMAX, DASUM, DDOT, DLAMCH, DLANGE
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL DAXPY, DSCAL, DTRSV, XERBLA
|
||||
|
|
@ -304,6 +304,7 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
|
|
@ -311,7 +312,6 @@
|
|||
*
|
||||
SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
@ -343,8 +343,67 @@
|
|||
IF( TMAX.LE.BIGNUM ) THEN
|
||||
TSCAL = ONE
|
||||
ELSE
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
CALL DSCAL( N, TSCAL, CNORM, 1 )
|
||||
*
|
||||
* Avoid NaN generation if entries in CNORM exceed the
|
||||
* overflow threshold
|
||||
*
|
||||
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
|
||||
* Case 1: All entries in CNORM are valid floating-point numbers
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
CALL DSCAL( N, TSCAL, CNORM, 1 )
|
||||
ELSE
|
||||
* Case 2: At least one column norm of A cannot be represented
|
||||
* as floating-point number. Find the offdiagonal entry A( I, J )
|
||||
* with the largest absolute value. If this entry is not +/- Infinity,
|
||||
* use this value as TSCAL.
|
||||
TMAX = ZERO
|
||||
IF( UPPER ) THEN
|
||||
*
|
||||
* A is upper triangular.
|
||||
*
|
||||
DO J = 2, N
|
||||
TMAX = MAX( DLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
|
||||
$ TMAX )
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* A is lower triangular.
|
||||
*
|
||||
DO J = 1, N - 1
|
||||
TMAX = MAX( DLANGE( 'M', N-J, 1, A( J+1, J ), 1,
|
||||
$ SUMJ ), TMAX )
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
DO J = 1, N
|
||||
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
|
||||
CNORM( J ) = CNORM( J )*TSCAL
|
||||
ELSE
|
||||
* Recompute the 1-norm without introducing Infinity
|
||||
* in the summation
|
||||
CNORM( J ) = ZERO
|
||||
IF( UPPER ) THEN
|
||||
DO I = 1, J - 1
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * ABS( A( I, J ) )
|
||||
END DO
|
||||
ELSE
|
||||
DO I = J + 1, N
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * ABS( A( I, J ) )
|
||||
END DO
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
ELSE
|
||||
* At least one entry of A is not a valid floating-point entry.
|
||||
* Rely on TRSV to propagate Inf and NaN.
|
||||
CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
|
||||
RETURN
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
* Compute a bound on the computed solution vector to see if the
|
||||
|
|
|
|||
|
|
@ -310,6 +310,7 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
|
|
@ -317,7 +318,6 @@
|
|||
*
|
||||
SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
|
|||
|
|
@ -264,8 +264,8 @@
|
|||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
INTEGER ISAMAX
|
||||
REAL SASUM, SDOT, SLAMCH
|
||||
EXTERNAL LSAME, ISAMAX, SASUM, SDOT, SLAMCH
|
||||
REAL SASUM, SDOT, SLAMCH, SLANGE
|
||||
EXTERNAL LSAME, ISAMAX, SASUM, SDOT, SLAMCH, SLANGE
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL SAXPY, SSCAL, STRSV, XERBLA
|
||||
|
|
@ -304,6 +304,7 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
|
|
@ -311,7 +312,6 @@
|
|||
*
|
||||
SMLNUM = SLAMCH( 'Safe minimum' ) / SLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
@ -343,8 +343,67 @@
|
|||
IF( TMAX.LE.BIGNUM ) THEN
|
||||
TSCAL = ONE
|
||||
ELSE
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
CALL SSCAL( N, TSCAL, CNORM, 1 )
|
||||
*
|
||||
* Avoid NaN generation if entries in CNORM exceed the
|
||||
* overflow threshold
|
||||
*
|
||||
IF ( TMAX.LE.SLAMCH('Overflow') ) THEN
|
||||
* Case 1: All entries in CNORM are valid floating-point numbers
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
CALL SSCAL( N, TSCAL, CNORM, 1 )
|
||||
ELSE
|
||||
* Case 2: At least one column norm of A cannot be represented
|
||||
* as floating-point number. Find the offdiagonal entry A( I, J )
|
||||
* with the largest absolute value. If this entry is not +/- Infinity,
|
||||
* use this value as TSCAL.
|
||||
TMAX = ZERO
|
||||
IF( UPPER ) THEN
|
||||
*
|
||||
* A is upper triangular.
|
||||
*
|
||||
DO J = 2, N
|
||||
TMAX = MAX( SLANGE( 'M', J-1, 1, A( 1, J ), 1, SUMJ ),
|
||||
$ TMAX )
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* A is lower triangular.
|
||||
*
|
||||
DO J = 1, N - 1
|
||||
TMAX = MAX( SLANGE( 'M', N-J, 1, A( J+1, J ), 1,
|
||||
$ SUMJ ), TMAX )
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
IF( TMAX.LE.SLAMCH('Overflow') ) THEN
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
DO J = 1, N
|
||||
IF( CNORM( J ).LE.SLAMCH('Overflow') ) THEN
|
||||
CNORM( J ) = CNORM( J )*TSCAL
|
||||
ELSE
|
||||
* Recompute the 1-norm without introducing Infinity
|
||||
* in the summation
|
||||
CNORM( J ) = ZERO
|
||||
IF( UPPER ) THEN
|
||||
DO I = 1, J - 1
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * ABS( A( I, J ) )
|
||||
END DO
|
||||
ELSE
|
||||
DO I = J + 1, N
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * ABS( A( I, J ) )
|
||||
END DO
|
||||
END IF
|
||||
END IF
|
||||
END DO
|
||||
ELSE
|
||||
* At least one entry of A is not a valid floating-point entry.
|
||||
* Rely on TRSV to propagate Inf and NaN.
|
||||
CALL STRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
|
||||
RETURN
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
* Compute a bound on the computed solution vector to see if the
|
||||
|
|
|
|||
|
|
@ -278,7 +278,7 @@
|
|||
$ ZDOTU, ZLADIV
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV, DLABAD
|
||||
EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTBSV
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
|
||||
|
|
@ -324,17 +324,14 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
* Determine machine dependent parameters to control overflow.
|
||||
*
|
||||
SMLNUM = DLAMCH( 'Safe minimum' )
|
||||
SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
CALL DLABAD( SMLNUM, BIGNUM )
|
||||
SMLNUM = SMLNUM / DLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
|
|||
|
|
@ -274,7 +274,7 @@
|
|||
$ ZDOTU, ZLADIV
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV, DLABAD
|
||||
EXTERNAL DSCAL, XERBLA, ZAXPY, ZDSCAL, ZTRSV
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
|
||||
|
|
@ -318,17 +318,14 @@
|
|||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
SCALE = ONE
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
* Determine machine dependent parameters to control overflow.
|
||||
*
|
||||
SMLNUM = DLAMCH( 'Safe minimum' )
|
||||
SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
CALL DLABAD( SMLNUM, BIGNUM )
|
||||
SMLNUM = SMLNUM / DLAMCH( 'Precision' )
|
||||
BIGNUM = ONE / SMLNUM
|
||||
SCALE = ONE
|
||||
*
|
||||
IF( LSAME( NORMIN, 'N' ) ) THEN
|
||||
*
|
||||
|
|
@ -360,8 +357,74 @@
|
|||
IF( TMAX.LE.BIGNUM*HALF ) THEN
|
||||
TSCAL = ONE
|
||||
ELSE
|
||||
TSCAL = HALF / ( SMLNUM*TMAX )
|
||||
CALL DSCAL( N, TSCAL, CNORM, 1 )
|
||||
*
|
||||
* Avoid NaN generation if entries in CNORM exceed the
|
||||
* overflow threshold
|
||||
*
|
||||
IF ( TMAX.LE.DLAMCH('Overflow') ) THEN
|
||||
* Case 1: All entries in CNORM are valid floating-point numbers
|
||||
TSCAL = HALF / ( SMLNUM*TMAX )
|
||||
CALL DSCAL( N, TSCAL, CNORM, 1 )
|
||||
ELSE
|
||||
* Case 2: At least one column norm of A cannot be
|
||||
* represented as a floating-point number. Find the
|
||||
* maximum offdiagonal absolute value
|
||||
* max( |Re(A(I,J))|, |Im(A(I,J)| ). If this entry is
|
||||
* not +/- Infinity, use this value as TSCAL.
|
||||
TMAX = ZERO
|
||||
IF( UPPER ) THEN
|
||||
*
|
||||
* A is upper triangular.
|
||||
*
|
||||
DO J = 2, N
|
||||
DO I = 1, J - 1
|
||||
TMAX = MAX( TMAX, ABS( DBLE( A( I, J ) ) ),
|
||||
$ ABS( DIMAG(A ( I, J ) ) ) )
|
||||
END DO
|
||||
END DO
|
||||
ELSE
|
||||
*
|
||||
* A is lower triangular.
|
||||
*
|
||||
DO J = 1, N - 1
|
||||
DO I = J + 1, N
|
||||
TMAX = MAX( TMAX, ABS( DBLE( A( I, J ) ) ),
|
||||
$ ABS( DIMAG(A ( I, J ) ) ) )
|
||||
END DO
|
||||
END DO
|
||||
END IF
|
||||
*
|
||||
IF( TMAX.LE.DLAMCH('Overflow') ) THEN
|
||||
TSCAL = ONE / ( SMLNUM*TMAX )
|
||||
DO J = 1, N
|
||||
IF( CNORM( J ).LE.DLAMCH('Overflow') ) THEN
|
||||
CNORM( J ) = CNORM( J )*TSCAL
|
||||
ELSE
|
||||
* Recompute the 1-norm of each column without
|
||||
* introducing Infinity in the summation.
|
||||
TSCAL = TWO * TSCAL
|
||||
CNORM( J ) = ZERO
|
||||
IF( UPPER ) THEN
|
||||
DO I = 1, J - 1
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * CABS2( A( I, J ) )
|
||||
END DO
|
||||
ELSE
|
||||
DO I = J + 1, N
|
||||
CNORM( J ) = CNORM( J ) +
|
||||
$ TSCAL * CABS2( A( I, J ) )
|
||||
END DO
|
||||
END IF
|
||||
TSCAL = TSCAL * HALF
|
||||
END IF
|
||||
END DO
|
||||
ELSE
|
||||
* At least one entry of A is not a valid floating-point
|
||||
* entry. Rely on TRSV to propagate Inf and NaN.
|
||||
CALL ZTRSV( UPLO, TRANS, DIAG, N, A, LDA, X, 1 )
|
||||
RETURN
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
* Compute a bound on the computed solution vector to see if the
|
||||
|
|
|
|||
Loading…
Reference in New Issue