!
!  RWORK(3)        HMIN, minimum allowable predictor stepsize.  If failures of
!                  the Newton correction force the stepsize down to this level,
!                  then the program will give up.  The default value is
!                  SQRT(EPMACH).
!
!  RWORK(4)        HMAX, maximum allowable predictor step.  Too generous a value
!                  may cause erratic behavior of the program.  The default
!                  value is SQRT(NVAR).
!
!  RWORK(5)        HTAN,  the predictor stepsize.  On first call, it should be
!                  set by the user.  Thereafter it is set by the program.
!                  RWORK(5) should be positive.  In order to travel in the
!                  negative direction, see RWORK(6).
!                  The default initial value equals 0.5*(RWORK(3)+RWORK(4)).
!
!  RWORK(6)        The local continuation direction, which is either +1.0
!                  or -1.0 .  This asserts that the program is moving in the
!                  direction of increasing or decreasing values of the local
!                  continuation variable, whose index is in IWORK(2).  On first
!                  call, the user must choose IWORK(2).  Therefore, by setting
!                  RWORK(6), the user may also specify whether the program is
!                  to move initially to increase or decrease the variable whose
!                  index is IWORK(2).
!                  RWORK(6) defaults to +1.
!
!  RWORK(7)        A target value.  It is only used if a target index
!                  has been specified through IWORK(5).  In that case, solution
!                  points with the IWORK(5) component equal to RWORK(7) are
!                  to be computed. The code will return each time it finds such
!                  a point.  RWORK(7) does not have a default value.  The
!                  program does not set it, and it is not referenced unless
!                  IWORK(5) has been set.
!
!  RWORK(8)        EPMACH, the value of the machine precision.  The computer
!                  can distinguish 1.0+EPMACH from 1.0, but it cannot
!                  distinguish 1.0+(EPMACH/2) from 1.0. This number is used
!                  when estimating a reasonable accuracy request on a given
!                  computer.  PITCON computes a value for EPMACH internally.
!
!  RWORK(9)        STEPX, the size, using the maximum-norm, of the last
!                  step of Newton correction used on the most recently
!                  computed point, whether a starting point, continuation
!                  point, limit point or target point.
!
!  RWORK(10)       A minimum angle used in the steplength computation,
!                  equal to 2.0*ARCCOS(1-EPMACH).
!
!  RWORK(11)       Estimate of the angle between the tangent vectors at the
!                  last two continuation points.
!
!  RWORK(12)       The pseudo-arclength coordinate of the previous continuation
!                  pointl; that is, the sum of the Euclidean distances between
!                  all computed continuation points beginning with the start
!                  point.  Thus each new point should have a larger coordinate,
!                  except for target and limit points which lie between the two
!                  most recent continuation points.
!
!  RWORK(13)       Estimate of the pseudo-arclength coordinate of the current
!                  continuation point.
!
!  RWORK(14)       Estimate of the pseudoarclength coordinate of the current
!                  limit or target point, if any.
!
!  RWORK(15)       Size of the correction of the most recent continuation
!                  point; that is, the maximum norm of the distance between the
!                  predicted point and the accepted corrected point.
!
!  RWORK(16)       Estimate of the curvature between the last two
!                  continuation points.
!
!  RWORK(17)       Sign of the determinant of the augmented matrix at the
!                  last continuation point whose tangent vector has been
!                  calculated.
!
!  RWORK(18)       This quantity is only used if the jacobian matrix is to
!                  be estimated using finite differences.  In that case,
!                  this value determines the size of the increments and
!                  decrements made to the current solution values, according
!                  to the following formula:
!
!                    Delta X(J) = RWORK(18) * (1.0 + ABS(X(J))).
!
!                  The value of every entry of the approximate jacobian could
!                  be extremely sensitive to RWORK(18).  Obviously, no value
!                  is perfect.  Values too small will surely cause singular
!                  jacobians, and values too large will surely cause inaccuracy.
!                  Little more is certain.  However, for many purposes, it
!                  is suitable to set RWORK(18) to the square root of the
!                  machine epsilon, or perhaps to the third or fourth root,
!                  if singularity seems to be occuring.
!
!                  RWORK(18) defaults to SQRT(SQRT(RWORK(8))).
!
!                  The user may set RWORK(18).  If it has a nonzero value on
!                  input, that value will be used.  Otherwise, the default
!                  value will be used.
!
!  RWORK(19)       Not currently used.
!
!  RWORK(20)       Maximum growth factor for the predictor stepsize based
!                  on the previous secant stepsize.  The stepsize algorithm
!                  will produce a suggested step that is no less that the
!                  previous secant step divided by this factor, and no greater
!                  than the previous secant step multiplied by that factor.
!                  RWORK(20) defaults to 3.
!
!  RWORK(21)       The (Euclidean) secant distance between the last two
!                  computed continuation points.
!
!  RWORK(22)       The previous value of RWORK(21).
!
!  RWORK(23)       A number judging the quality of the Newton corrector
!                  convergence at the last continuation point.
!
!  RWORK(24)       Value of the component of the current tangent vector
!                  corresponding to the current continuation index.
!
!  RWORK(25)       Value of the component of the previous tangent vector
!                  corresponding to the current continuation index.
!
!  RWORK(26)       Value of the component of the current tangent vector
!                  corresponding to the limit index in IWORK(6).
!
!  RWORK(27)       Value of the component of the previous tangent vector
!                  corresponding to the limit index in IWORK(6).
!
!  RWORK(28)       Value of RWORK(7) when the last target point was
!                  computed.
!
!  RWORK(29)       Sign of the determinant of the augmented matrix at the
!                  previous continuation point whose tangent vector has been
!                  calculated.
!
!  RWORK(30)       through RWORK(30+4*NVAR-1) are used by the program to hold
!                  an old and new continuation point, a tangent vector and a
!                  work vector.  Subsequent entries of RWORK are used by the
!                  linear solver.
!
!
!  I) Programming Notes:
!  --------------------
!
!  The minimal input and user routines required to apply the program are:
!
!    Write a function routine FX of the form described above.
!    Use DENSLV as the linear equation solver by setting SLVNAM to DENSLV.
!    Skip writing a Jacobian routine by using the finite difference option.
!    Pass the name of FX as the Jacobian name as well.
!    Declare the name of the function FX as external.
!    Set NVAR in accordance with your problem.
!
!  Then:
!
!    Dimension the vector IWORK to the size LIW = 29+NVAR.
!    Dimension the vector RWORK to the size LRW = 29+NVAR*(NVAR+6).
!    Dimension IPAR(1) and FPAR(1) as required in the function routine.
!    Dimension XR(NVAR) and set it to an approximate solution of F(XR) = 0.
!
!  Set the work arrays as follows:
!
!    Initialize IWORK to 0 and RWORK to 0.0.
!
!    Set IWORK(1) = 0 (Problem startup)
!    Set IWORK(7) = 3 (Maximum internally generated output)
!    Set IWORK(9) = 1 (Forward difference Jacobian)
!
!  Now call the program repeatedly, and never change any of its arguments.
!  Check IERROR to decide whether the code is working satisfactorily.
!  Print out the vector XR to see the current solution point.
!
!
!  The most obvious input to try to set appropriately after some practice
!  would be the error tolerances RWORK(1) and RWORK(2), the minimum, maximum
!  and initial stepsizes RWORK(3), RWORK(4) and RWORK(5), and the initial
!  continuation index IWORK(2).
!
!  For speed and efficiency, a Jacobian routine should be written. It can be
!  checked by comparing its results with the finite difference Jacobian.
!
!  For a particular problem, the target and limit point input can be very
!  useful.  For instance, in the case of a discretized boundary value problem
!  with a real parameter it may be desirable to compare the computed solutions
!  for different discretization-dimensions and equal values of the parameter.
!  For this the target option can be used with the prescribed values of the
!  parameter. Limit points usually are of importance in connection with
!  stability considerations.
!
!
!  The routine REPS attempts to compute the machine precision, which in practice
!  is simply the smallest power of two that can be added to 1 to produce a
!  result greater than 1.  If the REPS routine misbehaves, you can replace
!  it by code that assigns a constant precomputed value, or by a call to
!  the PORT/SLATEC routine R1MACH.  REPS is called just once, in the
!  PITCON routine, and the value returned is stored into RWORK(8).
!
!  In subroutines DENSLV and BANSLV, the parameter "EPS" is set to RWORK(18)
!  and used in estimating jacobian matrices via finite differences.  If EPS is
!  too large, the jacobian will be very inaccurate.  Unfortunately, if EPS is
!  too small, the "finite" differences may become "infinitesmal" differences.
!  That is, the difference of two function values at close points may be
!  zero.  This is a very common problem, and occurs even with a function
!  like F(X) = X*X.  A singular jacobian is much worse than an inaccurate one,
!  so we have tried setting the default value of RWORK(18) to
!  SQRT(SQRT(RWORK(8)).  Such a value attempts to ward off singularity at the
!  expense of accuracy.  You may find for a particular problem or machine that
!  this value is too large and should be adjusted.  It is an utterly arbitrary
!  value.
!
!
!  J) References:
!  -------------
!
!  1.
!  Werner Rheinboldt,
!  Solution Field of Nonlinear Equations and Continuation Methods,
!  SIAM Journal of Numerical Analysis,
!  Volume 17, 1980, pages 221-237.
!
!  2.
!  Cor den Heijer and Werner Rheinboldt,
!  On Steplength Algorithms for a Class of Continuation Methods,
!  SIAM Journal of Numerical Analysis,
!  Volume 18, 1981, pages 925-947.
!
!  3.
!  Werner Rheinboldt,
!  Numerical Analysis of Parametrized Nonlinear Equations
!  John Wiley and Sons, New York, 1986
!
!  4.
!  Werner Rheinboldt and John Burkardt,
!  A Locally Parameterized Continuation Process,
!  ACM Transactions on Mathematical Software,
!  Volume 9, Number 2, June 1983, pages 215-235.
!
!  5.
!  Werner Rheinboldt and John Burkardt,
!  Algorithm 596, A Program for a Locally Parameterized Continuation Process,
!  ACM Transactions on Mathematical Software,
!  Volume 9, Number 2, June 1983, Pages 236-241.
!
!  6.
!  J J Dongarra, J R Bunch, C B Moler and G W Stewart,
!  LINPACK User's Guide,
!  Society for Industrial and Applied Mathematics,
!  Philadelphia, 1979.
!
!  7.
!  Richard Brent,
!  Algorithms for Minimization without Derivatives,
!  Prentice Hall, 1973.
!
!  8.
!  Tony Chan,
!  Deflated Decomposition of Solutions of Nearly Singular Systems,
!  Technical Report 225,
!  Computer Science Department,
!  Yale University,
!  New Haven, Connecticut, 06520,
!  1982.
!
!
!  L)  Sample programs:
!  -------------------
!
!
!  A number of sample problems are included with PITCON.  There are several
!  examples of the Freudenstein-Roth function, which is a nice example
!  with only a few variables, and nice whole number starting and stopping
!  points.  Other sample problems demonstrate or test various options in
!  PITCON.
!
!
!  1:
!  pcprb1.f
!  The Freudenstein-Roth function.  3 variables.
!
!  This is a simple problem, whose solution curve has some severe bends.
!  This file solves the problem with a minimum of fuss.
!
!
!  2:
!  pcprb2.f
!  The Aircraft Stability problem.  8 variables.
!
!  This is a mildly nonlinear problem, whose solution curve has some
!  limit points that are difficult to track.
!
!
!  3:
!  pcprb3.f
!  A boundary value problem.  22 variables.
!
!  This problem has a limit point in the LAMBDA parameter, which we seek.  We
!  solve this problem 6 times, illustrating the use of full and banded
!  jacobians, and of user-generated, or forward or central difference
!  approximated jacobian matrices.
!
!  The first 21 variables represent the values of a function along a grid
!  of 21 points between 0 and 1.  By increasing the number of grid points,
!  the problem can be set up with arbitrarily many variables.  This change
!  requires changing a single parameter in the main program.
!
!