|
21 | 21 | *> \verbatim
|
22 | 22 | *>
|
23 | 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 | 24 | *> This routine guarantees it is rounded up instead of down by
|
26 |
| -*> multiplying LWORK by 1+eps, where eps is the relative machine precision. |
| 25 | +*> multiplying LWORK by 1+eps when it is necessary, where eps is the relative machine precision. |
| 26 | +*> E.g., |
27 | 27 | *>
|
28 | 28 | *> float( 16777217 ) == 16777216
|
29 |
| -*> float( 16777217 * (1.+eps) ) == 16777218 |
| 29 | +*> float( 16777217 ) * (1.+eps) == 16777218 |
30 | 30 | *>
|
31 | 31 | *> \return SROUNDUP_LWORK
|
32 | 32 | *> \verbatim
|
33 | 33 | *> SROUNDUP_LWORK >= LWORK.
|
| 34 | +*> SROUNDUP_LWORK is guaranteed to have zero decimal part. |
34 | 35 | *> \endverbatim
|
35 | 36 | *
|
36 | 37 | * Arguments:
|
@@ -67,14 +68,14 @@ REAL FUNCTION SROUNDUP_LWORK( LWORK )
|
67 | 68 | * =====================================================================
|
68 | 69 | * ..
|
69 | 70 | * .. Intrinsic Functions ..
|
70 |
| - INTRINSIC DIGITS, RADIX, EPSILON |
| 71 | + INTRINSIC EPSILON, REAL, INT |
71 | 72 | * ..
|
72 | 73 | * .. Executable Statements ..
|
73 | 74 | * ..
|
74 |
| - SROUNDUP_LWORK = LWORK |
| 75 | + SROUNDUP_LWORK = REAL( LWORK ) |
75 | 76 | *
|
76 |
| - IF( SROUNDUP_LWORK .GE. REAL(RADIX(0.0E+0))**DIGITS(0.0E+0) ) THEN |
77 |
| -* If LWORK can't be represented exactly in single precision |
| 77 | + IF( INT( SROUNDUP_LWORK ) .LT. LWORK ) THEN |
| 78 | +* Force round up of LWORK |
78 | 79 | SROUNDUP_LWORK = SROUNDUP_LWORK * ( 1.0E+0 + EPSILON(0.0E+0) )
|
79 | 80 | ENDIF
|
80 | 81 | *
|
|
0 commit comments