dqed_prb.f90

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

program dqed_prb
!
!***********************************************************************
!
!! DQED_PRB tests DQED.
!
!  This routine illustrates the use of DQED, the Hanson-Krogh nonlinear least
!  squares solver, by solving a certain heart dipole moment equation.
!
!  Reference: 
!
!    John Dennis, David Gay, and Phuong Vu.
!    A new nonlinear equations test problem, 
!    Rice University Department of Mathematics 
!    Report 83-16, 6/83, 6/85
!
  implicit none
!
  integer, parameter :: ldfj = 8
  integer, parameter :: liwa = 150
  integer, parameter :: lwa = 2000
  integer, parameter :: nvars = 8
!
  double precision bl(nvars+2)
  double precision bu(nvars+2)
  double precision df(nvars)
  double precision fj(ldfj,nvars+1)
  integer i
  integer ind(nvars+2)
  integer iopt(28)
  integer ios
  integer iwa(liwa)
  integer mode
  double precision ropt(1)
  double precision sigma(nvars)
  character ( len = 20 ) title
  double precision wa(lwa)
  double precision x(nvars)
  double precision xsave(nvars)
  double precision y(nvars)
!
  common /sigma/ sigma
  common /mode/  mode
!
  external dqedhd
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DQED_PRB'
  write ( *, '(a)' ) '  A set of tests for DQED, which can solve'
  write ( *, '(a)' ) '  bounded and constrained linear least squares problems'
  write ( *, '(a)' ) '  and systems of nonlinear equations.'

  do
!
!  Read the problem title.
!
    read ( *, '(a)', iostat = ios ) title

    if ( ios /= 0 ) then
      exit
    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) trim ( title )
!
!  Read in the sums of values and initial parameter values.
!
    read ( *, *, iostat = ios ) sigma(1:nvars), xsave(1:nvars)

    if ( ios /= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Unexpected end of file or other error.'
      exit
    end if

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Input SIGMA:'
    write ( *, '(a)' ) ' '
    do i = 1, nvars
      write ( *, '(i6,g14.6)' ) i, sigma(i) 
    end do

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Input X:'
    write ( *, '(a)' ) ' '
    do i = 1, nvars
      write ( *, '(i6,g14.6)' ) i, xsave(i) 
    end do
!
!  Test the partial derivative computation.
!
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'Test the partial derivative computation:'
    write ( *, '(a)' ) ' '

    call dpchek ( df, dqedhd, fj, iopt, ldfj, nvars, ropt, xsave, y )
!
!  Test 1, with all equations normal.
!
    mode = 0

    x(1:nvars) = xsave(1:nvars)

    call test01 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, &
      ropt, sigma, wa, x )
!
!  Test 1, with first two linear equations used as constraints.
!
    mode = 1

    x(1:nvars) = xsave(1:nvars)

    call test01 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, &
      ropt, sigma, wa, x )
!
!  Test 2, with all equations normal.
!
    mode = 0

    x(1:nvars) = xsave(1:nvars)

    call test02 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, &
      ropt, sigma, wa, x )
!
!  Test 2, with first two linear equations used as constraints.
!
    mode = 1

    x(1:nvars) = xsave(1:nvars)

    call test02 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, nvars, &
      ropt, sigma, wa, x )

  end do
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DQED_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, &
  nvars, ropt, sigma, wa, x )
!
!***********************************************************************
!
!! TEST01 uses an analytic jacobian.
!
!  The linear equations will not be constraints if MODE = 0.
!  Convergence is normally faster if the linear equations are used
!  as constraints, MODE = 1.
!
  implicit none
!
  integer ldfj
  integer liwa
  integer lwa
  integer nvars
!
  double precision bl(nvars+2)
  double precision bu(nvars+2)
  external dqedhd
  double precision fj(ldfj,nvars+1)
  double precision fnorm
  integer i
  integer igo
  integer ind(nvars+2)
  integer iopt(28)
  integer iwa(liwa)
  integer mcon
  integer mequa
  integer mode
  double precision ropt(1)
  double precision sigma(nvars)
  double precision wa(lwa)
  double precision x(nvars)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  Use an analytic jacobian.'
  write ( *, '(a,i6)' ) '  MODE = ', mode
!
!  Tell how much storage the solver has.
!
  iwa(1) = lwa
  iwa(2) = liwa
!
!  Set up print option.  (not used here).
!
  iopt(1) = -1
  iopt(2) = 1
!
!  Set up for linear model without any quadratic terms. (not used).
!
  iopt(3) = -14
  iopt(4) = 1
!
!  Do not allow convergence to be claimed on small steps.
!
  iopt(5) = 17
  iopt(6) = 1
  iopt(7) = 15
!
!  Allow up to NVARS quadratic model terms.
!
  iopt(8) = nvars
!
!  Change condition number for quadratic model degree.
!
  iopt(9) = 10
  iopt(10) = 1
  ropt(1) = 1.0D+04
!
!  No more options.
!
  iopt(11) = 99
!
!  MODE = 0, no constraints.
!
  if ( mode == 0 ) then

    mcon = 0
!
!  MODE = 1, there are constraints.
!
!  (The first two equations are linear, and can be used as constraints).
!
  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 ( dqedhd, 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,g14.6)' ) 'Residual after the fit = ',fnorm
  write ( *, '(a,i6)' ) 'QED output flag IGO = ',igo
 
  return
end
subroutine test02 ( bl, bu, fj, ind, iopt, iwa, ldfj, liwa, lwa, mode, &
  nvars, ropt, sigma, wa, x )
!
!***********************************************************************
!
!! TEST02 uses a finite-difference approximated jacobian.
!
!  The linear equations will not be constraints if MODE = 0.
!  Convergence is normally faster if the linear equations are used
!  as constraints, MODE = 1.
!
  implicit none
!
  integer ldfj
  integer liwa
  integer lwa
  integer nvars
!
  double precision bl(nvars+2)
  double precision bu(nvars+2)
  double precision fj(ldfj,nvars+1)
  external fjaprx
  double precision fnorm
  integer i
  integer igo
  integer ind(nvars+2)
  integer iopt(28)
  integer iwa(liwa)
  integer mcon
  integer mequa
  integer mode
  double precision ropt(1)
  double precision sigma(nvars)
  double precision wa(lwa)
  double precision x(nvars)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Use an approximate jacobian.'
  write ( *, '(a,i6)' ) '  MODE = ', mode
!
!  Tell how much storage the solver has.
!
  iwa(1) = lwa
  iwa(2) = liwa
!
!  Set up print option.  (not used here).
!
  iopt(1) = -1
  iopt(2) = 1
!
!  Set up for linear model without any quadratic terms. (not used).
!
  iopt(3) = -14
  iopt(4) = 1
!
!  Do not allow convergence to be claimed on small steps.
!
  iopt(5) = 17
  iopt(6) = 1
  iopt(7) = 15
!
!  Allow up to NVARS quadratic model terms.
!
  iopt(8) = nvars
!
!  Change condition number for quadratic model degree.
!
  iopt(9) = 10
  iopt(10) = 1
  ropt(1) = 10000.0D+00
!
!  No more options.
!
  iopt(11) = 99
!
!  MODE = 0, no constraints.
!
  if ( mode == 0 ) then

    mcon = 0