pitcon.f90

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

!  This is a two point boundary value problem, representing the behavior
!  of a rod under a load.  The shape of the rod is represented by piecewise
!  polynomials, whose form and continuity are specifiable as options.
!
!  The problem as programmed can allow up to 71 variables, but a simple
!  change of a few parameters will allow the problem to be arbitrarily
!  large.
!
!  This is a complicated program.  It is intended to demonstrate the
!  advanced kinds of problems that can be solved with PITCON, as opposed
!  to the "toy" problems like the Freudenstein-Roth function.
!
!
!  8:
!  pcprb8.f
!  The Freudenstein-Roth function.  3 variables.
!
!  This version of the problem tests the jacobian approximation options.
!  We run the problem 3 times, first with the user jacobian, second with
!  the forward difference approximation, and third with the central
!  difference approximation.
!
!
!  Recent Changes:
!  --------------
!
!
!  Changes made for version 7.0:
!
!    Inserted simple versions of LAPACK routines.
!
!  Changes made for version 6.6:
!
!    Some calculations in COQUAL, for ITOP and IBOT, could involve integers
!    to negative powers.  These were replaced by real calculations for TOP
!    and BOT.
!
!    There is a stringent convergence test in CORRECTOR which severely slows
!    down the traversal of the Freudenstein Roth curve, because it forces
!    a last very small step, which causes the computation of the stepsize
!    to be skewed.  I temporarily turned this off.
!
!  Changes made for version 6.5:
!
!    Spun off "STAJAC" from START.
!
!    Code written in lower case.
!
!    Replaced all labeled DO statements with DO/ENDDO loops.
!
!  Changes made for version 6.4:
!
!    Calls to LINPACK replaced by calls to LAPACK
!    Added routines SGEDET and SGBDET to compute determinants.
!
!    Call to SGBTRF, SGBTRS, SGBDET had incorrect value of LDA,
!    which was corrected 21 January 1994.
!
!    Dropped LOUNIT.  All WRITE's are to unit * now.
!
!
!  Changes made for version 6.3:
!
!    27 July 1992: Printout of "Possible bifurcation" message will
!    now occur if IWRITE.GE.1.
!
!    HMIN was reset to MAX(HMIN,SQRT(EPS)) before.  Now it is only
!    set to SQRT(EPS) if HMIN is nonpositive.
!
!    30 September 1991: SKALE was not declared in LIMIT, causing
!    problems for the double precision code.
!
!    26 September 1991: corrected the formulation of the limit
!    point calculation.
!
!    12 September 1991: Jyun-Ming Chen reported that the program did not
!    actually fix the parameter to a user specified value when requested to do
!    so.  That is, when IWORK(3) = 1, the program is never supposed to alter the
!    input value of IWORK(2).  However, the program was doing so.  This error
!    has been corrected now.  The only exception occurs when the initial
!    input value of IWORK(2) is less than 1 or greater than NVAR, in which
!    case the program sets IWORK(2) to NVAR, no matter what the value
!    of IWORK(3).
!
!    The target and limit point calculations, and many other operations,
!    were "spun off" into subroutines.
!
!
!  Changes made for version 6.2:
!
!    The internal documentation was corrected.  It originally stated that
!    the user input portions of IWORK were entries 1 through 8, 17 and 29.
!    This was corrected to 1 through 9 and 17.
!
!    The entry RWORK(18) was previously unused.  It is now used to allow
!    the user to set the value of the finite difference increment used
!    for the approximation of the jacobian.
!
!    Modified DENSLV and BANSLV to allow more user control over the
!    Jacobian comparison or printout.  IWORK(1) may now have the values -1
!    through -10.  The user thus chooses to use forward or central
!    differences, and minimal or maximal printout.
!
!    Modified "Programming notes" section to explain the role of the
!    REPS routine, the new use of IWORK(1), and the new parameter RWORK(18).
!
!
!  Changes made for version 6.0:
!
!    1) When computing a starting point, the Newton convergence
!    criterion was relaxed to require only that the function
!    norm decrease, but not to require that the step norm
!    decrease.
!
!    2) In the Newton correction routine, the MAX norm was replaced
!    by the L2 norm when computing the size of the step and the
!    residual.  This was an attempt to control the 'jerkiness'
!    of the convergence for poorly scaled problems.
!
!    3) Added option to check Jacobian.
!
!    4) Apparently, a programming error in the previous version meant
!    that IWORK(29) was set to zero on startup.  This meant that the
!    program would not realize that the user was using a band storage
!    scheme.  This has been corrected.
!
  external df
  external fx
  external slname
!
  double precision, parameter :: one = 1.0D+00
!
  integer liw
  integer lrw
  integer nvar
!
  double precision dets
  double precision fpar(*)
  integer i
  integer ierror
  integer ifound
  integer ipar(*)
  integer iwork(liw)
  integer iwrite
  integer job
  integer ltc
  integer lwk
  integer lxc
  integer lxf
  integer modnew
  double precision qual
  double precision rwork(lrw)
  double precision xr(nvar)
