update BLAS/LAPACK to version 3.11.0 from 22 Nov 2022

lib/linalg/dgelss.f

@@ -254,11 +254,11 @@
 *
 *              Compute space needed for DGEQRF
                CALL DGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
-               LWORK_DGEQRF=DUM(1)
+               LWORK_DGEQRF = INT( DUM(1) )
 *              Compute space needed for DORMQR
                CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, DUM(1), B,
      $                   LDB, DUM(1), -1, INFO )
-               LWORK_DORMQR=DUM(1)
+               LWORK_DORMQR = INT( DUM(1) )
                MM = N
                MAXWRK = MAX( MAXWRK, N + LWORK_DGEQRF )
                MAXWRK = MAX( MAXWRK, N + LWORK_DORMQR )
@@ -273,15 +273,15 @@
 *              Compute space needed for DGEBRD
                CALL DGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1),
      $                      DUM(1), DUM(1), -1, INFO )
-               LWORK_DGEBRD=DUM(1)
+               LWORK_DGEBRD = INT( DUM(1) )
 *              Compute space needed for DORMBR
                CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, DUM(1),
      $                B, LDB, DUM(1), -1, INFO )
-               LWORK_DORMBR=DUM(1)
+               LWORK_DORMBR = INT( DUM(1) )
 *              Compute space needed for DORGBR
                CALL DORGBR( 'P', N, N, N, A, LDA, DUM(1),
      $                   DUM(1), -1, INFO )
-               LWORK_DORGBR=DUM(1)
+               LWORK_DORGBR = INT( DUM(1) )
 *              Compute total workspace needed
                MAXWRK = MAX( MAXWRK, 3*N + LWORK_DGEBRD )
                MAXWRK = MAX( MAXWRK, 3*N + LWORK_DORMBR )
@@ -305,23 +305,23 @@
 *                 Compute space needed for DGELQF
                   CALL DGELQF( M, N, A, LDA, DUM(1), DUM(1),
      $                -1, INFO )
-                  LWORK_DGELQF=DUM(1)
+                  LWORK_DGELQF = INT( DUM(1) )
 *                 Compute space needed for DGEBRD
                   CALL DGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
      $                      DUM(1), DUM(1), -1, INFO )
-                  LWORK_DGEBRD=DUM(1)
+                  LWORK_DGEBRD = INT( DUM(1) )
 *                 Compute space needed for DORMBR
                   CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA,
      $                DUM(1), B, LDB, DUM(1), -1, INFO )
-                  LWORK_DORMBR=DUM(1)
+                  LWORK_DORMBR = INT( DUM(1) )
 *                 Compute space needed for DORGBR
                   CALL DORGBR( 'P', M, M, M, A, LDA, DUM(1),
      $                   DUM(1), -1, INFO )
-                  LWORK_DORGBR=DUM(1)
+                  LWORK_DORGBR = INT( DUM(1) )
 *                 Compute space needed for DORMLQ
                   CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, DUM(1),
      $                 B, LDB, DUM(1), -1, INFO )
-                  LWORK_DORMLQ=DUM(1)
+                  LWORK_DORMLQ = INT( DUM(1) )
 *                 Compute total workspace needed
                   MAXWRK = M + LWORK_DGELQF
                   MAXWRK = MAX( MAXWRK, M*M + 4*M + LWORK_DGEBRD )
@@ -341,15 +341,15 @@
 *                 Compute space needed for DGEBRD
                   CALL DGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
      $                      DUM(1), DUM(1), -1, INFO )
-                  LWORK_DGEBRD=DUM(1)
+                  LWORK_DGEBRD = INT( DUM(1) )
 *                 Compute space needed for DORMBR
                   CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, A, LDA,
      $                DUM(1), B, LDB, DUM(1), -1, INFO )
-                  LWORK_DORMBR=DUM(1)
+                  LWORK_DORMBR = INT( DUM(1) )
 *                 Compute space needed for DORGBR
                   CALL DORGBR( 'P', M, N, M, A, LDA, DUM(1),
      $                   DUM(1), -1, INFO )
-                  LWORK_DORGBR=DUM(1)
+                  LWORK_DORGBR = INT( DUM(1) )
                   MAXWRK = 3*M + LWORK_DGEBRD
                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORMBR )
                   MAXWRK = MAX( MAXWRK, 3*M + LWORK_DORGBR )

lib/linalg/dlaed4.f

@@ -328,9 +328,12 @@
          IF( C.LT.ZERO )
      $      C = ABS( C )
          IF( C.EQ.ZERO ) THEN
-*          ETA = B/A
+*           ETA = B/A
 *           ETA = RHO - TAU
-            ETA = DLTUB - TAU
+*           ETA = DLTUB - TAU
+*
+*           Update proposed by Li, Ren-Cang:
+            ETA = -W / ( DPSI+DPHI )
          ELSE IF( A.GE.ZERO ) THEN
             ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
          ELSE

lib/linalg/dlascl.f

@@ -272,6 +272,8 @@
          ELSE
             MUL = CTOC / CFROMC
             DONE = .TRUE.
+            IF (MUL .EQ. ONE)
+     $         RETURN
          END IF
       END IF
 *

lib/linalg/dlatrs.f

@@ -264,8 +264,8 @@
 *     .. External Functions ..
       LOGICAL            LSAME
       INTEGER            IDAMAX
