dqed_prb.f90

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

!
!  MODE = 1, there are constraints.
!  (The first two equations are linear, and can be used as
!  constraints instead).
!
  else

    mcon = 2
    ind(nvars+1) = 3
    ind(nvars+2) = 3
    bl(nvars+1) = sigma(1)
    bu(nvars+1) = sigma(1)
    bl(nvars+2) = sigma(2)
    bu(nvars+2) = sigma(2)

  end if
 
  mequa = nvars - mcon
!
!  All variables are otherwise free.
!
  ind(1:nvars) = 4

  call dqed ( fjaprx, mequa, nvars, mcon, ind, bl, bu, x, fj, ldfj, fnorm, &
    igo, iopt, ropt, iwa, wa )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Computed minimizing X:'
  write ( *, '(a)' ) ' '
  do i = 1, nvars
    write ( *, '(g14.6)' ) x(i)
  end do
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) 'Residual after the fit = ',fnorm
  write ( *, '(a,i6)' ) 'QED output flag IGO = ',igo
 
  return
end
subroutine dqedhd ( x, fj, ldfj, igo, iopt, ropt )
!
!***********************************************************************
!
!! DQEDHD evaluates functions and derivatives for DQED.
!
!
!  Discussion:
!
!    The user problem has MCON constraint functions,
!    MEQUA least squares equations, and involves NVARS
!    unknown variables.
!
!    When this subprogram is entered, the general (near)
!    linear constraint partial derivatives, the derivatives
!    for the least squares equations, and the associated
!    function values are placed into the array FJ.
!    all partials and functions are evaluated at the point
!    in X.  then the subprogram returns to the calling
!    program unit. typically one could do the following
!    steps:
!
!    if ( igo /= 0 ) then
!      place the partials of the i-th constraint function with respect to 
!      variable j in the array fj(i,j), i = 1,...,mcon, j=1,...,nvars.
!
!    place the values of the i-th constraint equation into fj(i,nvars+1).
!
!    if ( igo /= 0 ) then
!      place the partials of the i-th least squares equation with respect 
!      to variable j in the array fj(i,j), i = 1,...,mequa, j = 1,...,nvars.
!
!    place the value of the i-th least squares equation into fj(i,nvars+1).
!
  implicit none
!
  integer ldfj
  integer, parameter :: nvars = 8
!
  double precision fj(ldfj,nvars+1)
  integer igo
  integer iopt(*)
  integer mcon
  integer mequa
  integer mode
  double precision ropt(*)
  double precision x(nvars)
!
  common /mode/  mode
!
  if ( mode == 0 ) then
    mcon = 0
  else
    mcon = 2
  end if

  mequa = nvars - mcon

  call jack ( fj, ldfj, nvars, x )

  call func ( fj(1,9), iopt, mcon, mequa, nvars, ropt, x )

  return
end
subroutine fjaprx ( x, fj, ldfj, igo, iopt, ropt )
!
!***********************************************************************
!
!! FJAPRX uses DIFCEN to approximate the jacobian matrix.
!
  implicit none
!
  integer ldfj
  integer, parameter :: nvars = 8
!
  double precision fj(ldfj,nvars+1)
  double precision fx(8)
  integer igo
  integer iopt(*)
  integer mcon
  integer mequa
  integer mode
  double precision ropt(*)
  double precision x(nvars)
!
  external func
!
  common /mode/  mode
!
  if ( mode == 0 ) then
    mcon = 0
  else
    mcon = 2
  end if

  mequa = nvars - mcon

  call difcen ( fj, func, fx, iopt, ldfj, mcon, mequa, nvars, ropt, x )

  call func ( fj(1,nvars+1), iopt, mcon, mequa, nvars, ropt, x )

  return
end
subroutine func ( fx, iopt, mcon, mequa, nvars, ropt, x )
!
!***********************************************************************
!
!! FUNC evaluates the functions.
!
!    x(1)+x(2) = sigma(1)
!
!    x(3)+x(4) = sigma(2)
!
!    x(1)*x(5)+x(2)*x(6)-x(3)*x(7)-x(4)*x(8) = sigma(3)
!
!    x(1)*x(7)+x(2)*x(8)+x(3)*x(5)+x(4)*x(6) = sigma(4)
!
!    x(1)*(x(5)**2-x(7)**2)+x(2)*(x(6)**2-x(8)**2)+x(3)*(-2.0*x(5)*x(7))
!      +x(4)*(-2.0*x(6)*x(8)) = sigma(5)
!
!    x(1)*(2.0*x(5)*x(7))+x(2)*(2.0*x(6)*x(8))+x(3)*(x(5)**2-x(7)**2)
!      +x(4)*(x(6)**2-x(8)**2) = sigma(6)
!
!    x(1)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(2)*(x(6)*(x(6)**2-3.0*x(8)**2))
!      +x(3)*(x(7)*(x(7)**2-3.0*x(5)**2))+x(4)*(x(8)*(x(8)**2-3.0*x(6)**2))
!       = sigma(7)
!
!    x(1)*(-x(7)*(x(7)**2-3.0*x(5)**2))+x(2)*(-x(8)*(x(8)**2-3.0*x(6)**2))
!      +x(3)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(4)*(x(6)*(x(6)**2-3.0*x(8)**2))
!       = sigma(8)
!
  implicit none
