Merge 51eee01 into 3b26987

Reference-LAPACK · Aug 4, 2021 · fe1e922 · fe1e922
2 parents 3b26987 + 51eee01
commit fe1e922
Show file tree

Hide file tree

Showing 9 changed files with 199 additions and 20 deletions.
diff --git a/INSTALL/Makefile b/INSTALL/Makefile
@@ -48,3 +48,6 @@ cleantest:
 
 slamch.o: slamch.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
 dlamch.o: dlamch.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
+
+sroundup_lwork.o: sroundup_lwork.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
+droundup_lwork.o: droundup_lwork.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
diff --git a/INSTALL/droundup_lwork.f b/INSTALL/droundup_lwork.f
@@ -0,0 +1,86 @@
+*> \brief \b DROUNDUP_LWORK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*      DOUBLE PRECISION FUNCTION DROUNDUP_LWORK( LWORK )
+*
+*     .. Scalar Arguments ..
+*      INTEGER          LWORK
+*     ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DROUNDUP_LWORK deals with a subtle bug with returning LWORK as a Float.
+*> This routine guarantees it is rounded up instead of down by
+*> multiplying LWORK by 1+eps when it is necessary, where eps is the relative machine precision.
+*> E.g.,
+*>
+*>        float( 9007199254740993            ) == 9007199254740992
+*>        float( 9007199254740993 ) * (1.+eps) == 9007199254740994
+*>
+*> \return DROUNDUP_LWORK
+*> \verbatim
+*>         DROUNDUP_LWORK >= LWORK.
+*>         DROUNDUP_LWORK is guaranteed to have zero decimal part.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] LWORK Workspace size.
+*
+*  Authors:
+*  ========
+*
+*> \author Weslley Pereira, University of Colorado Denver, USA
+*
+*> \ingroup auxOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>  This routine was inspired in the method `magma_zmake_lwork` from MAGMA.
+*>  \see https://bitbucket.org/icl/magma/src/master/control/magma_zauxiliary.cpp
+*> \endverbatim
+*
+*  =====================================================================
+      DOUBLE PRECISION FUNCTION DROUNDUP_LWORK( LWORK )
+*
+*  -- LAPACK auxiliary routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      INTEGER          LWORK
+*     ..
+*
+* =====================================================================
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC         EPSILON, DBLE, INT
+*     ..
+*     .. Executable Statements ..
+*     ..
+      DROUNDUP_LWORK = DBLE( LWORK )
+*
+      IF( INT( DROUNDUP_LWORK ) .LT. LWORK ) THEN
+*         Force round up of LWORK
+          DROUNDUP_LWORK = DROUNDUP_LWORK * ( 1.0D+0 + EPSILON(0.0D+0) )
+      ENDIF
+*
+      RETURN
+*
+*     End of DROUNDUP_LWORK
+*
+      END
diff --git a/INSTALL/sroundup_lwork.f b/INSTALL/sroundup_lwork.f
@@ -0,0 +1,86 @@
+*> \brief \b SROUNDUP_LWORK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*  Definition:
+*  ===========
+*
+*      REAL             FUNCTION SROUNDUP_LWORK( LWORK )
+*
+*     .. Scalar Arguments ..
+*      INTEGER          LWORK
+*     ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SROUNDUP_LWORK deals with a subtle bug with returning LWORK as a Float.
+*> This routine guarantees it is rounded up instead of down by
+*> multiplying LWORK by 1+eps when it is necessary, where eps is the relative machine precision.
+*> E.g.,
+*>
+*>        float( 16777217            ) == 16777216
+*>        float( 16777217 ) * (1.+eps) == 16777218
+*>
+*> \return SROUNDUP_LWORK
+*> \verbatim
+*>         SROUNDUP_LWORK >= LWORK.
+*>         SROUNDUP_LWORK is guaranteed to have zero decimal part.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] LWORK Workspace size.
+*
+*  Authors:
+*  ========
+*
+*> \author Weslley Pereira, University of Colorado Denver, USA
+*
+*> \ingroup auxOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>  This routine was inspired in the method `magma_zmake_lwork` from MAGMA.
+*>  \see https://bitbucket.org/icl/magma/src/master/control/magma_zauxiliary.cpp
+*> \endverbatim
+*
+*  =====================================================================
+      REAL             FUNCTION SROUNDUP_LWORK( LWORK )
+*
+*  -- LAPACK auxiliary routine --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      INTEGER          LWORK
+*     ..
+*
+* =====================================================================
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC         EPSILON, REAL, INT
+*     ..
+*     .. Executable Statements ..
+*     ..
+      SROUNDUP_LWORK = REAL( LWORK )
+*
+      IF( INT( SROUNDUP_LWORK ) .LT. LWORK ) THEN
+*         Force round up of LWORK
+          SROUNDUP_LWORK = SROUNDUP_LWORK * ( 1.0E+0 + EPSILON(0.0E+0) )
+      ENDIF
+*
+      RETURN
+*
+*     End of SROUNDUP_LWORK
+*
+      END
diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt
@@ -56,7 +56,7 @@ set(SCLAUX
    slaset.f slasq1.f slasq2.f slasq3.f slasq4.f slasq5.f slasq6.f
    slasr.f  slasrt.f slassq.f90 slasv2.f spttrf.f sstebz.f sstedc.f
    ssteqr.f ssterf.f slaisnan.f sisnan.f
-   slartgp.f slartgs.f
+   slartgp.f slartgs.f ../INSTALL/sroundup_lwork.f
    ${SECOND_SRC})
 
 set(DZLAUX
@@ -75,7 +75,7 @@ set(DZLAUX
    dlaset.f dlasq1.f dlasq2.f dlasq3.f dlasq4.f dlasq5.f dlasq6.f
    dlasr.f  dlasrt.f dlassq.f90 dlasv2.f dpttrf.f dstebz.f dstedc.f
    dsteqr.f dsterf.f dlaisnan.f disnan.f
-   dlartgp.f dlartgs.f
+   dlartgp.f dlartgs.f ../INSTALL/droundup_lwork.f
    ../INSTALL/dlamch.f ${DSECOND_SRC})
 
 set(SLASRC

diff --git a/SRC/Makefile b/SRC/Makefile
@@ -90,7 +90,7 @@ SCLAUX = \
    slaset.o slasq1.o slasq2.o slasq3.o slasq4.o slasq5.o slasq6.o \
    slasr.o  slasrt.o slassq.o slasv2.o spttrf.o sstebz.o sstedc.o \
    ssteqr.o ssterf.o slaisnan.o sisnan.o \
-   slartgp.o slartgs.o \
+   slartgp.o slartgs.o ../INSTALL/sroundup_lwork.o \
    ../INSTALL/second_$(TIMER).o
 
 DZLAUX = \
@@ -109,7 +109,7 @@ DZLAUX = \
    dlaset.o dlasq1.o dlasq2.o dlasq3.o dlasq4.o dlasq5.o dlasq6.o \
    dlasr.o  dlasrt.o dlassq.o dlasv2.o dpttrf.o dstebz.o dstedc.o \
    dsteqr.o dsterf.o dlaisnan.o disnan.o \
-   dlartgp.o dlartgs.o \
+   dlartgp.o dlartgs.o ../INSTALL/droundup_lwork.o \
    ../INSTALL/dlamch.o ../INSTALL/dsecnd_$(TIMER).o
 
 SLASRC = \

diff --git a/SRC/cgesdd.f b/SRC/cgesdd.f
@@ -280,8 +280,9 @@ SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME, SISNAN
-      REAL               SLAMCH, CLANGE
-      EXTERNAL           LSAME, SLAMCH, CLANGE, SISNAN
+      REAL               SLAMCH, CLANGE, SROUNDUP_LWORK
+      EXTERNAL           LSAME, SLAMCH, CLANGE, SISNAN, 
+     $                   SROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, SQRT
@@ -617,7 +618,7 @@ SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
          MAXWRK = MAX( MAXWRK, MINWRK )
       END IF
       IF( INFO.EQ.0 ) THEN
-         WORK( 1 ) = MAXWRK
+         WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
          IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
             INFO = -12
          END IF
@@ -2213,7 +2214,7 @@ SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *
 *     Return optimal workspace in WORK(1)
 *
-      WORK( 1 ) = MAXWRK
+      WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
 *
       RETURN
 *

diff --git a/SRC/dgesdd.f b/SRC/dgesdd.f
@@ -266,8 +266,9 @@ SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME, DISNAN
-      DOUBLE PRECISION   DLAMCH, DLANGE
-      EXTERNAL           DLAMCH, DLANGE, LSAME, DISNAN
+      DOUBLE PRECISION   DLAMCH, DLANGE, DROUNDUP_LWORK
+      EXTERNAL           DLAMCH, DLANGE, LSAME, DISNAN, 
+     $                   DROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, SQRT
@@ -568,7 +569,7 @@ SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
          END IF
 
          MAXWRK = MAX( MAXWRK, MINWRK )
-         WORK( 1 ) = MAXWRK
+         WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
 *
          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
             INFO = -12
@@ -1541,7 +1542,7 @@ SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *
 *     Return optimal workspace in WORK(1)
 *
-      WORK( 1 ) = MAXWRK
+      WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
 *
       RETURN
 *

diff --git a/SRC/sgesdd.f b/SRC/sgesdd.f
@@ -266,8 +266,9 @@ SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME, SISNAN
-      REAL               SLAMCH, SLANGE
-      EXTERNAL           SLAMCH, SLANGE, LSAME, SISNAN
+      REAL               SLAMCH, SLANGE, SROUNDUP_LWORK
+      EXTERNAL           SLAMCH, SLANGE, LSAME, SISNAN, 
+     $                   SROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, SQRT
@@ -568,7 +569,7 @@ SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
          END IF
 
          MAXWRK = MAX( MAXWRK, MINWRK )
-         WORK( 1 ) = MAXWRK
+         WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
 *
          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
             INFO = -12
@@ -1541,7 +1542,7 @@ SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *
 *     Return optimal workspace in WORK(1)
 *
-      WORK( 1 ) = MAXWRK
+      WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
 *
       RETURN
 *

diff --git a/SRC/zgesdd.f b/SRC/zgesdd.f
@@ -280,8 +280,9 @@ SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME, DISNAN
-      DOUBLE PRECISION   DLAMCH, ZLANGE
-      EXTERNAL           LSAME, DLAMCH, ZLANGE, DISNAN
+      DOUBLE PRECISION   DLAMCH, ZLANGE, DROUNDUP_LWORK
+      EXTERNAL           LSAME, DLAMCH, ZLANGE, DISNAN, 
+     $                   DROUNDUP_LWORK
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          INT, MAX, MIN, SQRT
@@ -617,7 +618,7 @@ SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
          MAXWRK = MAX( MAXWRK, MINWRK )
       END IF
       IF( INFO.EQ.0 ) THEN
-         WORK( 1 ) = MAXWRK
+         WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
          IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
             INFO = -12
          END IF
@@ -2213,7 +2214,7 @@ SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
 *
 *     Return optimal workspace in WORK(1)
 *
-      WORK( 1 ) = MAXWRK
+      WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
 *
       RETURN
 *