1
SUBROUTINE STEPNF(N,NFE,IFLAG,START,CRASH,HOLD,H,RELERR,
2
& ABSERR,S,Y,YP,YOLD,YPOLD,A,QR,ALPHA,TZ,PIVOT,W,WP,
3
& Z0,Z1,SSPAR,PAR,IPAR)
5
C STEPNF TAKES ONE STEP ALONG THE ZERO CURVE OF THE HOMOTOPY MAP
6
C USING A PREDICTOR-CORRECTOR ALGORITHM. THE PREDICTOR USES A HERMITE
7
C CUBIC INTERPOLANT, AND THE CORRECTOR RETURNS TO THE ZERO CURVE ALONG
8
C THE FLOW NORMAL TO THE DAVIDENKO FLOW. STEPNF ALSO ESTIMATES A
9
C STEP SIZE H FOR THE NEXT STEP ALONG THE ZERO CURVE. NORMALLY
10
C STEPNF IS USED INDIRECTLY THROUGH FIXPNF , AND SHOULD BE CALLED
11
C DIRECTLY ONLY IF IT IS NECESSARY TO MODIFY THE STEPPING ALGORITHM'S
16
C N = DIMENSION OF X AND THE HOMOTOPY MAP.
18
C NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS.
20
C IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE.
22
C START = .TRUE. ON FIRST CALL TO STEPNF , .FALSE. OTHERWISE.
24
C HOLD = ||Y - YOLD||; SHOULD NOT BE MODIFIED BY THE USER.
26
C H = UPPER LIMIT ON LENGTH OF STEP THAT WILL BE ATTEMPTED. H MUST BE
27
C SET TO A POSITIVE NUMBER ON THE FIRST CALL TO STEPNF .
28
C THEREAFTER STEPNF CALCULATES AN OPTIMAL VALUE FOR H , AND H
29
C SHOULD NOT BE MODIFIED BY THE USER.
31
C RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES. THE ITERATION IS
32
C CONSIDERED TO HAVE CONVERGED WHEN A POINT W=(LAMBDA,X) IS FOUND
35
C ||Z|| <= RELERR*||W|| + ABSERR , WHERE
37
C Z IS THE NEWTON STEP TO W=(LAMBDA,X).
39
C S = (APPROXIMATE) ARC LENGTH ALONG THE HOMOTOPY ZERO CURVE UP TO
40
C Y(S) = (LAMBDA(S), X(S)).
42
C Y(1:N+1) = PREVIOUS POINT (LAMBDA(S), X(S)) FOUND ON THE ZERO CURVE OF
45
C YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP
48
C YOLD(1:N+1) = A POINT BEFORE Y ON THE ZERO CURVE OF THE HOMOTOPY MAP.
50
C YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY
53
C A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP.
55
C QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1),
56
C WP(1:N+1) ARE WORK ARRAYS USED FOR THE QR FACTORIZATION (IN THE
57
C NEWTON STEP CALCULATION) AND THE INTERPOLATION.
59
C Z0(1:N+1), Z1(1:N+1) ARE WORK ARRAYS USED FOR THE ESTIMATION OF THE
62
C SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) IS
63
C A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION.
65
C PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
66
C WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
71
C N , A , SSPAR ARE UNCHANGED.
73
C NFE HAS BEEN UPDATED.
76
C = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN.
78
C = 4 IF A JACOBIAN MATRIX WITH RANK < N HAS OCCURRED. THE
79
C ITERATION WAS NOT COMPLETED.
81
C = 6 IF THE ITERATION FAILED TO CONVERGE. W CONTAINS THE LAST
84
C START = .FALSE. ON A NORMAL RETURN.
87
C = .FALSE. ON A NORMAL RETURN.
89
C = .TRUE. IF THE STEP SIZE H WAS TOO SMALL. H HAS BEEN
90
C INCREASED TO AN ACCEPTABLE VALUE, WITH WHICH STEPNF MAY BE
93
C = .TRUE. IF RELERR AND/OR ABSERR WERE TOO SMALL. THEY HAVE
94
C BEEN INCREASED TO ACCEPTABLE VALUES, WITH WHICH STEPNF MAY
97
C HOLD = ||Y - YOLD||.
99
C H = OPTIMAL VALUE FOR NEXT STEP TO BE ATTEMPTED. NORMALLY H SHOULD
100
C NOT BE MODIFIED BY THE USER.
102
C RELERR, ABSERR ARE UNCHANGED ON A NORMAL RETURN.
104
C S = (APPROXIMATE) ARC LENGTH ALONG THE ZERO CURVE OF THE HOMOTOPY MAP
105
C UP TO THE LATEST POINT FOUND, WHICH IS RETURNED IN Y .
107
C Y, YP, YOLD, YPOLD CONTAIN THE TWO MOST RECENT POINTS AND TANGENT
108
C VECTORS FOUND ON THE ZERO CURVE OF THE HOMOTOPY MAP.
111
C CALLS D1MACH , DNRM2 , TANGNF .
113
DOUBLE PRECISION ABSERR,D1MACH,DCALC,DD001,DD0011,DD01,
114
& DD011,DELS,DNRM2,F0,F1,FOURU,FP0,FP1,H,HFAIL,HOLD,HT,
115
& LCALC,QOFS,RCALC,RELERR,RHOLEN,S,TEMP,TWOU
116
INTEGER IFLAG,ITNUM,J,JUDY,LITFH,N,NFE,NP1
117
LOGICAL CRASH,FAIL,START
119
C ***** ARRAY DECLARATIONS. *****
121
DOUBLE PRECISION Y(N+1),YP(N+1),YOLD(N+1),YPOLD(N+1),A(N),
122
& QR(N,N+2),ALPHA(N),TZ(N+1),W(N+1),WP(N+1),Z0(N+1),
123
& Z1(N+1),SSPAR(8),PAR(1)
124
INTEGER PIVOT(N+1),IPAR(1)
126
C ***** END OF DIMENSIONAL INFORMATION. *****
128
C THE LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED BEFORE REDUCING
129
C THE STEP SIZE H MAY BE CHANGED BY CHANGING THE FOLLOWING PARAMETER
133
C DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES.
135
DD01(F0,F1,DELS)=(F1-F0)/DELS
136
DD001(F0,FP0,F1,DELS)=(DD01(F0,F1,DELS)-FP0)/DELS
137
DD011(F0,F1,FP1,DELS)=(FP1-DD01(F0,F1,DELS))/DELS
138
DD0011(F0,FP0,F1,FP1,DELS)=(DD011(F0,F1,FP1,DELS) -
139
& DD001(F0,FP0,F1,DELS))/DELS
140
QOFS(F0,FP0,F1,FP1,DELS,S)=((DD0011(F0,FP0,F1,FP1,DELS)*(S-DELS) +
141
& DD001(F0,FP0,F1,DELS))*S + FP0)*S + F0
148
C THE ARCLENGTH S MUST BE NONNEGATIVE.
149
IF (S .LT. 0.0) RETURN
150
C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE.
151
IF (H .LT. FOURU*(1.0+S)) THEN
155
C IF ERROR TOLERANCES ARE TOO SMALL, INCREASE THEM TO ACCEPTABLE VALUES.
157
IF (.5*(RELERR*TEMP+ABSERR) .GE. TWOU*TEMP) GO TO 40
158
IF (RELERR .NE. 0.0) THEN
159
RELERR=FOURU*(1.0+FOURU)
160
ABSERR=MAX(ABSERR,0.0D0)
166
IF (.NOT. START) GO TO 300
168
C ***** STARTUP SECTION(FIRST STEP ALONG ZERO CURVE. *****
172
C DETERMINE SUITABLE INITIAL STEP SIZE.
173
H=MIN(H, .10D0, SQRT(SQRT(RELERR*TEMP+ABSERR)))
174
C USE LINEAR PREDICTOR ALONG TANGENT DIRECTION TO START NEWTON ITERATION.
179
CALL TANGNF(S,Y,YP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG,
181
IF (IFLAG .GT. 0) RETURN
183
TEMP=Y(J) + H * YP(J)
189
C CALCULATE THE NEWTON STEP TZ AT THE CURRENT POINT W .
190
CALL TANGNF(RHOLEN,W,WP,YPOLD,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG,
192
IF (IFLAG .GT. 0) RETURN
194
C TAKE NEWTON STEP AND CHECK CONVERGENCE.
199
C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
200
IF (JUDY .EQ. 1) THEN
201
LCALC=DNRM2(NP1,TZ,1)
206
ELSE IF (JUDY .EQ. 2) THEN
207
LCALC=DNRM2(NP1,TZ,1)/LCALC
210
C GO TO MOP-UP SECTION AFTER CONVERGENCE.
211
IF (DNRM2(NP1,TZ,1) .LE. RELERR*DNRM2(NP1,W,1)+ABSERR)
216
C NO CONVERGENCE IN LITFH ITERATIONS. REDUCE H AND TRY AGAIN.
217
IF (H .LE. FOURU*(1.0 + S)) THEN
224
C ***** END OF STARTUP SECTION. *****
226
C ***** PREDICTOR SECTION. *****
229
C COMPUTE POINT PREDICTED BY HERMITE INTERPOLANT. USE STEP SIZE H
230
C COMPUTED ON LAST CALL TO STEPNF .
232
TEMP=QOFS(YOLD(J),YPOLD(J),Y(J),YP(J),HOLD,HOLD+H)
237
C ***** END OF PREDICTOR SECTION. *****
239
C ***** CORRECTOR SECTION. *****
243
C CALCULATE THE NEWTON STEP TZ AT THE CURRENT POINT W .
244
CALL TANGNF(RHOLEN,W,WP,YP,A,QR,ALPHA,TZ,PIVOT,NFE,N,IFLAG,
246
IF (IFLAG .GT. 0) RETURN
248
C TAKE NEWTON STEP AND CHECK CONVERGENCE.
253
C COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION.
254
IF (JUDY .EQ. 1) THEN
255
LCALC=DNRM2(NP1,TZ,1)
260
ELSE IF (JUDY .EQ. 2) THEN
261
LCALC=DNRM2(NP1,TZ,1)/LCALC
264
C GO TO MOP-UP SECTION AFTER CONVERGENCE.
265
IF (DNRM2(NP1,TZ,1) .LE. RELERR*DNRM2(NP1,W,1)+ABSERR)
270
C NO CONVERGENCE IN LITFH ITERATIONS. RECORD FAILURE AT CALCULATED H ,
271
C SAVE THIS STEP SIZE, REDUCE H AND TRY AGAIN.
274
IF (H .LE. FOURU*(1.0 + S)) THEN
281
C ***** END OF CORRECTOR SECTION. *****
283
C ***** MOP-UP SECTION. *****
285
C YOLD AND Y ALWAYS CONTAIN THE LAST TWO POINTS FOUND ON THE ZERO
286
C CURVE OF THE HOMOTOPY MAP. YPOLD AND YP CONTAIN THE TANGENT
287
C VECTORS TO THE ZERO CURVE AT YOLD AND Y , RESPECTIVELY.
300
C ***** END OF MOP-UP SECTION. *****
302
C ***** OPTIMAL STEP SIZE ESTIMATION SECTION. *****
304
C CALCULATE THE DISTANCE FACTOR DCALC .
309
DCALC=DNRM2(NP1,TZ,1)
310
IF (DCALC .NE. 0.0) DCALC=DNRM2(NP1,W,1)/DCALC
312
C THE OPTIMAL STEP SIZE HBAR IS DEFINED BY
314
C HT=HOLD * [MIN(LIDEAL/LCALC, RIDEAL/RCALC, DIDEAL/DCALC)]**(1/P)
316
C HBAR = MIN [ MAX(HT, BMIN*HOLD, HMIN), BMAX*HOLD, HMAX ]
318
C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, SET THE CONTRACTION
319
C FACTOR LCALC TO ZERO.
320
IF (ITNUM .EQ. 1) LCALC = 0.0
321
C FORMULA FOR OPTIMAL STEP SIZE.
322
IF (LCALC+RCALC+DCALC .EQ. 0.0) THEN
325
HT = (1.0/MAX(LCALC/SSPAR(1), RCALC/SSPAR(2), DCALC/SSPAR(3)))
326
& **(1.0/SSPAR(8)) * HOLD
328
C HT CONTAINS THE ESTIMATED OPTIMAL STEP SIZE. NOW PUT IT WITHIN
330
H=MIN(MAX(HT,SSPAR(6)*HOLD,SSPAR(4)), SSPAR(7)*HOLD, SSPAR(5))
331
IF (ITNUM .EQ. 1) THEN
332
C IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, DON'T DECREASE H .
334
ELSE IF (ITNUM .EQ. LITFH) THEN
335
C IF CONVERGENCE REQUIRED THE MAXIMUM LITFH ITERATIONS, DON'T
339
C IF CONVERGENCE DID NOT OCCUR IN LITFH ITERATIONS FOR A PARTICULAR
340
C H = HFAIL , DON'T CHOOSE THE NEW STEP SIZE LARGER THAN HFAIL .
341
IF (FAIL) H=MIN(H,HFAIL)
added
removed
packages/hompack/stepnf.f