!
  integer mcon
  integer mequa
  integer nvars
!
  double precision fx(mequa+mcon)
  integer iopt(*)
  integer mode
  double precision ropt(*)
  double precision sigma(8)
  double precision x(nvars)
!
  common /sigma/ sigma
  common /mode/ mode
!
  if ( mode == 0 ) then
    fx(1) = x(1) + x(2) - sigma(1)
  else
    fx(1) = x(1) + x(2)
  end if

  if ( mode == 0 ) then
    fx(2) = x(3) + x(4) - sigma(2)
  else
    fx(2) = x(3) + x(4)
  end if

  fx(3) = x(1)*x(5) + x(2)*x(6) - x(3)*x(7) - x(4)*x(8) - sigma(3)

  fx(4) = x(1)*x(7) + x(2)*x(8) + x(3)*x(5) + x(4)*x(6) - sigma(4)

  fx(5) = x(1)*(x(5)**2-x(7)**2) + x(2)*(x(6)**2-x(8)**2) &
    +x(3)*(-2.0D+00*x(5)*x(7)) + x(4)*(-2.0D+00*x(6)*x(8)) - sigma(5)

  fx(6) = x(1)*(2.0D+00*x(5)*x(7)) + x(2)*(2.0D+00*x(6)*x(8)) &
    + x(3)*(x(5)**2-x(7)**2) &
    + x(4)*(x(6)**2-x(8)**2) - sigma(6)

  fx(7) = x(1)*(x(5)*(x(5)**2-3.0*x(7)**2)) &
    + x(2)*(x(6)*(x(6)**2-3.0D+00*x(8)**2)) &
    + x(3)*(x(7)*(x(7)**2-3.0D+00*x(5)**2)) &
    + x(4)*(x(8)*(x(8)**2-3.0D+00*x(6)**2)) &
    - sigma(7)

  fx(8) = x(1)*(-x(7)*(x(7)**2-3.0D+00*x(5)**2)) &
    + x(2)*(-x(8)*(x(8)**2-3.0D+00*x(6)**2)) &
    + x(3)*(x(5)*(x(5)**2-3.0D+00*x(7)**2)) &
    + x(4)*(x(6)*(x(6)**2-3.0D+00*x(8)**2)) - sigma(8)

  return
end
subroutine jack ( fj, ldfj, nvars, x )
!
!***********************************************************************
!
!! JACK evaluates the partial derivatives of the functions.
!
  implicit none
!
  integer ldfj
  integer nvars
!
  double precision fj(ldfj,nvars)
  double precision x(nvars)
!
!  Equation #1: 
!
!    x(1)+x(2) = sigma(1)
!
  fj(1,1) = 1.0D+00
  fj(1,2) = 1.0D+00
  fj(1,3) = 0.0D+00
  fj(1,4) = 0.0D+00
  fj(1,5) = 0.0D+00
  fj(1,6) = 0.0D+00
  fj(1,7) = 0.0D+00
  fj(1,8) = 0.0D+00
!
!  Equation #2:
!
!    x(3)+x(4) = sigma(2)
!
  fj(2,1) = 0.0D+00
  fj(2,2) = 0.0D+00
  fj(2,3) = 1.0D+00
  fj(2,4) = 1.0D+00
  fj(2,5) = 0.0D+00
  fj(2,6) = 0.0D+00
  fj(2,7) = 0.0D+00
  fj(2,8) = 0.0D+00
!
!  Equation #3:
!    x(1)*x(5)+x(2)*x(6)-x(3)*x(7)-x(4)*x(8) = sigma(3)
!
  fj(3,1) = x(5)
  fj(3,2) = x(6)
  fj(3,3) = -x(7)
  fj(3,4) = -x(8)
  fj(3,5) = x(1)
  fj(3,6) = x(2)
  fj(3,7) = -x(3)
  fj(3,8) = -x(4)
!
!  Equation #4:
!    x(1)*x(7)+x(2)*x(8)+x(3)*x(5)+x(4)*x(6) = sigma(4)
!
  fj(4,1) = x(7)
  fj(4,2) = x(8)
  fj(4,3) = x(5)
  fj(4,4) = x(6)
  fj(4,5) = x(3)
  fj(4,6) = x(4)
  fj(4,7) = x(1)
  fj(4,8) = x(2)
