--- rpl/lapack/lapack/dlatrs.f	2018/05/29 07:18:01	1.19
+++ rpl/lapack/lapack/dlatrs.f	2023/08/07 08:39:00	1.20
@@ -158,8 +158,6 @@
 *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.
 *
-*> \date December 2016
-*
 *> \ingroup doubleOTHERauxiliary
 *
 *> \par Further Details:
@@ -238,10 +236,9 @@
       SUBROUTINE DLATRS( UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE,
      $                   CNORM, INFO )
 *
-*  -- LAPACK auxiliary routine (version 3.7.0) --
+*  -- LAPACK auxiliary routine --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     December 2016
 *
 *     .. Scalar Arguments ..
       CHARACTER          DIAG, NORMIN, TRANS, UPLO
@@ -267,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
@@ -307,6 +304,7 @@
 *
 *     Quick return if possible
 *
+      SCALE = ONE
       IF( N.EQ.0 )
      $   RETURN
 *
@@ -314,7 +312,6 @@
 *
       SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
       BIGNUM = ONE / SMLNUM
-      SCALE = ONE
 *
       IF( LSAME( NORMIN, 'N' ) ) THEN
 *
@@ -346,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