dqed.f90

求解非线性方程组的一个高效算法

subroutine difcen ( fj, func, fx, iopt, ldfj, mcon, mequa, nvars, &
  ropt, x )
!
!***********************************************************************
!
!! DIFCEN estimates a jacobian using central differences.
!
!
!  Modified:
!
!    18 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, double precision FJ(LDFJ,NVARS), the estimated jacobian.
!
!    Input, external FUNC, the name of the user written
!    function evaluation routine.  FUNC should have the form:
!
!      subroutine func(fx,iopt,mcon,mequa,nvars,ropt,x)
!
!    and should accept X as input, and return in FX the value
!    of the MEQUA+MCON functions.
!
!    Workspace, double precision FX(MEQUA+MCON).
!
!    Throughput, integer IOPT(*), parameters to be passed to FUNC.
!
!    Input, integer LDFJ, the leading dimension of FJ, which must
!    be at least MEQUA+MCON.
!
!    Input, integer MCON, the number of constraints.
!
!    Input, integer MEQUA, the number of nonlinear functions.
!
!    Input, integer NVARS, the number of variables.
!
!    Throughput, double precision ROPT(*), parameters to be passed to FUNC.
!
!    Input, double precision X(NVARS), the point at which the
!    jacobian should be evaluated.
!
  implicit none
!
  integer ldfj
  integer mcon
  integer mequa
  integer nvars
!
  double precision dxj
  double precision eps
  double precision fj(ldfj,nvars)
  external func
  double precision fx(mequa+mcon)
  integer i
  integer iopt(*)
  integer j
  double precision ropt(*)
  double precision x(nvars)
  double precision xsave
!
!  Get the square root of the machine precision.
!
  eps = sqrt ( epsilon ( eps ) )
!
!  Consider each component X(J) of the set of variables.
!
  do j = 1, nvars
!
!  Set the appropriate increment DXJ to X(J).
!
    dxj = eps * ( abs ( x(j) ) + 1.0D+00 )
!
!  Make a copy XP of X, with X(J) incremented by DXJ.
!
    xsave = x(j)

    x(j) = xsave + dxj
!
!  Evaluate F(XP).
!
    call func ( fx, iopt, mcon, mequa, nvars, ropt, x )
!
!  Save F(XP).
!
    fj(1:mequa+mcon,j) = fx(1:mequa+mcon)
!
!  Make a copy XM of X, with X(J) decremented by DXJ.
!
    x(j) = xsave - dxj
!
!  Evaluate F(XM).
!
    call func ( fx, iopt, mcon, mequa, nvars, ropt, x )
!
!  Estimate the partial derivative d F/d X(J) by (F(XP)-F(XM))/(2*DXJ)
!
    fj(1:mequa+mcon,j) = ( fj(1:mequa+mcon,j) - fx(1:mequa+mcon) ) &
      / ( 2.0D+00 * dxj )
!
!  Restore the value of X(J).
!
    x(j) = xsave

  end do

  return
end
subroutine diffor ( fj, func, fx, iopt, ldfj, mcon, mequa, nvars, &
  ropt, x )
!
!***********************************************************************
!
!! DIFFOR estimates a jacobian using forward differences.
!
!
!  Modified:
!
!    18 February 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, double precision FJ(LDFJ,NVARS), the estimated jacobian.
!
!    Input, external FUNC, the name of the user written
!    function evaluation routine.  FUNC should have the form:
!
!      subroutine func(fx,iopt,mcon,mequa,nvars,ropt,x)
!
!    and should accept X as input, and return in FX the value
!    of the MEQUA+MCON functions.
!
!    Workspace, double precision FX(MEQUA+MCON).
!
!    Throughput, integer IOPT(*), parameters to be passed
!    to FUNC.
!
!    Input, integer LDFJ, the leading dimension of FJ, which must
!    be at least MEQUA+MCON.
!
!    Input, integer MCON, the number of constraints.
!
!    Input, integer MEQUA, the number of nonlinear functions.
!
!    Input, integer NVARS, the number of variables.
!
!    Throughput, double precision ROPT(*), parameters to be passed to FUNC.
!
!    Input, double precision X(NVARS), the point at which the
!    jacobian should be evaluated.
!
  implicit none
!
  integer ldfj
  integer mcon
  integer mequa
  integer nvars
!
  double precision dxj
  double precision eps
  double precision fj(ldfj,nvars)
  external func
  double precision fx(mequa+mcon)
  integer i
  integer iopt(*)
  integer j
  double precision ropt(*)
  double precision x(nvars)
  double precision xsave
!
!  Evaluate F(X) and save it in FX.
!
  call func ( fx, iopt, mcon, mequa, nvars, ropt, x )
!
!  Get the square root of the machine precision.
!
  eps = sqrt ( epsilon ( eps ) )
!
!  Consider each component X(J) of the set of variables.
!
  do j = 1, nvars
!
!  Set the appropriate increment DXJ to X(J).
!
    dxj = eps * ( abs ( x(j) ) + 1.0D+00 )
