dmain.f90

牛顿优化算法源fortran代码

hx = (r-l)/DBLE(nx+1)
hy = (t-b)/DBLE(ny+1)

DO  j = 1, 4
  IF (j == 1) THEN
    yt = b
    xt = l
    limit = nx + 2
  ELSE IF (j == 2) THEN
    yt = t
    xt = l
    limit = nx + 2
  ELSE IF (j == 3) THEN
    yt = b
    xt = l
    limit = ny + 2
  ELSE IF (j == 4) THEN
    yt = b
    xt = r
    limit = ny + 2
  END IF
  
!        Use Newton's method to solve xt = u + u*(v**2) - (u**3)/3,
!        yt = -v - (u**2)*v + (v**3)/3.
  
  DO  i = 1, limit
    u(1) = xt
    u(2) = -yt
    DO  k = 1, maxit
      nf(1) = u(1) + u(1)*u(2)**2 - u(1)**3/three - xt
      nf(2) = -u(2) - u(1)**2*u(2) + u(2)**3/three - yt
      fnorm = SQRT(nf(1)*nf(1)+nf(2)*nf(2))
      IF (fnorm <= tol) EXIT
      njac(1,1) = one + u(2)**2 - u(1)**2
      njac(1,2) = two*u(1)*u(2)
      njac(2,1) = -two*u(1)*u(2)
      njac(2,2) = -one - u(1)**2 + u(2)**2
      det = njac(1,1)*njac(2,2) - njac(1,2)*njac(2,1)
      u(1) = u(1) - (njac(2,2)*nf(1)-njac(1,2)*nf(2))/det
      u(2) = u(2) - (njac(1,1)*nf(2)-njac(2,1)*nf(1))/det
    END DO
    
    IF (j == 1) THEN
      bottom(i) = u(1)*u(1) - u(2)*u(2)
      xt = xt + hx
    ELSE IF (j == 2) THEN
      top(i) = u(1)*u(1) - u(2)*u(2)
      xt = xt + hx
    ELSE IF (j == 3) THEN
      left(i) = u(1)*u(1) - u(2)*u(2)
      yt = yt + hy
    ELSE IF (j == 4) THEN
      right(i) = u(1)*u(1) - u(2)*u(2)
      yt = yt + hy
    END IF
  END DO
END DO

RETURN
END SUBROUTINE dmsabc



SUBROUTINE dmsafg(nx, ny, x, f, fgrad, task, bottom, top, left, right)
 
! 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)          :: bottom(:)   ! bottom(nx+2)
REAL (dp), INTENT(IN)          :: top(:)      ! top(nx+2)
REAL (dp), INTENT(IN)          :: left(:)     ! left(ny+2)
REAL (dp), INTENT(IN)          :: right(:)    ! right(ny+2)

!  **********

!  Subroutine dmsafg

!  This subroutine computes the function and gradient of the
!  minimal surface area problem.

!  The subroutine statement is

!    subroutine dmsafg(nx, ny, x, f, fgrad, task, 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.

!    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.

!      On exit task is unchanged.

!    bottom is a REAL (dp) array of dimension nx + 2.
!      On entry bottom must contain boundary data beginning
!         with the lower left corner of the domain.
!      On exit bottom is unchanged.

!    top is a REAL (dp) array of dimension nx + 2.
!      On entry top must contain boundary data beginning with
!         the upper left corner of the domain.
!      On exit top is unchanged.

!    left is a REAL (dp) array of dimension ny + 2.
!      On entry left must contain boundary data beginning with
!         the lower left corner of the domain.
!      On exit left is unchanged.

!    right is a REAL (dp) array of dimension ny + 2.
!      On entry right must contain boundary data beginning with
!         the lower right corner of the domain.
!      On exit right is unchanged.

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

!  **********

REAL (dp), PARAMETER :: p5=0.5_dp, two=2.0_dp

LOGICAL   :: feval, geval
INTEGER   :: i, j, k
REAL (dp) :: alphaj, area, betai, dvdx, dvdy, fl, fu, hx, hy,  &
             v, vb, vl, vr, vt, xline, yline

!     Initialize.

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

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

IF (task(1:2) == 'XS') THEN
  DO  j = 1, ny
    alphaj = DBLE(j)*hy
    DO  i = 1, nx
      k = nx*(j-1) + i
      betai = DBLE(i)*hx
      yline = alphaj*top(i+1) + (one-alphaj)*bottom(i+1)
      xline = betai*right(j+1) + (one-betai)*left(j+1)
      x(k) = (yline+xline)/two
    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) f = zero
IF (geval) THEN
  DO  k = 1, nx*ny
    fgrad(k) = zero
  END DO
END IF

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

