dmain.f90

牛顿优化算法源fortran代码

!  Brett M. Averick and Jorge J. More'.

!  **********

REAL (dp), PARAMETER :: p5=0.5_dp, three=3.0_dp

LOGICAL   :: feval, geval
INTEGER   :: i, j, k
REAL (dp) :: area, cdiv3, dvdx, dvdy, flin, fquad, hx, hy,  &
             temp, temp1, v, vb, vl, vr, vt

hx = one/DBLE(nx+1)
hy = one/DBLE(ny+1)
area = p5*hx*hy
cdiv3 = c/three

!     Compute a lower bound for x if task = 'XL' or
!     an upper bound if task = 'XU'.

IF (task(1:2) == 'XL' .OR. task(1:2) == 'XU') THEN
  IF (task(1:2) == 'XL') temp1 = -one
  IF (task(1:2) == 'XU') temp1 = one
  DO  j = 1, ny
    temp = DBLE(MIN(j,ny-j+1))*hy
    DO  i = 1, nx
      k = nx*(j-1) + i
      x(k) = SIGN(MIN(DBLE(MIN(i,nx-i+1))*hx,temp),temp1)
    END DO
  END DO
  
  RETURN
  
END IF

!     Compute the standard starting point if task = 'XS'.

IF (task(1:2) == 'XS') THEN
  DO  j = 1, ny
    temp = DBLE(MIN(j,ny-j+1))*hy
    DO  i = 1, nx
      k = nx*(j-1) + i
      x(k) = MIN(DBLE(MIN(i,nx-i+1))*hx,temp)
    END DO
  END DO
  
  RETURN
  
END IF

IF (task(1:1) == 'F' .OR. task(1:2) == 'FG') THEN
  feval = .true.
ELSE
  feval = .false.
END IF
IF (task(1:1) == 'G' .OR. task(1:2) == 'FG') THEN
  geval = .true.
ELSE
  geval = .false.
END IF

!     Evaluate the function if task = 'F', the gradient if task = 'G',
!     or both if task = 'FG'.

IF (feval) THEN
  fquad = zero
  flin = zero
END IF
IF (geval) THEN
  DO  k = 1, nx*ny
    fgrad(k) = zero
  END DO
END IF

!     Computation of the function and the gradient over the lower
!     triangular elements.

DO  j = 0, ny
  DO  i = 0, nx
    k = nx*(j-1) + i
    v = zero
    vr = zero
    vt = zero
    IF (i >= 1 .AND. j >= 1) v = x(k)
    IF (i < nx .AND. j > 0) vr = x(k+1)
    IF (i > 0 .AND. j < ny) vt = x(k+nx)
    dvdx = (vr-v)/hx
    dvdy = (vt-v)/hy
    IF (feval) THEN
      fquad = fquad + dvdx**2 + dvdy**2
      flin = flin - cdiv3*(v+vr+vt)
    END IF
    IF (geval) THEN
      IF (i /= 0 .AND. j /= 0)  &
          fgrad(k) = fgrad(k) - dvdx/hx - dvdy/hy - cdiv3
      IF (i /= nx .AND. j /= 0) fgrad(k+1) = fgrad(k+1) + dvdx/hx - cdiv3
      IF (i /= 0 .AND. j /= ny) fgrad(k+nx) = fgrad(k+nx) + dvdy/hy - cdiv3
    END IF
  END DO
END DO

!     Computation of the function and the gradient over the upper
!     triangular elements.