-      DOUBLE PRECISION   DASUM, DDOT, DLAMCH
+      DOUBLE PRECISION   DASUM, DDOT, DLAMCH, DLANGE
-      EXTERNAL           LSAME, IDAMAX, DASUM, DDOT, DLAMCH
+      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

lib/linalg/dorgbr.f

@@ -232,7 +232,7 @@
                END IF
             END IF
          END IF
-         LWKOPT = WORK( 1 )
+         LWKOPT = INT( WORK( 1 ) )
          LWKOPT = MAX (LWKOPT, MN)
       END IF
 *

lib/linalg/dscal.f

@@ -93,11 +93,14 @@
 *
 *     .. Local Scalars ..
       INTEGER I,M,MP1,NINCX
+*     .. Parameters ..
+      DOUBLE PRECISION ONE
+      PARAMETER (ONE=1.0D+0)
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC MOD
 *     ..
-      IF (N.LE.0 .OR. INCX.LE.0) RETURN
+      IF (N.LE.0 .OR. INCX.LE.0 .OR. DA.EQ.ONE) RETURN
       IF (INCX.EQ.1) THEN
 *
 *        code for increment equal to 1

lib/linalg/dsyevd.f

@@ -257,7 +257,7 @@
                LWMIN = 2*N + 1
             END IF
             LOPT = MAX( LWMIN, 2*N +
-     $                  ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
+     $                  N*ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
             LIOPT = LIWMIN
          END IF
          WORK( 1 ) = LOPT

lib/linalg/dsygvd.f

@@ -330,8 +330,8 @@
       CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
       CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK,
      $             INFO )
-      LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) )
+      LOPT = INT( MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) )
-      LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
+      LIOPT = INT( MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) )
 *
       IF( WANTZ .AND. INFO.EQ.0 ) THEN
 *

lib/linalg/ieeeck.f

@@ -41,7 +41,7 @@
 *> \param[in] ISPEC
 *> \verbatim
 *>          ISPEC is INTEGER
-*>          Specifies whether to test just for inifinity arithmetic
+*>          Specifies whether to test just for infinity arithmetic
 *>          or whether to test for infinity and NaN arithmetic.
 *>          = 0: Verify infinity arithmetic only.
 *>          = 1: Verify infinity and NaN arithmetic.

lib/linalg/ilaenv.f

@@ -469,6 +469,15 @@
             ELSE
                NB = 64
             END IF
+         ELSE IF( C3.EQ.'SYL' ) THEN
+*           The upper bound is to prevent overly aggressive scaling.
+            IF( SNAME ) THEN
+               NB = MIN( MAX( 48, INT( ( MIN( N1, N2 ) * 16 ) / 100) ),
+     $                   240 )
+            ELSE
+               NB = MIN( MAX( 24, INT( ( MIN( N1, N2 ) * 8 ) / 100) ),
+     $                   80 )
+            END IF
          END IF
       ELSE IF( C2.EQ.'LA' ) THEN
          IF( C3.EQ.'UUM' ) THEN
@@ -477,6 +486,12 @@
             ELSE
                NB = 64
             END IF
+         ELSE IF( C3.EQ.'TRS' ) THEN
+            IF( SNAME ) THEN
+               NB = 32
+            ELSE
+               NB = 32
+            END IF
          END IF
       ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
          IF( C3.EQ.'EBZ' ) THEN

lib/linalg/zdscal.f

@@ -92,17 +92,20 @@
 *
 *     .. Local Scalars ..
       INTEGER I,NINCX
+*     .. Parameters ..
+      DOUBLE PRECISION ONE
+      PARAMETER (ONE=1.0D+0)
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC DCMPLX
+      INTRINSIC DBLE, DCMPLX, DIMAG
 *     ..
-      IF (N.LE.0 .OR. INCX.LE.0) RETURN
+      IF (N.LE.0 .OR. INCX.LE.0 .OR. DA.EQ.ONE) RETURN
       IF (INCX.EQ.1) THEN
 *
 *        code for increment equal to 1
 *
          DO I = 1,N
-            ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
+            ZX(I) = DCMPLX(DA*DBLE(ZX(I)),DA*DIMAG(ZX(I)))
          END DO
       ELSE
 *
@@ -110,7 +113,7 @@
 *
          NINCX = N*INCX
          DO I = 1,NINCX,INCX
-            ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
+            ZX(I) = DCMPLX(DA*DBLE(ZX(I)),DA*DIMAG(ZX(I)))
          END DO
       END IF
       RETURN

lib/linalg/zheevd.f

@@ -284,7 +284,7 @@
                LIWMIN = 1
             END IF
             LOPT = MAX( LWMIN, N +
-     $                  ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
+     $                  N*ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
             LROPT = LRWMIN
             LIOPT = LIWMIN
          END IF

lib/linalg/zlascl.f

@@ -272,6 +272,8 @@
          ELSE
             MUL = CTOC / CFROMC
             DONE = .TRUE.
+            IF (MUL .EQ. ONE)
+     $         RETURN
          END IF
       END IF
 *

lib/linalg/zscal.f

@@ -93,7 +93,11 @@
 *     .. Local Scalars ..
       INTEGER I,NINCX
 *     ..
-      IF (N.LE.0 .OR. INCX.LE.0) RETURN
+*     .. Parameters ..
+      COMPLEX*16 ONE
+      PARAMETER (ONE= (1.0D+0,0.0D+0))
+*     ..
+      IF (N.LE.0 .OR. INCX.LE.0 .OR. ZA.EQ.ONE) RETURN
       IF (INCX.EQ.1) THEN
 *
 *        code for increment equal to 1