@@ -186,127 +186,6 @@ subroutine fixpnf( &
186186 ! ! \( [\partial \rho/\partial \lambda, \partial \rho/\partial x] \) evaluated
187187 ! ! at \( (a,\lambda,x) \).
188188
189- ! ON INPUT:
190- !
191- ! N is the dimension of X, F(X), and RHO(A,LAMBDA,X).
192- !
193- ! Y(:) is an array of length N + 1. (Y(2),...,Y(N+1)) = A is the
194- ! starting point for the zero curve for the fixed point and
195- ! zero finding problems. (Y(2),...,Y(N+1)) = X0 for the curve
196- ! tracking problem.
197- !
198- ! IFLAG can be -2, -1, 0, 2, or 3. IFLAG should be 0 on the
199- ! first call to FIXPNF for the problem X=F(X), -1 for the
200- ! problem F(X)=0, and -2 for the problem RHO(A,LAMBDA,X)=0.
201- ! In certain situations IFLAG is set to 2 or 3 by FIXPNF,
202- ! and FIXPNF can be called again without changing IFLAG.
203- !
204- ! ARCRE , ARCAE are the relative and absolute errors, respectively,
205- ! allowed the normal flow iteration along the zero curve. If
206- ! ARC?E <= 0.0 on input it is reset to .5*SQRT(ANS?E) .
207- ! Normally ARC?E should be considerably larger than ANS?E .
208- !
209- ! ANSRE , ANSAE are the relative and absolute error values used for
210- ! the answer at LAMBDA = 1. The accepted answer Y = (LAMBDA, X)
211- ! satisfies
212- !
213- ! |Y(1) - 1| <= ANSRE + ANSAE .AND.
214- !
215- ! ||Z|| <= ANSRE*||X|| + ANSAE where
216- !
217- ! (.,Z) is the Newton step to Y.
218- !
219- ! TRACE is an integer specifying the logical I/O unit for
220- ! intermediate output. If TRACE > 0 the points computed on
221- ! the zero curve are written to I/O unit TRACE .
222- !
223- ! A(:) contains the parameter vector A . For the fixed point
224- ! and zero finding problems, A need not be initialized by the
225- ! user, and is assumed to have length N. For the curve
226- ! tracking problem, A must be initialized by the user.
227- !
228- ! SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) is
229- ! a vector of parameters used for the optimal step size estimation.
230- ! If SSPAR(J) <= 0.0 on input, it is reset to a default value
231- ! by FIXPNF . Otherwise the input value of SSPAR(J) is used.
232- ! See the comments below and in STEPNF for more information about
233- ! these constants.
234- !
235- ! POLY_SWITCH is an optional logical variable used only by the driver
236- ! POLSYS1H for polynomial systems.
237- !
238- !
239- ! ON OUTPUT:
240- !
241- ! N , TRACE , A are unchanged.
242- !
243- ! Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, and Y is an approximate
244- ! zero of the homotopy map. Normally LAMBDA = 1 and X is a
245- ! fixed point(zero) of F(X). In abnormal situations LAMBDA
246- ! may only be near 1 and X is near a fixed point(zero).
247- !
248- ! IFLAG =
249- ! -2 causes FIXPNF to initialize everything for the problem
250- ! RHO(A,LAMBDA,X) = 0 (use on first call).
251- !
252- ! -1 causes FIXPNF to initialize everything for the problem
253- ! F(X) = 0 (use on first call).
254- !
255- ! 0 causes FIXPNF to initialize everything for the problem
256- ! X = F(X) (use on first call).
257- !
258- ! 1 Normal return.
259- !
260- ! 2 Specified error tolerance cannot be met. Some or all of
261- ! ARCRE , ARCAE , ANSRE , ANSAE have been increased to
262- ! suitable values. To continue, just call FIXPNF again
263- ! without changing any parameters.
264- !
265- ! 3 STEPNF has been called 1000 times. To continue, call
266- ! FIXPNF again without changing any parameters.
267- !
268- ! 4 Jacobian matrix does not have full rank. The algorithm
269- ! has failed (the zero curve of the homotopy map cannot be
270- ! followed any further).
271- !
272- ! 5 The tracking algorithm has lost the zero curve of the
273- ! homotopy map and is not making progress. The error tolerances
274- ! ARC?E and ANS?E were too lenient. The problem should be
275- ! restarted by calling FIXPNF with smaller error tolerances
276- ! and IFLAG = 0 (-1, -2).
277- !
278- ! 6 The normal flow Newton iteration in STEPNF or ROOTNF
279- ! failed to converge. The error tolerances ANS?E may be too
280- ! stringent.
281- !
282- ! 7 Illegal input parameters, a fatal error.
283- !
284- ! 8 Memory allocation error, fatal.
285- !
286- ! ARCRE , ARCAE , ANSRE , ANSAE are unchanged after a normal return
287- ! (IFLAG = 1). They are increased to appropriate values on the
288- ! return IFLAG = 2 .
289- !
290- ! NFE is the number of function evaluations (= number of
291- ! Jacobian matrix evaluations).
292- !
293- ! ARCLEN is the length of the path followed.
294- !
295- ! Allocatable and automatic work arrays:
296- !
297- ! YP(1:N+1) is a work array containing the tangent vector to
298- ! the zero curve at the current point Y .
299- !
300- ! YOLD(1:N+1) is a work array containing the previous point found
301- ! on the zero curve.
302- !
303- ! YPOLD(1:N+1) is a work array containing the tangent vector to
304- ! the zero curve at YOLD .
305- !
306- ! QR(1:N,1:N+2), ALPHA(1:3*N+3), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1),
307- ! WP(1:N+1), Z0(1:N+1), Z1(1:N+1) are all work arrays used by
308- ! STEPNF to calculate the tangent vectors and Newton steps.
309-
310189 use hompack_kinds, only: zero, one, eps64
311190 implicit none
312191
@@ -335,18 +214,21 @@ subroutine fixpnf( &
335214 ! ! increased and the routine may be called again.
336215 ! ! * `3` : iteration limit reached; call again to continue.
337216 ! ! * `4` : Jacobian matrix lost full rank.
338- ! ! * `5` : tracking algorithm lost the zero curve.
339- ! ! * `6` : Newton iteration failed to converge.
217+ ! ! * `5` : tracking algorithm lost the zero curve and is not making progress; the
218+ ! ! tolerances `arc?e` and `ans?e` were too lenient. Retry with smaller
219+ ! ! tolerances.
220+ ! ! * `6` : normal flow Newton iteration failed to converge; tolerances `ans?e` may be
221+ ! ! too stringent.
340222 ! ! * `7` : illegal input parameters.
341223 ! ! * `8` : memory allocation failure.
342224 real (dp), intent (inout ) :: arcre
343- ! ! Relative error tolerance for the normal-flow iteration used while
344- ! ! tracking the zero curve. If nonpositive on input, a default value is
345- ! ! chosen . May be increased when `iflag=2`.
225+ ! ! Relative error tolerance for the normal-flow iteration used while tracking the
226+ ! ! zero curve. If nonpositive on input, the default value `arcre=0.5*sqrt(ansre)` is
227+ ! ! used . May be increased when `iflag=2`.
346228 real (dp), intent (inout ) :: arcae
347- ! ! Absolute error tolerance for the normal-flow iteration used while
348- ! ! tracking the zero curve. If nonpositive on input, a default value is
349- ! ! chosen . May be increased when `iflag=2`.
229+ ! ! Absolute error tolerance for the normal-flow iteration used while tracking the
230+ ! ! zero curve. If nonpositive on input, the default value `arcae=0.5*sqrt(ansae)` is
231+ ! ! used . May be increased when `iflag=2`.
350232 real (dp), intent (inout ) :: ansre
351233 ! ! Relative error tolerance required of the final solution at `lambda=1`.
352234 ! ! May be increased when `iflag=2`.
@@ -356,14 +238,9 @@ subroutine fixpnf( &
356238 real (dp), intent (inout ), optional :: sspar(8 )
357239 ! ! Step-size control parameters:
358240 ! ! `(lideal, rideal, dideal, hmin, hmax, bmin, bmax, p)`.
359- ! ! Parameters used by the optimal step-size estimation algorithm.
360241 ! ! Elements that are nonpositive on input are replaced by default values.
361242 real (dp), intent (in ), optional :: a(:)
362- ! ! Parameter vector `a`.
363- ! ! For fixed-point and zero-finding problems, the array is assumed to have
364- ! ! length `n` and need not be initialized by the user.
365- ! ! For curve-tracking problems, it must be initialized on input.
366- ! ! Unchanged on output.
243+ ! ! Parameter vector `a` for curve-tracking problems.
367244 integer , intent (in ), optional :: lunit
368245 ! ! Logical I/O unit for intermediate output.
369246 ! ! * `0` : No output is printed (default).
@@ -797,9 +674,8 @@ subroutine stepnf(state, callbacks)
797674 integer :: itnum, j, judy, np1
798675 logical :: fail
799676
800- ! The limit on the number of Newton iterations allowed before reducing the step
801- ! size 'h' may be changed by changing the following parameter
802- integer , parameter :: litfh = 4
677+ ! Newton iterations allowed before reducing the step size 'h'
678+ integer , parameter :: max_newton = 4
803679
804680 np1 = state% n + 1
805681 state% crash = .true.
@@ -851,7 +727,7 @@ subroutine stepnf(state, callbacks)
851727 ws% w = state% y + state% h* state% yp
852728 ws% z0 = ws% w
853729
854- do judy = 1 , litfh
730+ do judy = 1 , max_newton
855731
856732 ! Calculate the Newton step 'tz' at the current point 'w'
857733 rholen = - one
@@ -882,7 +758,7 @@ subroutine stepnf(state, callbacks)
882758
883759 end do
884760
885- ! No convergence in litfh iterations. Reduce h and try again.
761+ ! No convergence in 'max_newton' iterations. Reduce h and try again.
886762 if (state% h <= fouru* (one + state% s)) then
887763 state% iflag = 6
888764 return
@@ -909,7 +785,7 @@ subroutine stepnf(state, callbacks)
909785 ws% z0 = ws% w
910786
911787 ! CORRECTOR SECTION
912- do judy = 1 , litfh
788+ do judy = 1 , max_newton
913789
914790 ! Calculate the Newton step 'tz' at the current point 'w'
915791 rholen = - one
@@ -939,7 +815,7 @@ subroutine stepnf(state, callbacks)
939815
940816 end do
941817
942- ! No convergence in 'litfh ' iterations. Record failure at calculated 'h'
818+ ! No convergence in 'max_newton ' iterations. Record failure at calculated 'h'
943819 ! Save this step size, reduce 'h' and try again.
944820
945821 fail = .true.
@@ -1005,13 +881,13 @@ subroutine stepnf(state, callbacks)
1005881 if (itnum == 1 ) then
1006882 ! If convergence had occurred after 1 iteration, don't decrease 'h'.
1007883 state% h = max (state% h, state% hold)
1008- else if (itnum == litfh ) then
1009- ! If convergence required the maximum 'litfh ' iterations, don't increase 'h'.
884+ else if (itnum == max_newton ) then
885+ ! If convergence required the maximum 'max_newton ' iterations, don't increase 'h'.
1010886 state% h = min (state% h, state% hold)
1011887 end if
1012888
1013- ! If convergence did not occur in 'litfh ' iterations for a particular 'h = hfail',
1014- ! don't choose the new step size larger than 'hfail'.
889+ ! If convergence did not occur in 'max_newton ' iterations for a particular
890+ ! 'h = hfail', don't choose the new step size larger than 'hfail'.
1015891 if (fail) state% h = min (state% h, hfail)
1016892
1017893 end associate
0 commit comments