Paul BIGNIER [Fri, 26 Jul 2013 09:36:53 +0000 (11:36 +0200)]
Follow-up from https://codereview.scilab.org/#/c/11858/

Bringing corrections to get closer to IDA,
in order to cleanly return the root flags.

These changes are checkable by comparing ida.c and ddaskr.f rootfinding features,
which are basically the same.

Change-Id: Id86d2361278d9f40134668af2dddf77acfa1cecb

index 210a8a2..04c10ae 100644 (file)
@@ -2640,10 +2640,10 @@ C If a root was found on the previous step, evaluate R0 = R(T0). -------
ZROOT = .FALSE.
DO 210 I = 1,NRT
IF (ABS(R0(I)) .EQ. ZERO) THEN
-          ZROOT = .TRUE.
c*********SCILAB ENTERPRISES INPUT
c******** Like with JOB = 1, simply initialize JROOT to 0,
c******** mask the ones that are null at the left endpoint
+c          ZROOT = .TRUE.
ELSE
JROOT(I) = 0
@@ -2687,7 +2687,7 @@ c**** Here, the calculaltion mode can save some computations
IF ((TOUT - TN)*H .GE. ZERO) THEN
T1 = TN
GO TO 330
-         ENDIF
+      ENDIF
T1 = TOUT
ELSE
T1 = TN
@@ -2843,25 +2843,24 @@ c****** If a root function is null at both endpoints, flag it as STUCK.
c        ZROOT = .TRUE.
GO TO 120
C At this point, R0(i) has been checked and cannot be zero. ------------
- 110    IF (SIGN(1.0D0,R0(I)) .EQ. SIGN(1.0D0,R1(I))) GO TO 120
c******** Here, test if some roots get UNSTUCK.
-          IF (JROOT(i).EQ.MASKED) IUNSTUCK = I
-          T2 = ABS(R1(I)/(R1(I)-R0(I)))
-          IF (T2 .LE. TMAX) GO TO 120
-            TMAX = T2
-            IMAX = I
+ 110    IF (JROOT(I).EQ.MASKED) IUNSTUCK = I
+        IF (R0(I)*R1(I) .GT. ZERO) GO TO 120
+        T2 = ABS(R1(I)/(R1(I)-R0(I)))
+        IF (T2 .LE. TMAX) GO TO 120
+        TMAX = T2
+        IMAX = I
120    CONTINUE
IF (IMAX .GT. 0) GO TO 130
c******* STUCK and UNSTUCK root functions count as sign changes.
IF (ISTUCK .GT. 0) THEN
-            SGNCHG = .TRUE.
IMAX = ISTUCK
+            GO TO 130
ELSEIF (IUNSTUCK .GT. 0) THEN
-            SGNCHG = .TRUE.
IMAX = IUNSTUCK
-         ELSE
-      SGNCHG = .FALSE.
+            GO TO 130
ENDIF
+      SGNCHG = .FALSE.
GO TO 140
130  SGNCHG = .TRUE.
140  IF (.NOT. SGNCHG) GO TO 400
@@ -2912,28 +2911,27 @@ C Check to see in which interval R changes sign. -----------------------
ZROOT = .FALSE.
DO 220 I = 1,NRT
IF (ABS(RX(I)) .GT. ZERO) GO TO 210
-        IF (ABS(RX(I)) .EQ. ZERO .AND. JROOT(i) .NE. MASKED) ISTUCK = I
+        IF (ABS(RX(I)).EQ.ZERO .AND. JROOT(I).NE.MASKED) ISTUCK = I
c        ZROOT = .TRUE.
GO TO 220
C Neither R0(i) nor RX(i) can be zero at this point. -------------------
- 210    IF (SIGN(1.0D0,R0(I)) .EQ. SIGN(1.0D0,RX(I))) GO TO 220
-          IF (JROOT(i).EQ.MASKED) IUNSTUCK = I
-          T2 = ABS(RX(I)/(RX(I) - R0(I)))
-          IF (T2 .LE. TMAX) GO TO 220
-            TMAX = T2
-            IMAX = I
+ 210    IF (JROOT(I).EQ.MASKED) IUNSTUCK = I
+        IF (R0(I)*RX(I) .GT. 0) GO TO 220
+        T2 = ABS(RX(I)/(RX(I) - R0(I)))
+        IF (T2 .LE. TMAX) GO TO 220
+          TMAX = T2
+          IMAX = I
220    CONTINUE
IF (IMAX .GT. 0) GO TO 230
IF (ISTUCK .GT. 0) THEN
-            SGNCHG = .TRUE.
IMAX = ISTUCK
+            GO TO 230
ELSEIF (IUNSTUCK .GT. 0) THEN
-            SGNCHG = .TRUE.
IMAX = IUNSTUCK
-         ELSE
+            GO TO 230
+         ENDIF
SGNCHG = .FALSE.
IMAX = IMXOLD
-         ENDIF
GO TO 240
230  SGNCHG = .TRUE.
240  NXLAST = LAST
C Return with X1 as the root.  Set JROOT.  Set X = X1 and RX = R1. -----
300  JFLAG = 2
X = X1
+      CALL DCOPY (NRT, R1, 1, RX, 1)
c**** The following part unmasks root functions if needed
c**** and gives final values to JROOT
UMROOT = .FALSE.
-      CALL DCOPY (NRT, R1, 1, RX, 1)
DO 320 I = 1,NRT
c        JROOT(I) = 0
IF (ABS(R1(I)) .EQ. ZERO) THEN
-            JROOT(I) = -SIGN(1.0D0,R0(I))
+            IF (R0(I).GT.ZERO) THEN
+               JROOT(I) = -1
+            ELSE
+               JROOT(I) = 1
+            ENDIF
ZROOT = .TRUE.
GO TO 320
ENDIF
-          IF (SIGN(1.0D0,R0(I)) .NE. SIGN(1.0D0,R1(I))) THEN
-            JROOT(I) = SIGN(1.0D0,R1(I) - R0(I))
+          IF (R0(I)*R1(I).LT.ZERO) THEN
+            JROOT(I) = SIGN(1.0D0,R1(I))
+            ZROOT = .TRUE.
ELSE
JROOT(I) = 0
ZROOT = .FALSE.
ENDIF
ELSE
IF (ABS(R1(I)) .NE. ZERO) THEN
-            JROOT(I) = SIGN(2.0D0,R1(I))
+            IF (R1(I) .GT. ZERO) THEN
+              JROOT(I) = 2
+            ELSE
+              JROOT(I) = -2
+            ENDIF
UMROOT = .TRUE.
-            ZROOT = .FALSE.
ELSE
JROOT(I) = 0
-            ZROOT = .FALSE.
ENDIF
+          ZROOT = .FALSE.
ENDIF
320    CONTINUE
IF (ZROOT) THEN