Skip to content

Commit 7a4c63b

Browse files
langoulangou
authored
Merge pull request #605 from weslleyspereira/try-workspace-query
Add (s,d)ROUNDUP_LWORK functions to fix the workspace computation
2 parents 53e8c3d + c7400e5 commit 7a4c63b

File tree

9 files changed

+199
-20
lines changed

9 files changed

+199
-20
lines changed

‎INSTALL/Makefile

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,3 +48,6 @@ cleantest:
4848

4949
slamch.o: slamch.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
5050
dlamch.o: dlamch.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
51+
52+
sroundup_lwork.o: sroundup_lwork.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<
53+
droundup_lwork.o: droundup_lwork.f ; $(FC) $(FFLAGS_NOOPT) -c -o $@ $<

‎INSTALL/droundup_lwork.f

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
*> \brief \b DROUNDUP_LWORK
2+
*
3+
* =========== DOCUMENTATION ===========
4+
*
5+
* Online html documentation available at
6+
* http://www.netlib.org/lapack/explore-html/
7+
*
8+
* Definition:
9+
* ===========
10+
*
11+
* DOUBLE PRECISION FUNCTION DROUNDUP_LWORK( LWORK )
12+
*
13+
* .. Scalar Arguments ..
14+
* INTEGER LWORK
15+
* ..
16+
*
17+
*
18+
*> \par Purpose:
19+
* =============
20+
*>
21+
*> \verbatim
22+
*>
23+
*> DROUNDUP_LWORK deals with a subtle bug with returning LWORK as a Float.
24+
*> This routine guarantees it is rounded up instead of down by
25+
*> multiplying LWORK by 1+eps when it is necessary, where eps is the relative machine precision.
26+
*> E.g.,
27+
*>
28+
*> float( 9007199254740993 ) == 9007199254740992
29+
*> float( 9007199254740993 ) * (1.+eps) == 9007199254740994
30+
*>
31+
*> \return DROUNDUP_LWORK
32+
*> \verbatim
33+
*> DROUNDUP_LWORK >= LWORK.
34+
*> DROUNDUP_LWORK is guaranteed to have zero decimal part.
35+
*> \endverbatim
36+
*
37+
* Arguments:
38+
* ==========
39+
*
40+
*> \param[in] LWORK Workspace size.
41+
*
42+
* Authors:
43+
* ========
44+
*
45+
*> \author Weslley Pereira, University of Colorado Denver, USA
46+
*
47+
*> \ingroup auxOTHERauxiliary
48+
*
49+
*> \par Further Details:
50+
* =====================
51+
*>
52+
*> \verbatim
53+
*> This routine was inspired in the method `magma_zmake_lwork` from MAGMA.
54+
*> \see https://bitbucket.org/icl/magma/src/master/control/magma_zauxiliary.cpp
55+
*> \endverbatim
56+
*
57+
* =====================================================================
58+
DOUBLE PRECISION FUNCTION DROUNDUP_LWORK( LWORK )
59+
*
60+
* -- LAPACK auxiliary routine --
61+
* -- LAPACK is a software package provided by Univ. of Tennessee, --
62+
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
63+
*
64+
* .. Scalar Arguments ..
65+
INTEGER LWORK
66+
* ..
67+
*
68+
* =====================================================================
69+
* ..
70+
* .. Intrinsic Functions ..
71+
INTRINSIC EPSILON, DBLE, INT
72+
* ..
73+
* .. Executable Statements ..
74+
* ..
75+
DROUNDUP_LWORK = DBLE( LWORK )
76+
*
77+
IF( INT( DROUNDUP_LWORK ) .LT. LWORK ) THEN
78+
* Force round up of LWORK
79+
DROUNDUP_LWORK = DROUNDUP_LWORK * ( 1.0D+0 + EPSILON(0.0D+0) )
80+
ENDIF
81+
*
82+
RETURN
83+
*
84+
* End of DROUNDUP_LWORK
85+
*
86+
END