!
!******************************************************************************
!
!  1.  Set a few parameters.
!
!******************************************************************************
!
  ierror = 0
  iwrite = iwork(7)

  lxc = 29            + 1
  lxf = 29 +     nvar + 1
  ltc = 29 + 2 * nvar + 1
  lwk = 29 + 3 * nvar + 1

  if ( iwork(1) < -10 ) then
    iwork(1) = - 2
  else if ( iwork(1) > 5 ) then
    iwork(1) = 0
  end if
!
!******************************************************************************
!
!  2.  Preparations.
!
!  Check entries of IWORK and RWORK, compute some constants.
!
!******************************************************************************
!
  if ( iwork(1) == 0 ) then

    if ( iwrite >= 1 ) then
      write ( *, * ) ' '
      write ( *, * ) 'PITCON 7.0'
      write ( *, * ) '  University of Pittsburgh continuation code'
      write ( *, * ) ' '
      write ( *, * ) '  Modified:        12 November 1999'
      write ( *, * ) '  Linear algebra:  LAPACK'
      write ( *, * ) '  Precision:       double'
      write ( *, * ) ' '
    end if

    call checkw ( ierror, iwork, liw, lrw, nvar, rwork )

    if ( ierror /= 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'PITCON - Fatal error!'
      write ( *, * ) '  An error was detected in the user input.'
      write ( *, * ) '  PITCON can not begin.'
      stop
    end if

  end if
!
!*******************************************************************************
!
!  3.  Check the user Jacobian routine against a finite difference 
!  approximation.
!
!*******************************************************************************
!
  if ( iwork(1) < 0 ) then

    job = 3

    call slname ( dets, fx, df, fpar, ierror, iwork(2), ipar, iwork, &
      liw, job, nvar, rwork, lrw, xr, rwork(lwk) )

    if ( ierror /= 0 ) then
      if ( iwrite > 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'PITCON - Warning!'
        write ( *, * ) '  An error occurred during the jacobian '
        write ( *, * ) '  check.'
      end if
    end if

    return

  end if
!
!*******************************************************************************
!
!  4.  Starting point check
!
!  On first call for a given problem, check that F(XR) is small
!  enough so that the starting point may be considered to lie on
!  the curve.
!
!  If this is not the case, call the corrector to try to enforce it.
!
!*******************************************************************************
!
  if ( iwork(1) == 0 ) then
!
!  See if the initial jacobian should be set up.
!
    if ( iwork(4) == 2 ) then

      call stajac ( df,fpar,fx,ierror,ipar,iwork(2),iwork,liw, &
        lrw,nvar,rwork,rwork(lwk),xr,slname)

      if ( ierror /= 0 ) then
        if ( iwrite > 0 ) then
          write ( *, * ) ' '
          write ( *, * ) 'PITCON - Fatal error!'
          write ( *, * ) '  An error occurred during the'
          write ( *, * ) '  initial jacobian setup.'
          write ( *, * ) ' '
          write ( *, * ) '  The program can not continue!'
          stop
        end if
        return
      end if

    end if
!
!  Check that the starting point satisfies the equations.
!
    call start ( df,fpar,fx,ierror,ipar,iwork(2),iwork,liw, &
      lrw,nvar,rwork,rwork(lwk),rwork(lxc),rwork(lxf),xr,slname)

    if ( ierror /= 0 ) then
      if ( iwrite > 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'PITCON - Fatal error!'
        write ( *, * ) '  An error occurred during the starting'
        write ( *, * ) '  point check.'
        write ( *, * ) ' '
        write ( *, * ) '  The starting point does not satisfy the'
        write ( *, * ) '  accuracy requirements, and PITCON could'
        write ( *, * ) '  not correct it.'
        write ( *, * ) ' '
        write ( *, * ) '  The program can not continue!'
      end if
    end if

    do i = 1, nvar
      rwork(ltc+i-1) = 0.0D+00
    end do

    rwork(ltc+iwork(2)-1) = 1.0D+00
    return

  end if
!
!*******************************************************************************
!
!  5.  Target point
!
!  If IWORK(5) is nonzero, target points are sought.  Check to see if
!  target component IWORK(5), also called "IT", has value lying between
!  XC(IT) and XF(IT).  If so, get linearly interpolated starting point,
!  and use Newton's method to get target point.
!
!*******************************************************************************
!
  call target ( df, fpar, fx, ierror, ifound, ipar, iwork, liw, lrw, &
    nvar, rwork, slname, rwork(lwk), rwork(lxc), rwork(lxf), xr )

  if ( ifound == 1 ) then
    if ( ierror /= 0 ) then
      ierror = 9
      if ( iwrite > 0 ) then
        write ( *, * ) ' '
        write ( *, * ) 'PITCON - Warning!'
        write ( *, * ) '  An error occurred during the'