Skip to content

Commit a620672

Browse files
weslleyspereiraweslleyspereira
committed
Add (s,d)ROUNDUP_LWORK functions to fix the workspace computation
1 parent aa631b4 commit a620672

File tree

9 files changed

+184
-20
lines changed

9 files changed

+184
-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: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
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+
*> If LWORK > 2**24, then it will get rounded as a Float.
25+
*> This routine guarantees it is rounded up instead of down by
26+
*> multiplying LWORK by 1+eps, where eps is the relative machine precision.
27+
*>
28+
*> \return DROUNDUP_LWORK
29+
*> \verbatim
30+
*> DROUNDUP_LWORK >= LWORK.
31+
*> \endverbatim
32+
*
33+
* Arguments:
34+
* ==========
35+
*
36+
*> \param[in] LWORK Workspace size.
37+
*
38+
* Authors:
39+
* ========
40+
*
41+
*> \author Weslley Pereira, University of Colorado Denver, USA
42+
*
43+
*> \ingroup auxOTHERauxiliary
44+
*
45+
*> \par Further Details:
46+
* =====================
47+
*>
48+
*> \verbatim
49+
*> This routine was inspired in the method `magma_zmake_lwork` from MAGMA.
50+
*> \see https://bitbucket.org/icl/magma/src/master/control/magma_zauxiliary.cpp
51+
*> \endverbatim
52+
*
53+
* =====================================================================
54+
DOUBLE PRECISION FUNCTION DROUNDUP_LWORK( LWORK )
55+
*
56+
* -- LAPACK auxiliary routine --
57+
* -- LAPACK is a software package provided by Univ. of Tennessee, --
58+
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
59+
*
60+
* .. Scalar Arguments ..
61+
INTEGER LWORK
62+
* ..
63+
*
64+
* =====================================================================
65+
* ..
66+
* .. External Functions ..
67+
DOUBLE PRECISION DLAMCH
68+
EXTERNAL DLAMCH
69+
* ..
70+
* .. Executable Statements ..
71+
* ..
72+
DROUNDUP_LWORK = LWORK * ( 1 + DLAMCH('EPS') )
73+
RETURN
74+
*
75+
* End of DROUNDUP_LWORK
76+
*
77+
END

‎INSTALL/sroundup_lwork.f

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
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+
*> If LWORK > 2**24, then it will get rounded as a Float.
25+
*> This routine guarantees it is rounded up instead of down by
26+
*> multiplying LWORK by 1+eps, where eps is the relative machine precision.
27+
*>
28+
*> float( 16777217 ) == 16777216
29+
*> float( 16777217 * (1.+eps) ) == 16777218
30+
*>
31+
*> \return SROUNDUP_LWORK
32+
*> \verbatim
33+
*> SROUNDUP_LWORK >= LWORK.
34+
*> \endverbatim
35+
*
36+
* Arguments:
37+
* ==========
38+
*
39+
*> \param[in] LWORK Workspace size.
40+
*
41+
* Authors:
42+
* ========
43+
*
44+
*> \author Weslley Pereira, University of Colorado Denver, USA
45+
*
46+
*> \ingroup auxOTHERauxiliary
47+
*
48+
*> \par Further Details:
49+
* =====================
50+
*>
51+
*> \verbatim
52+
*> This routine was inspired in the method `magma_zmake_lwork` from MAGMA.
53+
*> \see https://bitbucket.org/icl/magma/src/master/control/magma_zauxiliary.cpp
54+
*> \endverbatim
55+
*
56+
* =====================================================================
57+
REAL FUNCTION SROUNDUP_LWORK( LWORK )
58+
*
59+
* -- LAPACK auxiliary routine --
60+
* -- LAPACK is a software package provided by Univ. of Tennessee, --
61+
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
62+
*
63+
* .. Scalar Arguments ..
64+
INTEGER LWORK
65+
* ..
66+
*
67+
* =====================================================================
68+
* ..
69+
* .. External Functions ..
70+
REAL SLAMCH
71+
EXTERNAL SLAMCH
72+
* ..
73+
* .. Executable Statements ..
74+
* ..
75+
SROUNDUP_LWORK = LWORK * ( 1 + SLAMCH('EPS') )
76+
RETURN
77+
*
78+
* End of SROUNDUP_LWORK
79+
*
80+
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)