[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Octave-bug-tracker] [bug #52569] unexptected error occured when using d
From: |
Dan Sebald |
Subject: |
[Octave-bug-tracker] [bug #52569] unexptected error occured when using dynare with octave |
Date: |
Fri, 1 Dec 2017 13:43:11 -0500 (EST) |
User-agent: |
Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:55.0) Gecko/20100101 Firefox/55.0 |
Follow-up Comment #2, bug #52569 (project octave):
Thanks for the links.
I've looked into the code a bit, just to get a handle on what the issue is.
SLATEC D9LGIT is a very short routine, so not difficult to comprehend. SLATEC
is a collection of lesser libraries, of which d9lgit belongs to a library
called The Fullerton Function Library FN, which has been updated and ported to
C/C++ in some way but I didn't look at much of that:
https://people.sc.fsu.edu/~jburkardt/f77_src/fn/fn.html
https://people.sc.fsu.edu/~jburkardt/cpp_src/fn/fn.html
http://www.netlib.org/fn/d9lgit.f
The last reference is the actual, more up-to-date code. However, it isn't
much different, aside from using the double versions of some of the routines.
I'm attaching a diff file that makes Octave's d9lgit.f in line with the one in
the link above. (Don't apply this to the repository, though. I don't think
we want to start freely modifying third-party-sourced code.) I don't see any
difference in convergence behavior when changing the ABS to DABS.
The issue is that poisscdf (A,A+E) doesn't converge well when A grows large
and E shrinks, e.g., poisscdf(1000,1001) versus poisscdf(1000,1002). I don't
have time to think of the mathematical ramifications of this property (e.g.,
does this algorithm outright lose accuracy for the mathematical function as A
becomes large?).
That being the case, it appears that simply bumping up the allowable number of
iterations at least avoids the FORTRAN error. That is, I made this change:
diff --git a/liboctave/external/slatec-fn/d9lgit.f
b/liboctave/external/slatec-f
--- a/liboctave/external/slatec-fn/d9lgit.f
+++ b/liboctave/external/slatec-fn/d9lgit.f
@@ -46,7 +46,7 @@ C
R = 0.D0
P = 1.D0
S = P
- DO 20 K=1,200
+ DO 20 K=1,500
FK = K
T = (A+FK)*X*(1.D0+R)
R = T/((AX+FK)*(A1X+FK)-T)
@@ -55,7 +55,7 @@ C
IF (ABS(P).LT.EPS*S) GO TO 30
20 CONTINUE
CALL XERMSG ('SLATEC', 'D9LGIT',
- + 'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 3, 2)
+ + 'NO CONVERGENCE IN 500 TERMS OF CONTINUED FRACTION', 3, 2)
C
30 HSTAR = 1.0D0 - X*S/A1X
IF (HSTAR .LT. SQEPS) CALL XERMSG ('SLATEC', 'D9LGIT',
(END)
and can then get results for larger A. For example:
octave:5> poisscdf (4000,4000)
ans = 0.50421
octave:6> poisscdf (5000,5000)
***MESSAGE FROM ROUTINE D9LGIT IN LIBRARY SLATEC.
***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
* NO CONVERGENCE IN 500 TERMS OF CONTINUED FRACTION
* ERROR NUMBER = 3
*
***END OF MESSAGE
***JOB ABORT DUE TO FATAL ERROR.
0 ERROR MESSAGE SUMMARY
LIBRARY SUBROUTINE MESSAGE START NERR LEVEL COUNT
SLATEC D9LGIT NO CONVERGENCE IN 50 3 2 2
error: gammainc: exception encountered in Fortran subroutine xgammainc_
error: called from
poisscdf at line 60 column 12
octave:6>
So, obviously this algorithm is going to fall apart at some point. Perhaps
the original submitter could give us some insight as to reasonable
expectations with regard to applications.
This "error: Due to a bug in Octave" isn't quite a bug, it would seem. It's
more like a limitation on numerical methods for approximating the function
which perhaps we could stretch out but never totally eliminate.
(file #42546)
_______________________________________________________
Additional Item Attachment:
File name: octave-d9lgit_convergence-djs2017dec01.diff Size:1 KB
_______________________________________________________
Reply to this item at:
<http://savannah.gnu.org/bugs/?52569>
_______________________________________________
Message sent via/by Savannah
http://savannah.gnu.org/