DO  j = 0, ny
  DO  i = 0, nx
    k = nx*(j-1) + i
    IF (i >= 1 .AND. j >= 1) THEN
      v = x(k)
    ELSE
      IF (j == 0) v = bottom(i+1)
      IF (i == 0) v = left(j+1)
    END IF
    IF (i < nx .AND. j > 0) THEN
      vr = x(k+1)
    ELSE
      IF (i == nx) vr = right(j+1)
      IF (j == 0) vr = bottom(i+2)
    END IF
    IF (i > 0 .AND. j < ny) THEN
      vt = x(k+nx)
    ELSE
      IF (i == 0) vt = left(j+2)
      IF (j == ny) vt = top(i+1)
    END IF
    dvdx = (vr-v)/hx
    dvdy = (vt-v)/hy
    fl = SQRT(one+dvdx**2+dvdy**2)
    IF (feval) f = f + fl
    IF (geval) THEN
      IF (i >= 1 .AND. j >= 1) fgrad(k) = fgrad(k) - (dvdx/hx+dvdy/hy)/fl
      IF (i < nx .AND. j > 0) fgrad(k+1) = fgrad(k+1) + (dvdx/hx)/fl
      IF (i > 0 .AND. j < ny) fgrad(k+nx) = fgrad(k+nx) + (dvdy/hy)/fl
    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
    IF (i <= nx .AND. j > 1) THEN
      vb = x(k-nx)
    ELSE
      IF (j == 1) vb = bottom(i+1)
      IF (i == nx+1) vb = right(j)
    END IF
    IF (i > 1 .AND. j <= ny) THEN
      vl = x(k-1)
    ELSE
      IF (j == ny+1) vl = top(i)
      IF (i == 1) vl = left(j+1)
    END IF
    IF (i <= nx .AND. j <= ny) THEN
      v = x(k)
    ELSE
      IF (i == nx+1) v = right(j+1)
      IF (j == ny+1) v = top(i+1)
    END IF
    dvdx = (v-vl)/hx
    dvdy = (v-vb)/hy
    fu = SQRT(one+dvdx**2+dvdy**2)
    IF (feval) f = f + fu
    IF (geval) THEN
      IF (i <= nx .AND. j > 1) fgrad(k-nx) = fgrad(k-nx) - (dvdy/hy)/fu
      IF (i > 1 .AND. j <= ny) fgrad(k-1) = fgrad(k-1) - (dvdx/hx)/fu
      IF (i <= nx .AND. j <= ny) fgrad(k) = fgrad(k) + (dvdx/hx+dvdy/hy)/fu
    END IF
  END DO
END DO

!     Scale the function and the gradient.

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

RETURN
END SUBROUTINE dmsafg



SUBROUTINE dodcfg(nx, ny, x, f, fgrad, task, lambda)
 
! 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)          :: lambda

!  **********

!  Subroutine dodcfg

!  This subroutine computes the function and gradient of the
!  optimal design with composite materials problem.

!  The subroutine statement is

!    subroutine dodcfg(nx, ny, x, f, fgrad, task, lambda)

!  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.

!      On exit task is unchanged.

!    lambda is a REAL (dp) variable.
!      On entry lambda is the Lagrange multiplier.
!      On exit lambda is unchanged.

!  Subprograms called

!    MINPACK-supplied   ...   dodcps

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

!  **********

REAL (dp), PARAMETER :: p5=0.5_dp, two=2.0_dp, mu1=one, mu2=two

LOGICAL   :: feval, geval
INTEGER   :: i, j, k
REAL (dp) :: area, dpsi, dpsip, dvdx, dvdy, gradv, hx, hxhy,  &
             hy, temp, t1, t2, v, vb, vl, vr, vt

! EXTERNAL dodcps

!     Initialization.

hx = one/DBLE(nx+1)
hy = one/DBLE(ny+1)
hxhy = hx*hy
area = p5*hxhy

!     Compute the break points.

t1 = SQRT(two*lambda*mu1/mu2)
t2 = SQRT(two*lambda*mu2/mu1)

!     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))**2
    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) f = zero
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 (j >= 1 .AND. i >= 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
    gradv = dvdx**2 + dvdy**2
    IF (feval) THEN
      CALL dodcps(gradv, mu1, mu2, t1, t2, dpsi, 0, lambda)
      f = f + dpsi
    END IF
    IF (geval) THEN
      CALL dodcps(gradv, mu1, mu2, t1, t2, dpsip, 1, lambda)
      IF (i >= 1 .AND. j >= 1)  &
          fgrad(k) = fgrad(k) - two*(dvdx/hx+dvdy/hy)*dpsip
      IF (i < nx .AND. j > 0) fgrad(k+1) = fgrad(k+1) + two*(dvdx/hx)*dpsip
      IF (i > 0 .AND. j <