From 223882bcd6fa70ee1ae93064670e6ffe011d413d Mon Sep 17 00:00:00 2001
From: Phillip Berndt <phillip.berndt@googlemail.com>
Date: Wed, 29 Jul 2015 17:44:52 +0200
Subject: [PATCH] Redirect all error output into the Python exception
---
lib/radau.f | 70 ++++++++++++++++++++++++++++++-----------------------
pyradau13.c | 17 ++++++++++---
test.py | 5 ++++
3 files changed, 58 insertions(+), 34 deletions(-)
diff --git a/lib/radau.f b/lib/radau.f
index f2203e0..17e2745 100644
--- a/lib/radau.f
+++ b/lib/radau.f
@@ -3,7 +3,8 @@
& JAC ,IJAC,MLJAC,MUJAC,
& MAS ,IMAS,MLMAS,MUMAS,
& SOLOUT,IOUT,
- & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
+ & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID,
+ & ERRBUF)
C ----------------------------------------------------------
C NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
C SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS
@@ -394,6 +395,8 @@ C *** *** *** *** *** *** *** *** *** *** *** *** ***
DIMENSION RPAR(*),IPAR(*)
LOGICAL IMPLCT,JBAND,ARRET,STARTN,PRED
EXTERNAL FCN,JAC,MAS,SOLOUT
+C -- pberndt XXX Added ERRBUF parameter to store WRITE output
+ CHARACTER(1024) ERRBUF
C *** *** *** *** *** *** ***
C SETTING THE PARAMETERS
C *** *** *** *** *** *** ***
@@ -428,7 +431,7 @@ C -------- NUMBER MAXIMAL AND MINIMAL OF STAGES NS
ELSE
NSUS=IWORK(13)
IF (NSUS.LE.0.OR.NS.GE.8.OR.NS.EQ.2.OR.NS.EQ.4.OR.NS.EQ.6) THEN
- WRITE(6,*)' WRONG INPUT IWORK(13)=',IWORK(13)
+ WRITE(ERRBUF,*)' WRONG INPUT IWORK(13)=',IWORK(13)
ARRET=.TRUE.
END IF
END IF
@@ -438,7 +441,7 @@ C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
ELSE
NMAX=IWORK(2)
IF (NMAX.LE.0) THEN
- WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
+ WRITE(ERRBUF,*)' WRONG INPUT IWORK(2)=',IWORK(2)
ARRET=.TRUE.
END IF
END IF
@@ -448,7 +451,7 @@ C -------- NIT MAXIMAL NUMBER OF NEWTON ITERATIONS
ELSE
NIT=IWORK(3)
IF (NIT.LE.0.OR.NIT.GT.50) THEN
- WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
+ WRITE(ERRBUF,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
ARRET=.TRUE.
END IF
END IF
@@ -464,7 +467,8 @@ C -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS
NIND3=IWORK(7)
IF (NIND1.EQ.0) NIND1=N
IF (NIND1+NIND2+NIND3.NE.N) THEN
- WRITE(6,*)' CURIOUS INPUT FOR IWORK(5,6,7)=',NIND1,NIND2,NIND3
+ WRITE(ERRBUF,*)' CURIOUS INPUT FOR IWORK(5,6,7)=',
+ & NIND1,NIND2,NIND3
ARRET=.TRUE.
END IF
C -------- PRED STEP SIZE CONTROL
@@ -480,7 +484,7 @@ C -------- PARAMETER FOR SECOND ORDER EQUATIONS
IF (M1.EQ.0) M2=N
IF (M2.EQ.0) M2=M1
IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN
- WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2
+ WRITE(ERRBUF,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2
ARRET=.TRUE.
END IF
C -------- UROUND SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0
@@ -489,20 +493,21 @@ C -------- UROUND SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0
ELSE
UROUND=WORK(1)
IF (UROUND.LE.1.0D-19.OR.UROUND.GE.1.0D0) THEN
- WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
+ WRITE(ERRBUF,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',
+ & WORK(1)
ARRET=.TRUE.
END IF
END IF
C --------- CHECK IF TOLERANCES ARE O.K.
IF (ITOL.EQ.0) THEN
IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN
- WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
+ WRITE(ERRBUF,*) ' TOLERANCES ARE TOO SMALL'
ARRET=.TRUE.
END IF
ELSE
DO I=1,N
IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN
- WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
+ WRITE(ERRBUF,*) ' TOLERANCES(',I,') ARE TOO SMALL'
ARRET=.TRUE.
END IF
END DO
@@ -513,7 +518,7 @@ C --------- SAFE SAFETY FACTOR IN STEP SIZE PREDICTION
ELSE
SAFE=WORK(2)
IF (SAFE.LE.0.001D0.OR.SAFE.GE.1.0D0) THEN
- WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
+ WRITE(ERRBUF,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
ARRET=.TRUE.
END IF
END IF
@@ -523,7 +528,7 @@ C ------ THET DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
ELSE
THET=WORK(3)
IF (THET.GE.1.0D0) THEN
- WRITE(6,*)' CURIOUS INPUT FOR WORK(3)=',WORK(3)
+ WRITE(ERRBUF,*)' CURIOUS INPUT FOR WORK(3)=',WORK(3)
ARRET=.TRUE.
END IF
END IF
@@ -539,7 +544,7 @@ C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
QUOT2=WORK(6)
END IF
IF (QUOT1.GT.1.0D0.OR.QUOT2.LT.1.0D0) THEN
- WRITE(6,*)' CURIOUS INPUT FOR WORK(5,6)=',QUOT1,QUOT2
+ WRITE(ERRBUF,*)' CURIOUS INPUT FOR WORK(5,6)=',QUOT1,QUOT2
ARRET=.TRUE.
END IF
C -------- MAXIMAL STEP SIZE
@@ -560,7 +565,7 @@ C ------- FACL,FACR PARAMETERS FOR STEP SIZE SELECTION
FACR=1.D0/WORK(9)
END IF
IF (FACL.LT.1.0D0.OR.FACR.GT.1.0D0) THEN
- WRITE(6,*)' CURIOUS INPUT WORK(8,9)=',WORK(8),WORK(9)
+ WRITE(ERRBUF,*)' CURIOUS INPUT WORK(8,9)=',WORK(8),WORK(9)
ARRET=.TRUE.
END IF
C -------- PARAMETERS FOR ORDER SELECTION STRATEGY
@@ -616,8 +621,8 @@ C -- MASS MATRIX
END IF
C ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC"
IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
- WRITE (6,*) 'BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF
- & "JAC"'
+ WRITE(ERRBUF,*) 'BANDWITH OF "MAS" NOT SMALLER THAN
+ & BANDWITH OF "JAC"'
ARRET=.TRUE.
END IF
ELSE
@@ -632,8 +637,8 @@ C ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC"
LDMAS2=MAX(1,LDMAS)
C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN
- WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH
- &FULL JACOBIAN'
+ WRITE(ERRBUF,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS
+ &WITH FULL JACOBIAN'
ARRET=.TRUE.
END IF
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
@@ -652,7 +657,8 @@ C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
C ------ TOTAL STORAGE REQUIREMENT -----------
ISTORE=IEE+NMEE*LDE1-1
IF(ISTORE.GT.LWORK)THEN
- WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
+ WRITE(ERRBUF,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',
+ & ISTORE
ARRET=.TRUE.
END IF
C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
@@ -662,7 +668,8 @@ C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
C --------- TOTAL REQUIREMENT ---------------
ISTORE=IEIPH+NM1-1
IF (ISTORE.GT.LIWORK) THEN
- WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
+ WRITE(ERRBUF,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',
+ & ISTORE
ARRET=.TRUE.
END IF
C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
@@ -679,7 +686,7 @@ C -------- CALL TO CORE INTEGRATOR ------------
& WORK(IESCAL),WORK(IEFF),WORK(IEJAC),WORK(IEE1),WORK(IEE),
& WORK(IEMAS),WORK(IECON),IWORK(IEIP1),IWORK(IEIP2),
& IWORK(IEIPH),VITU,VITD,HHOU,HHOD,NFCN,NJAC,NSTEP,NACCPT,
- & NREJCT,NDEC,NSOL,RPAR,IPAR)
+ & NREJCT,NDEC,NSOL,RPAR,IPAR,ERRBUF)
IWORK(13)=NSUS
IWORK(14)=NFCN
IWORK(15)=NJAC
@@ -702,7 +709,7 @@ C
& NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1,NSMIN,NSMAX,NNMS,NM1NS,
& NMEE,IMPLCT,BANDED,LDJAC,LDE1,LDMAS,ZZ,Y0,SCAL,FF,FJAC,
& E1,EE2,FMAS,CONT,IP1,IP2,IPHES,VITU,VITD,HHOU,HHOD,
- & NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR)
+ & NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR,ERRBUF)
C ----------------------------------------------------------
C CORE INTEGRATOR FOR RADAU
C PARAMETERS SAME AS IN RADAU WITH WORKSPACE ADDED
@@ -750,6 +757,8 @@ C --- THIS PARAMETER HAS TO BE CHANGED IF NUMBER OF STAGES IS >=7
LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,STARTN,CALHES
LOGICAL INDEX1,INDEX2,INDEX3,LAST,PRED,CHANGE,UNEXP,UNEXN,VARIAB
EXTERNAL FCN,JAC
+C -- pberndt XXX Added ERRBUF parameter to store WRITE output
+ CHARACTER(1024) ERRBUF
C *** *** *** *** *** *** ***
C INITIALISATIONS
C *** *** *** *** *** *** ***
@@ -1701,29 +1710,30 @@ C --- UNEXPECTED STEP-REJECTION
IF (CALJAC) GOTO 20
GOTO 10
C --- FAIL EXIT
+C --- pberndt XXX Replaced the two lines from previously with a single
+C one for the messages
175 CONTINUE
- WRITE(9,979)X
- WRITE(6,*) ' RHS EVALUATES TO NAN'
+ WRITE(ERRBUF,*) 'EXIT OF RADAU AT X=', X, ', RHS EVALUATES TO NAN'
IDID=-42
RETURN
176 CONTINUE
- WRITE(6,979)X
- WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER
+ WRITE(ERRBUF,*) 'EXIT OF RADAU AT X=', X, ', MATRIX IS REPEATEDLY
+ & SINGULAR, IER=',IER
IDID=-4
RETURN
177 CONTINUE
- WRITE(6,979)X
- WRITE(6,*) ' STEP SIZE T0O SMALL, H=',H
+ WRITE(ERRBUF,*) 'EXIT OF RADAU AT X=', X, ', STEP SIZE T0O SMALL,
+ & H=',H
IDID=-3
RETURN
178 CONTINUE
- WRITE(6,979)X
- WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED'
+ WRITE(ERRBUF,*) 'EXIT OF RADAU AT X=', X, ', MORE THAN NMAX = ',
+ & NMAX,'STEPS ARE NEEDED'
IDID=-2
RETURN
C --- EXIT CAUSED BY SOLOUT
179 CONTINUE
-C WRITE(6,979)X
+C WRITE(ERRBUF,979)X
979 FORMAT(' EXIT OF RADAU AT X=',E18.4)
IDID=2
RETURN
diff --git a/pyradau13.c b/pyradau13.c
index 2224e68..1d926de 100644
--- a/pyradau13.c
+++ b/pyradau13.c
@@ -41,13 +41,14 @@ extern "C" {
int *LIWORK, // (2+(NSMAX-1)/2)*N+20
float *RPAR, // User-supplied RHS arguments
int *IPAR, // See RPAR
- int *IDID // Return value
+ int *IDID, // Return value
// IDID= 1 COMPUTATION SUCCESSFUL,
// IDID= 2 COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT)
// IDID=-1 INPUT IS NOT CONSISTENT,
// IDID=-2 LARGER NMAX IS NEEDED,
// IDID=-3 STEP SIZE BECOMES TOO SMALL,
// IDID=-4 MATRIX IS REPEATEDLY SINGULAR.
+ char *ERRBUF // Buffer to hold error messages, at least 1024 bytes
);
double contra_(
@@ -344,6 +345,9 @@ static PyObject *radau13(PyObject *self, PyObject *args, PyObject *kwargs) {
PyObject *list_retval = PyList_New(0);
+ char errbuf[1024];
+ memset(errbuf, 0, sizeof(errbuf));
+
for(current_level = 0; current_level < time_levels; current_level++) {
t_final = ((double *)PyArray_DATA(times_array))[current_level];
@@ -380,13 +384,14 @@ static PyObject *radau13(PyObject *self, PyObject *args, PyObject *kwargs) {
&liwork, // (2+(NSMAX-1)/2)*N+20
NULL, // User-supplied RHS arguments
(int*)&options, // See RPAR
- &idid // Return value
+ &idid, // Return value
// IDID= 1 COMPUTATION SUCCESSFUL,
// IDID= 2 COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT)
// IDID=-1 INPUT IS NOT CONSISTENT,
// IDID=-2 LARGER NMAX IS NEEDED,
// IDID=-3 STEP SIZE BECOMES TOO SMALL,
// IDID=-4 MATRIX IS REPEATEDLY SINGULAR.
+ errbuf
);
t_0 = t_final;
@@ -407,8 +412,12 @@ static PyObject *radau13(PyObject *self, PyObject *args, PyObject *kwargs) {
if(idid != 1 || y_out == NULL) {
if(PyErr_Occurred() == NULL) {
- char err[255];
- sprintf(err, "radau failed with idid = %d (%s)", idid, idid_error_strings[2 - idid]);
+ int i;
+ for(i=sizeof(errbuf)-1; i>0 && (errbuf[i] == ' ' || errbuf[i] == '\n' || errbuf[i] == 0 || errbuf[i] == '\t'); i--) {
+ errbuf[i] = 0;
+ }
+ char err[1024];
+ snprintf(err, sizeof(err), "radau failed with idid = %d (%s)%s%s", idid, idid_error_strings[2 - idid], errbuf ? "\n" : "", errbuf);
PyErr_SetString(PyExc_RuntimeError, err);
}
if(time_levels > 1) {
diff --git a/test.py b/test.py
index f255280..fae9025 100755
--- a/test.py
+++ b/test.py
@@ -90,6 +90,11 @@ class TestIntegration(unittest.TestCase):
retval = radau13(rhs, [ 1, 2 ], 10, mass_matrix=mass_matrix)[0]
self.assertAlmostEqual(retval, np.exp(2 * 10), 2)
+ def test_failure(self):
+ "This problem *is* unsolvable"
+ with self.assertRaises(RuntimeError):
+ radau13(lambda t, x: x**2, 1, 1)
+
if __name__ == '__main__':
unittest.main()
--
GitLab