DO  j = 1, ny + 1
  DO  i = 1, nx + 1
    k = nx*(j-1) + i
    vb = zero
    vl = zero
    v = zero
    IF (i <= nx .AND. j > 1) vb = x(k-nx)
    IF (i > 1 .AND. j <= ny) vl = x(k-1)
    IF (i <= nx .AND. j <= ny) v = x(k)
    dvdx = (v-vl)/hx
    dvdy = (v-vb)/hy
    IF (feval) THEN
      fquad = fquad + dvdx**2 + dvdy**2
      flin = flin - cdiv3*(vb+vl+v)
    END IF
    IF (geval) THEN
      IF (i /= nx+1 .AND. j /= 1) fgrad(k-nx) = fgrad(k-nx) - dvdy/hy - cdiv3
      IF (i /= 1 .AND. j /= ny+1) fgrad(k-1) = fgrad(k-1) - dvdx/hx - cdiv3
      IF (i /= nx+1 .AND. j /= ny+1)  &
          fgrad(k) = fgrad(k) + dvdx/hx + dvdy/hy - cdiv3
    END IF
  END DO
END DO

!     Scale the result.

IF (feval) f = area*(p5*fquad+flin)
IF (geval) THEN
  DO  k = 1, nx*ny
    fgrad(k) = area*fgrad(k)
  END DO
END IF

RETURN
END SUBROUTINE deptfg



SUBROUTINE dpjbfg(nx, ny, x, f, fgrad, task, ecc, b)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:49:08

INTEGER, INTENT(IN)            :: nx
INTEGER, INTENT(IN)            :: ny
REAL (dp), INTENT(OUT)         :: x(:)    ! x(nx*ny)
REAL (dp), INTENT(OUT)         :: f
REAL (dp), INTENT(OUT)         :: fgrad(:)    ! fgrad(nx*ny)
CHARACTER (LEN=*), INTENT(IN)  :: task
REAL (dp), INTENT(IN)          :: ecc
REAL (dp), INTENT(IN)          :: b

!  **********

!  Subroutine dpjbfg

!  This subroutine computes the function and gradient of the
!  pressure distribution in a journal bearing problem.

!  The subroutine statement is

!    subroutine dpjbfg(nx, ny, x, f, fgrad, task, ecc, b)

!  where

!    nx is an integer variable.
!      On entry nx is the number of grid points in the first
!         coordinate direction.
!      On exit nx is unchanged.

!    ny is an integer variable.
!      On entry ny is the number of grid points in the second
!         coordinate direction.
!      On exit ny is unchanged.

!    x is a REAL (dp) array of dimension nx*ny.
!      On entry x specifies the vector x if task = 'F', 'G', or 'FG'.
!         Otherwise x need not be specified.
!      On exit x is unchanged if task = 'F', 'G', or 'FG'. Otherwise
!         x is set according to task.

!    f is a REAL (dp) variable.
!      On entry f need not be specified.
!      On exit f is set to the function evaluated at x if task = 'F' or 'FG'.

!    fgrad is a REAL (dp) array of dimension nx*ny.
!      On entry fgrad need not be specified.
!      On exit fgrad contains the gradient evaluated at x if task = 'G' or 'FG'.

!    task is a character*60 variable.
!      On entry task specifies the action of the subroutine:

!         task               action
!         ----               ------
!          'F'     Evaluate the function at x.
!          'G'     Evaluate the gradient at x.
!          'FG'    Evaluate the function and the gradient at x.
!          'XS'    Set x to the standard starting point xs.
!          'XL'    Set x to the lower bound xl.

!      On exit task is unchanged.

!    ecc is a REAL (dp) variable
!      On entry ecc is the eccentricity in (0,1).
!      On exit ecc is unchanged

!    b is a REAL (dp) variable
!      On entry b defines the domain as D = (0,2*pi) X (0,2*b).
!      On exit b is unchanged.

!  MINPACK-2 Project. March 1999.
!  Argonne National Laboratory and University of Minnesota.
!  Brett M. Averick and Jorge J. More'.

!  **********

REAL (dp), PARAMETER :: p5=0.5_dp
REAL (dp), PARAMETER :: two=2.0_dp
REAL (dp), PARAMETER :: four=4.0_dp
REAL (dp), PARAMETER :: six=6.0_dp