‎INSTALL/sroundup_lwork.f

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
*> \brief \b SROUNDUP_LWORK
2+
*
3+
* =========== DOCUMENTATION ===========
4+
*
5+
* Online html documentation available at
6+
* http://www.netlib.org/lapack/explore-html/
7+
*
8+
* Definition:
9+
* ===========
10+
*
11+
* REAL FUNCTION SROUNDUP_LWORK( LWORK )
12+
*
13+
* .. Scalar Arguments ..
14+
* INTEGER LWORK
15+
* ..
16+
*
17+
*
18+
*> \par Purpose:
19+
* =============
20+
*>
21+
*> \verbatim
22+
*>
23+
*> SROUNDUP_LWORK deals with a subtle bug with returning LWORK as a Float.
24+
*> This routine guarantees it is rounded up instead of down by
25+
*> multiplying LWORK by 1+eps when it is necessary, where eps is the relative machine precision.
26+
*> E.g.,
27+
*>
28+
*> float( 16777217 ) == 16777216
29+
*> float( 16777217 ) * (1.+eps) == 16777218
30+
*>
31+
*> \return SROUNDUP_LWORK
32+
*> \verbatim
33+
*> SROUNDUP_LWORK >= LWORK.
34+
*> SROUNDUP_LWORK is guaranteed to have zero decimal part.
35+
*> \endverbatim
36+
*
37+
* Arguments:
38+
* ==========
39+
*
40+
*> \param[in] LWORK Workspace size.
41+
*
42+
* Authors:
43+
* ========
44+
*
45+
*> \author Weslley Pereira, University of Colorado Denver, USA
46+
*
47+
*> \ingroup auxOTHERauxiliary
48+
*
49+
*> \par Further Details:
50+
* =====================
51+
*>
52+
*> \verbatim
53+
*> This routine was inspired in the method `magma_zmake_lwork` from MAGMA.
54+
*> \see https://bitbucket.org/icl/magma/src/master/control/magma_zauxiliary.cpp
55+
*> \endverbatim
56+
*
57+
* =====================================================================
58+
REAL FUNCTION SROUNDUP_LWORK( LWORK )
59+
*
60+
* -- LAPACK auxiliary routine --
61+
* -- LAPACK is a software package provided by Univ. of Tennessee, --
62+
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
63+
*
64+
* .. Scalar Arguments ..
65+
INTEGER LWORK
66+
* ..
67+
*
68+
* =====================================================================
69+
* ..
70+
* .. Intrinsic Functions ..
71+
INTRINSIC EPSILON, REAL, INT
72+
* ..
73+
* .. Executable Statements ..
74+
* ..
75+
SROUNDUP_LWORK = REAL( LWORK )
76+
*
77+
IF( INT( SROUNDUP_LWORK ) .LT. LWORK ) THEN
78+
* Force round up of LWORK
79+
SROUNDUP_LWORK = SROUNDUP_LWORK * ( 1.0E+0 + EPSILON(0.0E+0) )
80+
ENDIF
81+
*
82+
RETURN
83+
*
84+
* End of SROUNDUP_LWORK
85+
*
86+
END

‎SRC/CMakeLists.txt

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ set(SCLAUX
5656
slaset.f slasq1.f slasq2.f slasq3.f slasq4.f slasq5.f slasq6.f
5757
slasr.f slasrt.f slassq.f90 slasv2.f spttrf.f sstebz.f sstedc.f
5858
ssteqr.f ssterf.f slaisnan.f sisnan.f
59-
slartgp.f slartgs.f
59+
slartgp.f slartgs.f ../INSTALL/sroundup_lwork.f
6060
${SECOND_SRC})
6161

