pitcon.f90

求解非线性方程组的一组源代码,FORTRAN90.用于解决N个未知数,N-1个方程.

!         and the user routines FX and DF.  IPAR is declared in the PITCON
!         program and passed through it to FX and DF, but otherwise ignored.
!         Note, however, that if BANSLV is used for the solver routine, then
!         IPAR(1) must contain the lower bandwidth, and IPAR(2) the upper
!         bandwidth of the Jacobian matrix.
!
!  IWORK  Input/output, integer IWORK(LIW).  Communication and workspace array.
!
!         The specific allocation of IWORK is described in the section devoted
!         to the work arrays.  Some elements of IWORK MUST be set by the user,
!         others may be set to change the way PITCON works.
!
!  LIW    Input, integer LIW, the dimension of IWORK.
!
!         The minimum acceptable value of LIW depends on the solver chosen,
!         but for either DENSLV or BANSLV, setting LIW = 29+NVAR is sufficient.
!
!  NVAR   Input, integer NVAR, the number of variables, the dimension of X.
!
!         This is, of course, one greater than the number of equations or
!         functions.  NVAR must be at least 2.
!
!  RWORK  Input/output, real RWORK(LRW), work array.
!
!         The specific allocation of RWORK is described in the section
!         devoted to the work arrays.
!
!  LRW    Input, integer LRW, the dimension of RWORK.
!
!         The minimum acceptable value depends heavily on the solver options.
!         There is storage required for scalars, vectors, and the Jacobian
!         array.  The minimum acceptable value of LRW is the sum of three
!         corresponding numbers.
!
!         For DENSLV with user-supplied Jacobian,
!
!           LRW = 29 + 4*NVAR + NVAR*NVAR.
!
!         For DENSLV with internally approximated Jacobian,
!
!           LRW = 29 + 6*NVAR + NVAR*NVAR.
!
!         For BANSLV, with a Jacobian matrix with upper bandwidth MU and lower
!         bandwidth ML, and NBAND = 2*ML+MU+1, with user supplied Jacobian,
!
!           LRW = 29 + 6*NVAR + (NVAR-1)*NBAND.
!
!         For BANSLV with internally approximated Jacobian,
!
!           LRW = 29 + 9*NVAR + (NVAR-1)*NBAND.
!
!  XR     Input/output, real XR(NVAR), the current solution point.
!
!         On the first call, the user should set XR to a starting point which
!         at least approximately satisfies F(XR) = 0.  The user need never
!         update XR again.
!
!         Thereafter, on each return from the program with IERROR = 0, XR will
!         hold the most recently computed point, whether a continuation, target
!         or limit point.
!
!  SLVNAM Input, external SLVNAM, the name of the solver to use on linear
!         systems.
!
!         The linear systems have the form A*x = b, where A is the augmented
!         Jacobian matrix.  A will be square, and presumably nonsingular.
!         The routine must return the value of the solution x.
!
!         Two possible choices for SLVNAM are "DENSLV" and "BANSLV", which are
!         the names of routines provided with the package.  DENSLV is
!         appropriate for a full storage jacobian, and BANSLV for a jacobian
!         which is banded except for the last column.
!
!         The advanced user may study the source code for these two routines
!         and write an equivalent solver more suited to a given problem.
!
!
!  G) The Integer Work Array IWORK:
!  -------------------------------
!
!  Input to PITCON includes the setting of some of the entries in IWORK.
!  Some of this input is optional.  The user input section of IWORK involves
!  entries 1 through 9, and, possibly also 17.
!
!  IWORK(1) must be set by the user.  All other entries have default values.
!
!
!  IWORK(1)        On first call only, the user must set IWORK(1) = 0.
!                  Thereafter, the program sets IWORK(1) before return to
!                  explain what kind of point is being returned.  This return
!                  code is:
!
!                      1 return of corrected starting point.
!                      2 return of continuation point.
!                      3 return of target point.
!                      4 return of limit point.
!
!                  NOTE:  At any time, PITCON may be called with a negative
!                  value of IWORK(1). This requests a check of the
!                  jacobian routine against a finite difference approximation,
!                  or a printout of the jacobian or its approximation.
!
!                      -1, compare Jacobian against forward difference,
!                          print maximum discrepancy only.
!                      -2, compare Jacobian against central difference,
!                          print maximum discrepancy only.
!                      -3, compare Jacobian against forward difference,
!                          print out the entire discrepancy matrix.
!                      -4, compare Jacobian against central difference,
!                          print out the entire discrepancy matrix.
!                      -5, compute forward difference Jacobian,
!                          print maximum entry only.
!                      -6, compute central difference Jacobian,
!                          print maximum entry only.
!                      -7, compute forward difference Jacobian,
!                          print entire matrix.
!                      -8, compute central difference Jacobian,
!                          print entire matrix.
!                      -9, request user supplied Jacobian,
!                          print maximum entry only.
!                     -10, request user supplied Jacobian,
!                          print entire matrix.
!
!                  Before a call with negative IWORK(1), the current value of
!                  IWORK(1) should be saved, and then restored to the previous
!                  value after the call, in order to resume calculation.
!
!                  IWORK(1) does not have a default value.  The user MUST set
!                  it.
!
!  IWORK(2)        IPC, the component of the current continuation point XR
!                  which is to be used as the continuation parameter.  On first
!                  call, the program is willing to use the index NVAR as a
!                  default, but the user should set this value if better
!                  information is available.
!
!                  After the first call, the program sets this value for each
!                  step automatically unless the user prevents this by setting
!                  the parameterization option IWORK(3) to a non-zero valus.
!                  Note that a poor choice of IWORK(2) may cause the algorithm
!                  to fail.  IWORK(2) defaults to NVAR on the first step.
!
!  IWORK(3)        Parameterization option.  The program would prefer to be
!                  free to choose a new local parameter from step to step.
!                  The value of IWORK(3) allows or prohibits this action.
!                  IWORK(3) = 0 allows the program to vary the parameter,
!                  IWORK(3) = 1 forces the program to use whatever the contents
!                  of IWORK(2) are, which will not be changed from the user's
!                  input or the default.  The default is IWORK(3) = 0.
!
!  IWORK(4)        Newton iteration Jacobian update option.
!                  0, the Jacobian is reevaluated at every step of the
!                     Newton iteration.  This is costly, but may result in
!                     fewer Newton steps and fewer Newton iteration rejections.
!                  1, the Jacobian is evaluated only on the first and
!                     IWORK(17)-th steps of the Newton process.
!                  2, the Jacobian is evaluated only when absolutely
!                     necessary, namely, at the first step, and when the
!                     process fails. This option is most suitable for problems
!                     with mild nonlinearities.
!
!                  The default is IWORK(4) = 0.
!
!  IWORK(5)        IT, target point index.  If IWORK(5) is not zero, it is
!                  presumed to be the component index between 1 and NVAR for
!                  which target points are sought.  In this case, the value of
!                  RWORK(7) is assumed to be the target value.  The program
!                  will monitor every new continuation point, and if it finds
!                  that a target point may lie between the new point and the
!                  previous point, it will compute this target point and
!                  return.  This target point is defined by the property that
!                  its component with the index prescribed in IWORK(5) will
!                  have the value given in RWORK(7).  For a given problem there
!                  may be zero, one, or many target points requested.
!                  The default of IWORK(5) is 0.
!
!  IWORK(6)        LIM, the limit point index.  If IWORK(6) is nonzero, then
!                  the program will search for limit points with respect to
!                  the component with index IWORK(6); that is, of points for
!                  which the IWORK(6)-th variable has a local extremum, or
!                  equivalently where the IWORK(6)-th component of the tangent
!                  vector is zero.  The default of IWORK(6) is zero.
!
!  IWORK(7)        IWRITE, which controls the amount of output produced by the
!                  program. IWORK(7) may have a value between 0 and 3.
!                  For IWORK(7) = 0 there is almost no output while for
!                  IWORK(7) = 3 the most output is produced.
!                  The default is 1.
!
!  IWORK(9)        Control of the Jacobian option specifying whether the user
!                  has supplied a Jacobian routine, or wants the program
!                  to approximate the Jacobian.
!                  0, the user has supplied the Jacobian.
!                  1, program is to use forward difference approximation.
!                  2, program is to use central difference approximation.
!                  IWORK(9) defaults to 0.
!
!  IWORK(10)       State indicator of the progress of the program.
!                  The values are:
!                  0, start up with unchecked starting point.
!                  1, first step.  Corrected starting point available.
!                  2, two successive continuation points available, as well
!                     as the tangent vector at the oldest of them.
!                  3, two successive continuation points available, as well
!                     as the tangent vector at the newest of them.
!
!  IWORK(11)       Index of the last computed target point. This is used to
!                  avoid repeated computation of a target point.  If a target
!                  point has been found, then the target index IWORK(5) is
!                  copied into IWORK(11).
!
!  IWORK(12)       Second best choice for the local parameterization index.
!                  This index may be tried if the first choice causes poor
!                  performance in the Newton corrector.
!
!  IWORK(13)       Beginning location in IWORK of unused integer work space
!                  available for use by the solver.
!
!  IWORK(14)       LIW, the user declared dimension of the array IWORK.
!
!  IWORK(15)       Beginning location in RWORK of unused real work space
!                  available for use by the solver.
!
!  IWORK(16)       LRW, the user declared dimension of RWORK.
!
!  IWORK(17)       Maximum number of corrector steps allowed during one run
!                  of the Newton process in which the Jacobian is updated at
!                  every step.  If the Jacobian is only evaluated at
!                  the beginning of the Newton iteration then 2*IWORK(17) steps
!                  are allowed.
!                  IWORK(17) must be greater than 0.  It defaults to 10.
!
!  IWORK(18)       Number of stepsize reductions that were needed for
!                  producing the last continuation point.
!
!  IWORK(19)       Total number of calls to the user Jacobian routine DF.
!
!  IWORK(20)       Total number of calls to the matrix factorization routine.
!                  If DENSLV is the chose solver then factorization is done by
!                  the LINPACK routine SGEFA.  If BANSLV is the solver, the
!                  LINPACK routine SGBFA will be used.
!
!  IWORK(21)       Total number of calls to the back-substitution routine.
!                  If DENSLV is the chosen solver, then back substitution is
!                  done by the LINPACK routine SGESL.  If BANSLV is used, then
!                  the LINPACK routine SGBSL will be used.
!
!  IWORK(22)       Total number of calls to the user function routine FX.
!
!  IWORK(23)       Total number of steps taken in limit point iterations.
!                  Each step involves determining an approximate limit point
!                  and applying a Newton iteration to correct it.
!
!  IWORK(24)       Total number of Newton corrector steps used during the
!                  computation of target points.
!
!  IWORK(25)       Total number of Newton steps taken during the correction
!                  of a starting point or the continuation points.
!
!  IWORK(26)       Total number of predictor stepsize-reductions needed
!                  since the start of the continuation procesds.
!
!  IWORK(27)       Total number of calls to the program.  This also
!                  corresponds to the number of points computed.
!
!  IWORK(28)       Total number of Newton steps taken during current iteration.
!
!  IWORK(30)       and on are reserved for use by the linear equation solver,
!                  and typically are used for pivoting.
!
!
!  H) The Real Work Array RWORK:
!  ----------------------------
!
!  Input to PITCON includes the setting of some of the entries in RWORK.
!  Some of this input is optional.  The user input section of RWORK involves
!  entries 1 through 7 and possibly 18 and 20.
!
!  All entries of RWORK have default values.
!
!
!  RWORK(1)        ABSERR, absolute error tolerance.   This value is used
!                  mainly during the Newton iteration.  RWORK(1) defaults to
!                  SQRT(EPMACH) where EPMACH is the machine relative precision
!                  stored in RWORK(8).
!
!  RWORK(2)        RELERR, relative error tolerance.  This value is used mainly
!                  during the Newton iteration.  RWORK(2) defaults to
!                  SQRT(EPMACH) where EPMACH is the machine relative precision
!                  stored in RWORK(8).
!
!  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