LOGICAL   :: feval, geval
INTEGER   :: i, j, k
REAL (dp) :: dvdx, dvdy, ehxhy, flin, fquad, hx, hxhy, hy, pi,  &
             temp, trule, v, vb, vl, vr, vt, xi

! REAL (dp) :: p

! p(xi) = (1 + ecc*COS(xi))**3

!     Initialization.

pi = four*ATAN(one)
hx = two*pi/DBLE(nx+1)
hy = two*b/DBLE(ny+1)
hxhy = hx*hy
ehxhy = ecc*hxhy

!     Compute the lower bound xl for x if task = 'XL'.

IF (task(1:2) == 'XL') THEN
  DO  k = 1, nx*ny
    x(k) = zero
  END DO
  
  RETURN
  
END IF

!     Compute the standard starting point if task = 'XS'.

IF (task(1:2) == 'XS') THEN
  DO  i = 1, nx
    temp = MAX(SIN(DBLE(i)*hx), zero)
    DO  j = 1, ny
      k = nx*(j-1) + i
      x(k) = temp
    END DO
  END DO
  
  RETURN
  
END IF

IF (task(1:1) == 'F' .OR. task(1:2) == 'FG') THEN
  feval = .true.
ELSE
  feval = .false.
END IF
IF (task(1:1) == 'G' .OR. task(1:2) == 'FG') THEN
  geval = .true.
ELSE
  geval = .false.
END IF

!     Compute the function if task = 'F', the gradient if task = 'G', or
!     both if task = 'FG'.

IF (feval) THEN
  fquad = zero
  flin = zero
END IF
IF (geval) THEN
  DO  k = 1, nx*ny
    fgrad(k) = zero
  END DO
END IF

!     Computation of the quadratic part of the function and corresponding
!     components of the gradient over the lower triangular elements.

DO  i = 0, nx
  xi = DBLE(i)*hx
  trule = hxhy*(p(xi, ecc) + p(xi+hx, ecc) + p(xi, ecc))/six
  DO  j = 0, ny
    k = nx*(j-1) + i
    v = zero
    vr = zero
    vt = zero
    IF (i /= 0 .AND. j /= 0) v = x(k)
    IF (i /= nx .AND. j /= 0) vr = x(k+1)
    IF (i /= 0 .AND. j /= ny) vt = x(k+nx)
    dvdx = (vr-v)/hx
    dvdy = (vt-v)/hy
    IF (feval) fquad = fquad + trule*(dvdx**2+dvdy**2)
    IF (geval) THEN
      IF (i /= 0 .AND. j /= 0) fgrad(k) = fgrad(k) - trule*(dvdx/hx+dvdy/hy)
      IF (i /= nx .AND. j /= 0) fgrad(k+1) = fgrad(k+1) + trule*dvdx/hx
      IF (i /= 0 .AND. j /= ny) fgrad(k+nx) = fgrad(k+nx) + trule*dvdy/hy
    END IF
  END DO
END DO

!     Computation of the quadratic part of the function and corresponding
!     components of the gradient over the upper triangular elements.

DO  i = 1, nx + 1
  xi = DBLE(i)*hx
  trule = hxhy*(p(xi, ecc) + p(xi-hx, ecc) + p(xi, ecc))/six
  DO  j = 1, ny + 1
    k = nx*(j-1) + i
    vb = zero
    vl = zero
    v = zero
    IF (i /= nx+1 .AND. j /= 1) vb = x(k-nx)
    IF (i /= 1 .AND. j /= ny+1) vl = x(k-1)
    IF (i /= nx+1 .AND. j /= ny+1) v = x(k)
    dvdx = (v-vl)/hx
    dvdy = (v-vb)/hy
    IF (feval) fquad = fquad + trule*(dvdx**2+dvdy**2)
    IF (geval) THEN
      IF (i <= nx .AND. j > 1) fgrad(k-nx) = fgrad(k-nx) - trule*dvdy/hy
      IF (i > 1 .AND. j <= ny) fgrad(k-1) = fgrad(k-1) - trule*dvdx/hx
      IF (i <= nx .AND. j <= ny) fgrad(k) = fgrad(k) + trule*(dvdx/hx+dvdy/hy)
    END IF
  END DO