6262
set(DZLAUX
@@ -75,7 +75,7 @@ set(DZLAUX
7575
dlaset.f dlasq1.f dlasq2.f dlasq3.f dlasq4.f dlasq5.f dlasq6.f
7676
dlasr.f dlasrt.f dlassq.f90 dlasv2.f dpttrf.f dstebz.f dstedc.f
7777
dsteqr.f dsterf.f dlaisnan.f disnan.f
78-
dlartgp.f dlartgs.f
78+
dlartgp.f dlartgs.f ../INSTALL/droundup_lwork.f
7979
../INSTALL/dlamch.f ${DSECOND_SRC})
8080

8181
set(SLASRC

‎SRC/Makefile

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -90,7 +90,7 @@ SCLAUX = \
9090
slaset.o slasq1.o slasq2.o slasq3.o slasq4.o slasq5.o slasq6.o \
9191
slasr.o slasrt.o slassq.o slasv2.o spttrf.o sstebz.o sstedc.o \
9292
ssteqr.o ssterf.o slaisnan.o sisnan.o \
93-
slartgp.o slartgs.o \
93+
slartgp.o slartgs.o ../INSTALL/sroundup_lwork.o \
9494
../INSTALL/second_$(TIMER).o
9595

9696
DZLAUX = \
@@ -109,7 +109,7 @@ DZLAUX = \
109109
dlaset.o dlasq1.o dlasq2.o dlasq3.o dlasq4.o dlasq5.o dlasq6.o \
110110
dlasr.o dlasrt.o dlassq.o dlasv2.o dpttrf.o dstebz.o dstedc.o \
111111
dsteqr.o dsterf.o dlaisnan.o disnan.o \
112-
dlartgp.o dlartgs.o \
112+
dlartgp.o dlartgs.o ../INSTALL/droundup_lwork.o \
113113
../INSTALL/dlamch.o ../INSTALL/dsecnd_$(TIMER).o
114114

115115
SLASRC = \

‎SRC/cgesdd.f

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -280,8 +280,9 @@ SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
280280
* ..
281281
* .. External Functions ..
282282
LOGICAL LSAME, SISNAN
283-
REAL SLAMCH, CLANGE
284-
EXTERNAL LSAME, SLAMCH, CLANGE, SISNAN
283+
REAL SLAMCH, CLANGE, SROUNDUP_LWORK
284+
EXTERNAL LSAME, SLAMCH, CLANGE, SISNAN,
285+
$ SROUNDUP_LWORK
285286
* ..
286287
* .. Intrinsic Functions ..
287288
INTRINSIC INT, MAX, MIN, SQRT
@@ -617,7 +618,7 @@ SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
617618
MAXWRK = MAX( MAXWRK, MINWRK )
618619
END IF
619620
IF( INFO.EQ.0 ) THEN
620-
WORK( 1 ) = MAXWRK
621+
WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
621622
IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
622623
INFO = -12
623624
END IF
@@ -2213,7 +2214,7 @@ SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
22132214
*
22142215
* Return optimal workspace in WORK(1)
22152216
*
2216-
WORK( 1 ) = MAXWRK
2217+
WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
22172218
*
22182219
RETURN
22192220
*

‎SRC/dgesdd.f

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -266,8 +266,9 @@ SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
266266
* ..
267267
* .. External Functions ..
268268
LOGICAL LSAME, DISNAN
269-
DOUBLE PRECISION DLAMCH, DLANGE
270-
EXTERNAL DLAMCH, DLANGE, LSAME, DISNAN
269+
DOUBLE PRECISION DLAMCH, DLANGE, DROUNDUP_LWORK
270+
EXTERNAL DLAMCH, DLANGE, LSAME, DISNAN,
271+
$ DROUNDUP_LWORK
271272
* ..
272273
* .. Intrinsic Functions ..
273274
INTRINSIC INT, MAX, MIN, SQRT
@@ -568,7 +569,7 @@ SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
568569
END IF
569570

570571
MAXWRK = MAX( MAXWRK, MINWRK )
571-
WORK( 1 ) = MAXWRK
572+
WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
572573
*
573574
IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
574575
INFO = -12
@@ -1541,7 +1542,7 @@ SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
15411542
*
15421543
* Return optimal workspace in WORK(1)
15431544
*
1544-
WORK( 1 ) = MAXWRK
1545+
WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
15451546
*
15461547
RETURN
15471548
*

‎SRC/sgesdd.f

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -266,8 +266,9 @@ SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
266266
* ..
267267
* .. External Functions ..
268268
LOGICAL LSAME, SISNAN
269-
REAL SLAMCH, SLANGE
270-
EXTERNAL SLAMCH, SLANGE, LSAME, SISNAN
269+
REAL SLAMCH, SLANGE, SROUNDUP_LWORK
270+
EXTERNAL SLAMCH, SLANGE, LSAME, SISNAN,
271+
$ SROUNDUP_LWORK
271272
* ..
272273
* .. Intrinsic Functions ..
273274
INTRINSIC INT, MAX, MIN, SQRT
@@ -568,7 +569,7 @@ SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
568569
END IF
569570

570571
MAXWRK = MAX( MAXWRK, MINWRK )
571-
WORK( 1 ) = MAXWRK
572+
WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
572573
*
573574
IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
574575
INFO = -12
@@ -1541,7 +1542,7 @@ SUBROUTINE SGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
15411542
*
15421543
* Return optimal workspace in WORK(1)
15431544
*
1544-
WORK( 1 ) = MAXWRK
1545+
WORK( 1 ) = SROUNDUP_LWORK( MAXWRK )
15451546
*
15461547
RETURN
15471548
*

‎SRC/zgesdd.f

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -280,8 +280,9 @@ SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
280280
* ..
281281
* .. External Functions ..
282282
LOGICAL LSAME, DISNAN
283-
DOUBLE PRECISION DLAMCH, ZLANGE
284-
EXTERNAL LSAME, DLAMCH, ZLANGE, DISNAN
283+
DOUBLE PRECISION DLAMCH, ZLANGE, DROUNDUP_LWORK
284+
EXTERNAL LSAME, DLAMCH, ZLANGE, DISNAN,
285+
$ DROUNDUP_LWORK
285286
* ..
286287
* .. Intrinsic Functions ..
287288
INTRINSIC INT, MAX, MIN, SQRT
@@ -617,7 +618,7 @@ SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
617618
MAXWRK = MAX( MAXWRK, MINWRK )
618619
END IF
619620
IF( INFO.EQ.0 ) THEN
620-
WORK( 1 ) = MAXWRK
621+
WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
621622
IF( LWORK.LT.MINWRK .AND. .NOT. LQUERY ) THEN
622623
INFO = -12
623624
END IF
@@ -2213,7 +2214,7 @@ SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT,
22132214
*
22142215
* Return optimal workspace in WORK(1)
22152216
*
2216-
WORK( 1 ) = MAXWRK
2217+
WORK( 1 ) = DROUNDUP_LWORK( MAXWRK )
22172218
*
22182219
RETURN
22192220
*

0 commit comments

Comments
 (0)