!
!  Equation #5:
!    x(1)*(x(5)**2-x(7)**2)+x(2)*(x(6)**2-x(8)**2)+x(3)*(-2.0*x(5)*x(7))
!      +x(4)*(-2.0*x(6)*x(8)) = sigma(5)
!
  fj(5,1) = x(5)**2-x(7)**2
  fj(5,2) = x(6)**2-x(8)**2
  fj(5,3) = -2.0D+00 *x(5)*x(7)
  fj(5,4) = -2.0D+00 *x(6)*x(8)
  fj(5,5) = 2.0D+00 *(x(1)*x(5)-x(3)*x(7))
  fj(5,6) = 2.0D+00 *(x(2)*x(6)-x(4)*x(8))
  fj(5,7) = -2.0D+00 *(x(1)*x(7)+x(3)*x(5))
  fj(5,8) = -2.0D+00 *(x(2)*x(8)+x(4)*x(6))
!
!  Equation #6:
!    x(1)*(2.0*x(5)*x(7))+x(2)*(2.0*x(6)*x(8))+x(3)*(x(5)**2-x(7)**2)
!      +x(4)*(x(6)**2-x(8)**2) = sigma(6)
!
  fj(6,1) = 2.0D+00 *x(5)*x(7)
  fj(6,2) = 2.0D+00 *x(6)*x(8)
  fj(6,3) = x(5)**2-x(7)**2
  fj(6,4) = x(6)**2-x(8)**2
  fj(6,5) = 2.0D+00 *(x(1)*x(7)+x(3)*x(5))
  fj(6,6) = 2.0D+00 *(x(2)*x(8)+x(4)*x(6))
  fj(6,7) = 2.0D+00 *(x(1)*x(5)-x(3)*x(7))
  fj(6,8) = 2.0D+00 *(x(2)*x(6)-x(4)*x(8))
!
!  Equation #7:
!    x(1)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(2)*(x(6)*(x(6)**2-3.0*x(8)**2))
!      +x(3)*(x(7)*(x(7)**2-3.0*x(5)**2))+x(4)*(x(8)*(x(8)**2-3.0*x(6)**2))
!       = sigma(7)
!
  fj(7,1) = x(5)*(x(5)**2-3.0D+00*x(7)**2)
  fj(7,2) = x(6)*(x(6)**2-3.0D+00*x(8)**2)
  fj(7,3) = x(7)*(x(7)**2-3.0D+00*x(5)**2)
  fj(7,4) = x(8)*(x(8)**2-3.0D+00*x(6)**2)
  fj(7,5) = 3.0D+00*(x(1)*(x(5)**2-x(7)**2)+x(3)*(-2.0D+00*x(5)*x(7)))
  fj(7,6) = 3.0D+00*(x(2)*(x(6)**2-x(8)**2)+x(4)*(-2.0D+00*x(6)*x(8)))
  fj(7,7) = -3.0D+00*(x(1)*(2.0*x(5)*x(7))+x(3)*(x(5)**2-x(7)**2))
  fj(7,8) = -3.0D+00*(x(2)*(2.0*x(6)*x(8))+x(4)*(x(6)**2-x(8)**2))
!
!  Equation #8:
!    x(1)*(-x(7)*(x(7)**2-3.0*x(5)**2))+x(2)*(-x(8)*(x(8)**2-3.0*x(6)**2))
!      +x(3)*(x(5)*(x(5)**2-3.0*x(7)**2))+x(4)*(x(6)*(x(6)**2-3.0*x(8)**2))
!       = sigma(8)
!
  fj(8,1) = -x(7)*(x(7)**2-3.0D+00*x(5)**2)
  fj(8,2) = -x(8)*(x(8)**2-3.0D+00*x(6)**2)
  fj(8,3) = x(5)*(x(5)**2-3.0D+00*x(7)**2)
  fj(8,4) = x(6)*(x(6)**2-3.0D+00*x(8)**2)
  fj(8,5) = 3.0D+00*(x(1)*(2.0D+00*x(5)*x(7))+x(3)*(x(5)**2-x(7)**2))
  fj(8,6) = 3.0D+00*(x(2)*(2.0D+00*x(6)*x(8))+x(4)*(x(6)**2-x(8)**2))
  fj(8,7) = 3.0D+00*(x(1)*(x(5)**2-x(7)**2)+x(3)*(-2.0D+00*x(5)*x(7)))
  fj(8,8) = 3.0D+00*(x(2)*(x(6)**2-x(8)**2)+x(4)*(-2.0D+00*x(6)*x(8)))

  return
end