END DO

!     Computation of the linear part of the function and
!     corresponding components of the gradient.

DO  i = 1, nx
  temp = SIN(DBLE(i)*hx)
  IF (feval) THEN
    DO  j = 1, ny
      k = nx*(j-1) + i
      flin = flin + temp*x(k)
    END DO
  END IF
  IF (geval) THEN
    DO  j = 1, ny
      k = nx*(j-1) + i
      fgrad(k) = fgrad(k) - ehxhy*temp
    END DO
  END IF
END DO

!     Finish off the function.

IF (feval) f = p5*fquad - ehxhy*flin

RETURN
END SUBROUTINE dpjbfg




FUNCTION p(xi, ecc) RESULT(fn_val)
! This is the statement function which was part of routine DPJBFG.

REAL (dp), INTENT(IN) :: xi, ecc
REAL (dp)             :: fn_val

fn_val = (1 + ecc*COS(xi))**3
RETURN
END FUNCTION p




SUBROUTINE dmsabc(nx, ny, bottom, top, left, right)
 
! Code converted using TO_F90 by Alan Miller
! Date: 1999-06-29  Time: 15:49:10

INTEGER, INTENT(IN)     :: nx
INTEGER, INTENT(IN)     :: ny
REAL (dp), INTENT(OUT)  :: bottom(:)    ! bottom(nx+2)
REAL (dp), INTENT(OUT)  :: top(:)       ! top(nx+2)
REAL (dp), INTENT(OUT)  :: left(:)      ! left(ny+2)
REAL (dp), INTENT(OUT)  :: right(:)     ! right(ny+2)

!  **********

!  Subroutine dmsabc

!  This subroutine computes Enneper's boundary conditions for the minimal
!  surface area problem on the unit square centered at the origin.

!  The subroutine statement is

!    subroutine dmsabc(nx, ny, bottom, top, left, right)

!  where

!    nx is an integer variable.
!      On entry nx is the number of grid points in the first coordinate
!         direction.
!      On exit nx is unchanged.

!    ny is an integer variable.
!      On entry ny is the number of grid points in the second
!         coordinate direction.
!      On exit ny is unchanged.

!    bottom is a REAL (dp) array of dimension nx + 2.
!      On entry bottom need not be specified.
!      On exit bottom contains boundary values for the bottom
!         boundary of the domain.

!    top is a REAL (dp) array of dimension nx + 2.
!      On entry top need not be specified.
!      On exit top contains boundary values for the top boundary of the domain.
!    left is a REAL (dp) array of dimension ny + 2.
!      On entry left need not be specified.
!      On exit left contains boundary values for the left boundary
!         of the domain.

!    right is a REAL (dp) array of dimension ny + 2.
!      On entry right need not be specified.
!      On exit right contains boundary values for the right boundary
!         of the domain.

!  MINPACK-2 Project. March 1999.
!  Argonne National Laboratory and University of Minnesota.
!  Brett M. Averick.

!  **********

REAL (dp), PARAMETER :: two=2.0_dp, three=3.0_dp, tol=1.0D-10, b=-.50_dp, &
                        t=.50_dp, l=-.50_dp, r=.50_dp
INTEGER, PARAMETER   :: maxit=5

INTEGER   :: i, j, k, limit
REAL (dp) :: det, fnorm, hx, hy, xt, yt
REAL (dp) :: nf(2), njac(2,2), u(2)

!     Compute Enneper's boundary conditions: bottom, top, left, then
!     right.  Enneper's boundary values are obtained by defining
!     bv(x,y) = u**2 - v**2 where u and v are the unique solutions of
!     x = u + u*(v**2) - (u**3)/3, y = -v - (u**2)*v + (v**3)/3.