!
!  Make a copy XP of X, with X(J) incremented by DXJ.
!
    xsave = x(j)

    x(j) = xsave + dxj
!
!  Evaluate F(XP) and store it in column J.
!
    call func ( fj(1,j), iopt, mcon, mequa, nvars, ropt, x )
!
!  Estimate the partial derivative d F/d X(J) by (F(XP)-F(X))/DXJ
!
    fj(1:mequa+mcon,j) = ( fj(1:mequa+mcon,j) - fx(1:mequa+mcon) ) / dxj
!
!  Restore the value of X(J).
!
    x(j) = xsave

  end do

  return
end
function idamax ( n, x, incx )
!
!*******************************************************************************
!
!! IDAMAX finds the index of the vector element of maximum absolute value.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539,
!    ACM Transactions on Mathematical Software,
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, double precision X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of SX.
!
!    Output, integer IDAMAX, the index of the element of SX of maximum
!    absolute value.
!
  implicit none
!
  integer i
  integer incx
  integer idamax
  integer ix
  integer n
  double precision damax
  double precision x(*)
!
  if ( n <= 0 ) then

    idamax = 0

  else if ( n == 1 ) then

    idamax = 1

  else if ( incx == 1 ) then

    idamax = 1
    damax = abs ( x(1) )

    do i = 2, n

      if ( abs ( x(i) ) > damax ) then
        idamax = i
        damax = abs ( x(i) )
      end if

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    idamax = 1
    damax = abs ( x(ix) )

    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > damax ) then
        idamax = i
        damax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
subroutine ivout ( n, ix, title, idigit )
!
!***********************************************************************
!
!! IVOUT prints integer vectors.
!
!
!  Author:
!
!    John Wisniewski and Richard Hanson,
!    SANDIA LABS ALBUQUERQUE.
!
!  Parameters:
!
!    Input, integer N, the size of the vector IX.
!
!    Input, integer IX(N), the array to be printed.
!
!    Input, character ( len = * ) TITLE, a title to be printed.
!
!  IDIGIT  PRINT UP TO IABS(IDIGIT) DECIMAL DIGITS PER NUMBER.
!          THE SUBPROGRAM WILL CHOOSE THAT integer 4,6,10 OR 14
!          WHICH WILL PRINT AT LEAST IABS(IDIGIT) NUMBER OF
!          PLACES.  IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE UTILIZED
!          TO WRITE EACH LINE OF OUTPUT OF THE ARRAY IX(*). (THIS
!          CAN BE USED ON MOST TIME-SHARING TERMINALS). IF
!          IDIGIT.GE.0, 133 PRINTING COLUMNS ARE UTILIZED. (THIS CAN
!          BE USED ON MOST LINE PRINTERS).
!
  implicit none
!
  integer n
!
  integer i
  integer idigit
  integer ix(n)
  integer j
  integer k1
  integer k2
  integer ndigit
  character ( len = * ) title
!
  write ( *, '(a)' ) trim ( title )

  if ( n <= 0 ) then
    return
  end if

  ndigit = idigit

  if ( idigit == 0 ) then
    ndigit = 4
  end if

  if ( idigit >= 0 ) go to 80

  ndigit = -idigit

  if ( ndigit <= 4 )then

    do k1 = 1, n, 10
      k2 = min ( n, k1+9 )
      write(*,1000) k1, k2, ix(k1:k2)
    end do

    return

  end if

  if ( ndigit > 6) go to 40

  do k1=1,n,7
    k2 = min(n,k1+6)
    write(*,1001) k1,k2, ix(k1:k2)
  end do

  return

   40 continue
  if ( ndigit > 10) go to 60

  do k1=1,n,5
    k2=min(n,k1+4)
    write(*,1002) k1,k2, ix(k1:k2)
  end do

  return

   60 continue

  do k1=1,n,3
    k2 = min(n,k1+2)
    write(*,1003) k1,k2, ix(k1:k2)
  end do

  return

   80 continue

  if ( ndigit > 4) go to 100

  do k1=1,n,20
    k2 = min(n,k1+19)
    write(*,1000) k1,k2, ix(k1:k2)
  end do

  return

  100 continue

  if ( ndigit > 6) go to 120

  do k1=1,n,15
    k2 = min(n,k1+14)
    write(*,1001) k1,k2, ix(k1:k2)
  end do

  return

  120 continue

  if ( ndigit > 10) go to 140

  do k1=1,n,10
    k2 = min(n,k1+9)
    write(*,1002) k1,k2, ix(k1:k2)
  end do

  return

  140 continue

  do k1=1,n,7
    k2 = min(n,k1+6)
    write(*,1003) k1,k2, ix(k1:k2)
  end do

  return
 1000 format(1x,i4,' - ',i4,20(1x,i5))
 1001 format(1x,i4,' - ',i4,15(1x,i7))
 1002 format(1x,i4,' - ',i4,10(1x,i11))
 1003 format(1x,i4,' - ',i4,7(1x,i15))
end
function damax ( n, x, incx )
!
!*******************************************************************************
!
!! DAMAX returns the maximum absolute value of the entries in a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters: