pitcon.f90

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

        write ( *, * ) '  target point computation.'
        write ( *, * ) ' '
        write ( *, * ) '  The target point returned does'
        write ( *, * ) '  not satisfy the accuracy requirements.'
        write ( *, * ) ' '
        write ( *, * ) '  However, the code can continue.'
      end if
    end if
    return
  end if
!
!*******************************************************************************
!
!  6.  Tangent and local continuation parameter calculation.
!
!  Unless the tangent and limit point calculations were already
!  performed, (because the loop was interrupted for a limit point
!  calculation), set up and solve the equation for the tangent vector.
!
!  Force the tangent vector to be of unit length, and try to preserve
!  the "sense" or "direction" of the curve by forcing the IPL-th
!  component of the tangent vector to agree in sign with the IPL-th
!  component of the previous secant vector.  (On the first step, however,
!  we have to use the user's input direction to choose the sign).
!
!  Set the local continuation parameter IPC.
!
!  If IWORK(3) is 0, the program is free to vary IPC from step to step.
!  In that case, IPC is normally set to the index of the component of
!  greatest magnitude in the tangent vector.  However, if a limit point
!  appears to be coming in that direction, the index of the second
!  greatest magnitude component might be chosen instead.
!
!*******************************************************************************
!
  if ( iwork(10) /= 3 ) then

    call tanpar ( df,fpar,fx,ierror,ipar,iwork,liw,lrw,nvar,rwork,slname, &
      rwork(ltc),rwork(lwk),rwork(lxc),rwork(lxf),xr)

    if ( ierror/= 0 ) then
      if ( iwrite>0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'PITCON - Fatal error.'
        write ( *, * ) '  The computation failed while computing'
        write ( *, * ) '  the parameter and the tangent vector.'
        write ( *, * ) ' '
        write ( *, * ) '  The program can not proceed!'
      end if
      return
    end if
!
!*******************************************************************************
!
!  7.  Limit point check.
!
!  Skip this section if IWORK(6) = 0.
!
!  Otherwise, user has requested a search for limit points in a given
!  index by setting IWORK(6), also called "LIM", to a nonzero value.
!
!  Compare LIM-th components of previous and current tangent vectors.
!  If a sign difference occurs, we assume a limit point has been passed.
!  Attempt to compute a point XR between the previous and current points,
!  for which the LIM-th component of the tangent is zero.
!
!  This search will be guided by a rootfinder.  The scalar function
!  to be zeroed out is the LIM-th tangent vector component.
!
!*******************************************************************************
!
    if ( (iwork(6)/= 0) .and. (iwork(1)/=4).and. (iwork(10) == 3).and. &
       (sign(one,rwork(26))/= sign(one,rwork(27))) ) then

      call limit(df,fpar,fx,ierror,ipar,iwork,liw,lrw,nvar,rwork, &
        slname,rwork(ltc),rwork(lwk),rwork(lxc),rwork(lxf),xr)

      if ( ierror/= 0 ) then
        if ( iwrite>0 ) then
          write ( *, * ) ' '
          write ( *, * ) 'PITCON - Warning!'
          write ( *, * ) '  An error occurred during the'
          write ( *, * ) '  limit point computation.'
          write ( *, * ) ' '
          write ( *, * ) '  The computed limit point does not'
          write ( *, * ) '  satisfy the accuracy requirements.'
          write ( *, * ) ' '
          write ( *, * ) '  However, the code can continue.'
        end if
      end if

      return

    end if

  end if
!
!*******************************************************************************
!
!  8.  Compute next predictor step length, HTAN.
!
!*******************************************************************************
!
  if ( iwork(10) > 1 ) then
    call setstp ( iwork, liw, lrw, rwork )
  end if
!
!*******************************************************************************
!
!  9.  Continuation step
!
!  Our current data is the current point XC, its tangent vector TC, and
!  a steplength HTAN.  We predict the location of the next point on the
!  curve using the Euler approximation XR = XC+HTAN*TC.
!
!  Newton iteration is applied to this point, to force it to lie on the
!  curve.  In order to make the system square, an augmenting equation
!  is added to the system, specifying that XR(IPC) = XC(IPC)+HTAN*TC(IPC).
!  (The right hand side is a constant.)
!
!  If the Newton correction process fails, the stepsize is reduced and
!  prediction and correction retried.  Failure will most likely be
!  signaled by repeated step reductions, until the minimum allowable
!  stepsize is reached.  If this occurs, PITCON has failed, and cannot
!  proceed along the curve any more.
!
!*******************************************************************************
!
  call trystp ( df, fpar, fx, ierror, ipar, iwork, liw, lrw, nvar, rwork, &
    slname, rwork(ltc), rwork(lwk), rwork(lxf), xr )

  if ( ierror /= 0 ) then
    if ( iwrite > 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'PITCON - Fatal error.'
      write ( *, * ) 'The computation failed while trying'
      write ( *, * ) 'to compute the next point.'
      write ( *, * ) ' '
      write ( *, * ) 'The program can not proceed!'
    end if
    return
  end if
!
!*******************************************************************************
!
!  10.  Successful step.  Update information.
!
!*******************************************************************************
!
  call update ( iwork, liw, lrw, nvar, rwork, rwork(ltc), rwork(lxc), &
    rwork(lxf), xr )
!
!  Compute the convergence "quality", a factor between 1/8 and 8, which
!  tries to estimate how deeply inside the Newton attraction region we
!  are, and hence how much bolder or more timid we could be on the next
!  prediction.
!
  modnew = iwork(4)

  call coqual ( iwrite, modnew, qual, iwork, liw, rwork, lrw )

  rwork(23) = qual

  return
end
subroutine dgb_jac ( eps, fcol, fpar, fprime, frow, fx, ipar, ipc, iwork, &
  jac, liw, nband, nvar, x )
!
!******************************************************************************
!
!! DGB_JAC approximates a banded jacobian matrix.
!
!
!  Discussion:
!
!    DGB_JAC estimates the jacobian matrix FPRIME of the function FX,
!    using forward or central finite differences.  DGB_JAC is called by
!    DGB_SLV when the user has specified the jacobian option as 1 or 2.
!
!  Parameters:
!
!    Input, real EPS, a tolerance used for shifting the X values.
!    A value of the square root of the machine precision is
!    usually appropriate.
!
!    Output, real FCOL(NEQN), the last column of the approximate
!    jacobian, which is allowed to be "full".  This comprises
!    matrix entries FPRIME(1,NVAR) through FPRIME(NEQN,NVAR).
!
!    Input/output, real FPAR(*), user parameter vector, not
!    touched by this routine, but passed on to user routines.
!
!    Output, real FPRIME(NBAND,NEQN), is the array into which the
!    the banded portion of the computed jacobian will be stored.
!    The LAPACK general band format is used, assigning entry (I,J)
!    to FPRIME(I-J+ML+MU+1,J), where ML and MU are the lower and
!    upper half bandwidths respectively.
!
!    Output, real FROW(NVAR), storage for the last (augmenting) row
!    of the jacobian, which will be all zero except for a 1 in
!    location IPC.
!
!    Input, external FX, the name of the routine which evaluates the
!    function.
!
!    FX computes the value of the nonlinear function.  This name
!    must be declared external in the calling program.  FX should
!    evaluate the NVAR-1 function components at the input point X,
!    and store the result in the vector FVEC.  An augmenting
!    equation will be stored in entry NVAR of FVEC by the PITCON
!    program.
!
!    FX should have the following form:
!
!    subroutine fx ( nvar, fpar, ipar, x, fvec )
!
!      Input, integer NVAR, number of variables.
!
!      Input/output, real FPAR(*), array of user parameters.
!
!      Input/output, integer IPAR(*), array of user parameters.
!
!      Input, real X(NVAR), the point at which function evaluation is required.
!
!      Output, real FVEC(NVAR), the value of the function at point
!      X.  Only the first NVAR-1 entries of FVEC are to be set by
!      the routine.  PITCON sets the final value itself.
!
!    Input, integer IPAR(*), a user parameter vector passed to FX.
!    However, because this is a problem with a banded jacobian, entries
!    IPAR(1) and IPAR(2) are read by this routine.  IPAR(1) contains
!    ML, the lower half bandwidth of the jacobian, and IPAR(2) contains
!    MU, the upper half bandwidth of the jacobian.
!
!    Input, integer IPC, the index of the current continuation parameter,
!    which is needed to determine the form of FROW.
!
!    Input, integer IWORK(LIW), work and statistics vector.  Only
!    required here so that we can count the number of function
!    evaluations.
!
!    Input, integer JAC, the user requested jacobian option.  For
!    our purposes, the only two values of interest are:
!
!      1 = estimate jacobian with forward differences,
!      2 = estimate jacobian with central differences (twice the work)
!
!    Input, integer LIW, the dimension of IWORK.
!
!    Input, integer NBAND, the first dimension of the jacobian matrix
!    FPRIME, NBAND = ML+MU+1.
!
!    Input, integer NVAR, the number of variables.
!
!    Input, real X(NVAR), the point at which the jacobian is desired.
!
  external fx
!
  integer liw
  integer nband
  integer nvar
!
  double precision eps
  double precision fcol(nvar-1)
  double precision fpar(*)
  double precision fprime(nband,nvar-1)
  double precision frow(nvar)
  integer i
  integer ihi
  integer ilo
  integer ipar(2)
  integer ipc
  integer iwork(liw)
  integer j
  integer jac
  integer kcall
  integer mband
  integer ml
  integer mu
  double precision skale
  double precision x(nvar)
  double precision xjac
  double precision xtemp(nvar)
  double precision work1(nvar)
  double precision work2(nvar)
!
  ml = ipar(1)
  mu = ipar(2)
  mband = ml+mu+1

  if ( jac == 1 ) then
    call fx ( nvar, fpar, ipar, x, work2 )
    iwork(22) = iwork(22)+1
  end if

  if ( jac == 2 ) then
    xjac = 2.0
  else
    xjac = 1.0
  end if

  do kcall = 1, mband

    xtemp(1:nvar) = x(1:nvar)

    do j = kcall, nvar-1, mband
      xtemp(j) = x(j)+eps*(1.0+abs(x(j)))
    end do

    call fx ( nvar, fpar, ipar, xtemp, work1 )
    iwork(22) = iwork(22)+1

    if ( jac == 2 ) then

      xtemp(1:nvar) = x(1:nvar)

      do j = kcall, nvar - 1, mband