routines.f90

这是一个实数编码的遗传算法

RECURSIVE FUNCTION DIFFER(func_name,as,param_no,        &
                          eps_of_param_no) result(diff)
implicit none
real(8), dimension(:,:)            ::as
real(8), dimension(size(as,dim=1)) ::eps,diff,f1,f2,    &
                                     f3,f4,tmpas
real(8)  small,eps_of_param_no
integer  no_as,param_no

interface
  subroutine func_name(as,f)
     real(8), dimension(:,:)            ::as
     real(8), dimension(size(as,dim=1)) ::f
  end subroutine func_name
end interface

no_as = size(as,dim=1)
small = sqrt(epsilon(0.0d0))

if((param_no.gt.size(as,dim=2)).or.(param_no.lt.1)) then
   write(*,*)'invalid parameter number in DIFFER'
   write(*,*)'PROGRAM STOPPED'
   stop
end if
if(eps_of_param_no.gt.0.0d0) EPS(:)=max(                &
     eps_of_param_no*abs(as(:,param_no)),small)
if(eps_of_param_no.lt.0.0d0) EPS(:)=max(                &
     -eps_of_param_no,small)
if(eps_of_param_no.eq.0.0d0) EPS(:)=max(                &
     small*abs(as(:,param_no)),small)
tmpas=as(:,param_no)
as(:,param_no)=as(:,param_no)-eps-eps
call func_name(as,f1)
as(:,param_no)=as(:,param_no)+eps
call func_name(as,f2)
as(:,param_no)=as(:,param_no)+eps+eps
call func_name(as,f3)
as(:,param_no)=as(:,param_no)+eps
call func_name(as,f4)
as(:,param_no)=tmpas
Diff=( f1-f4 + 8.0d0*(f3-f2))/ (12.0d0*eps(:))

END FUNCTION DIFFER



FUNCTION DIFFERENTIATE(func_name,as,param_no,           &
                       eps_of_param_no)
implicit none
real(8), dimension(:,:)                             ::as
real(8), dimension(4*size(as,dim=1),size(as,dim=2)) ::as1
real(8), dimension(size(as,dim=1))    ::eps,differentiate
real(8), dimension(4*size(as,dim=1))  ::f
real(8)  small,eps_of_param_no
integer  k,no_as,param_no,point_1,point_2,point_4,point_5

interface
  subroutine func_name(x,f)
     real(8), dimension(:,:)           ::x
     real(8), dimension(size(x,dim=1)) ::f
  end subroutine func_name
end interface

no_as = size(as,dim=1)
small  = sqrt(epsilon(0.0d0))
point_5=no_as*4
point_4=no_as*3
point_2=no_as*2
point_1=no_as
if((param_no.gt.size(as,dim=2)).or.(param_no.lt.1)) then
   write(*,*)'invalid parameter number in DIFFERENTIATE'
   write(*,*)'PROGRAM STOPPED'
   stop
end if

do k=1,4
   as1((k-1)*no_as+1:k*no_as,:)=as
end do
if(eps_of_param_no.gt.0.0d0) EPS(:)=max(                &
             eps_of_param_no*abs(as(:,param_no)),small)
if(eps_of_param_no.lt.0.0d0) EPS(:)=max(                &
             -eps_of_param_no,small)
if(eps_of_param_no.eq.0.0d0) EPS(:)=max(                &
             small*abs(as(:,param_no)),small)
as1(1:point_1,param_no) = as(:,param_no)-eps-eps
as1(point_1+1:point_2,param_no) = as(:,param_no)-eps
as1(point_2+1:point_4,param_no) = as(:,param_no)+eps
as1(point_4+1:point_5,param_no) = as(:,param_no)+eps+eps

call func_name(as1,f)
Differentiate=(f(1:point_1)-f(point_4+1:point_5) +      &
               8.0d0*(f(point_2+1:point_4)-             &
               f(point_1+1:point_2)))/(12.0d0*eps(:))

END FUNCTION DIFFERENTIATE

RECURSIVE FUNCTION INTEG(func_name,as,param_no,lower,   &
                         upper,half_no_partitions)      &
                         result(intege)
implicit none
integer  half_no_partitions,param_no
real(8), dimension(:,:)           ::as
real(8), dimension(size(as,dim=1))::eps,intege,lower,   &
                                    upper,tmp,tmpas,f1,f2
integer  k,no_as,N

interface
  subroutine func_name(as,f)
     real(8), dimension(:,:)            ::as
     real(8), dimension(size(as,dim=1)) ::f
  end subroutine func_name
end interface

if((param_no.gt.size(as,dim=2)).or.(param_no.lt.1)) then
   write(*,*)'invalid parameter number in INTEG'
   WRITE(*,*)'PROGRAM STOPPED'
   stop
end if

no_as  = size(as,dim=1)
N      = 2*half_no_partitions+1
EPS(:) = (upper(:)-lower(:))/real(N-1,kind(0.0d0))
tmpas  = as(:,param_no)
tmp    = lower

as(:,param_no)=tmp
call func_name(as,f1)
as(:,param_no)=as(:,param_no)+eps
call func_name(as,f2)
Intege = f1 + 4.0d0*f2
do k=3,N-1,2
   as(:,param_no)=as(:,param_no)+eps
   call func_name(as,f1)
   as(:,param_no)=as(:,param_no)+eps
   call func_name(as,f2)
   Intege = Intege + 2.0d0*f1 + 4.0d0*f2
end do
as(:,param_no)=as(:,param_no)+eps
call func_name(as,f1)
Intege = (Intege + f1)*eps/3.0d0
as(:,param_no)=tmpas

END FUNCTION INTEG



FUNCTION INTEGRATE(func_name,as,param_no,lower,upper,   &
                   half_no_partitions)
implicit none
integer  half_no_partitions,param_no
real(8), dimension(:,:)                         ::as
real(8), dimension((2*half_no_partitions+1)*            &
         size(as,dim=1),size(as,dim=2))         ::as1
real(8), dimension(size(as,dim=1))    ::eps,integrate,  &
                                        lower,upper,tmp
real(8), dimension((2*half_no_partitions+1)*            &
         size(as,dim=1))                        ::f
integer, dimension(0:2*half_no_partitions+1)    ::point
integer  k,no_as,N

interface
  subroutine func_name(x,f)
     real(8), dimension(:,:)           ::x
     real(8), dimension(size(x,dim=1)) ::f
  end subroutine func_name
end interface

no_as  = size(as,dim=1)
N      = 2*half_no_partitions+1
if((param_no.gt.size(as,dim=2)).or.(param_no.lt.1)) then
   write(*,*)'invalid parameter number in INTEGRATE'
   write(*,*)'PROGRAM STOPPED'
   stop
end if

do k=0,N
   point(k)=k*no_as
end do
EPS(:)=(upper(:)-lower(:))/real(N-1,kind(0.0d0))
tmp=lower
do k=1,N
   as1(point(k-1)+1:point(k),:)=as
   as1(point(k-1)+1:point(k),param_no) = tmp
   tmp=tmp+eps
end do

call func_name(as1,f)
Integrate=f(1:point(1))+4.0d0*f(point(1)+1:point(2)) +  &
          f(point(N-1)+1:point(N))
do k=3,N-1,2
   Integrate=Integrate+2.0d0*f(point(k-1)+1:point(k)) + &
                       4.0d0*f(point(k)+1:point(k+1))
end do
Integrate=Integrate*eps/3.0d0
END FUNCTION INTEGRATE



RECURSIVE FUNCTION INTEG_OVER_PARTITION(func_name,as,   &
                              param_no,partition)       &
                              result(integ)
implicit none
integer, dimension(:)                          ::param_no
real(8), dimension(:,:)                        ::as
real(8), dimension(:,:)                        ::partition
real(8), dimension(size(as,dim=1))             ::integ,f

real(8), dimension(size(as,dim=1),size(param_no)) ::tmpas
integer  k,no_as,N,M,i

interface
  subroutine func_name(as,f)
     real(8), dimension(:,:)            ::as
     real(8), dimension(size(as,dim=1)) ::f
  end subroutine func_name
end interface

no_as = size(as,dim=1)
N     = size(partition,dim=2)
M     = size(param_no)
i     = size(as,dim=2)
if(N.eq.1) then  ! need at least two points for an area
   integ = 0.0d0                ! integral is zero
   goto 1
end if
if(M.ne.size(partition,dim=1)) then
   write(*,*)'SIZE(PARTITION,DIM=1)  .NE.  SIZE(PARAM_NO)'
   write(*,*)'IN ROUTINE:        INTEG_OVER_PARTITION'
   write(*,*)'PROGRAM STOPPED'
   stop
end if

do k=1,M
   if((param_no(k).gt.i).or.(param_no(k).lt.1)) then
      write(*,*)'invalid parameter number in'
      write(*,*)'INTEG_OVER_PARTITION'
      write(*,*)'PROGRAM STOPPED'
      stop
   end if
end do

tmpas(:,:)=as(:,param_no)
as(:,param_no)=spread(partition(1:M,1),dim=1,           &
                      ncopies=no_as)
call func_name(as,f)
Integ=f*(partition(1,2)-partition(1,1))
do k=2,N-1
   as(:,param_no)=spread(partition(1:M,k),dim=1,        &
                         ncopies=no_as)
   call func_name(as,f)
   Integ=Integ+f*(partition(1,k+1)-partition(1,k-1))
end do
as(:,param_no)=spread(partition(1:M,N),dim=1,           &
                      ncopies=no_as)
call func_name(as,f)
Integ=(Integ+f*(partition(1,N)-partition(1,N-1)))*0.5d0

1 END FUNCTION INTEG_OVER_PARTITION



FUNCTION INTEGRATE_OVER_PARTITION(func_name,as,         &
                                  param_no,partition)

implicit none
integer, dimension(:)                      ::param_no
real(8), dimension(:,:)                    ::as
real(8), dimension(:,:)                    ::partition
real(8), dimension(size(as,dim=1))         ::           &
                             integrate_over_partition
real(8), dimension(size(partition,dim=2)*               &
         size(as,dim=1),size(as,dim=2))    ::as1
real(8), dimension(size(partition,dim=2)*               &
         size(as,dim=1))                   ::f
integer, dimension(0:size(partition,dim=2))::point
integer  k,no_as,N,M,i

interface
  subroutine func_name(as,f)
     real(8), dimension(:,:)            ::as
     real(8), dimension(size(as,dim=1)) ::f
  end subroutine func_name
end interface

no_as = size(as,dim=1)
N     = size(partition,dim=2)
M     = size(param_no)

if(N.eq.1) then   ! need at least two points for an area
   integrate_over_partition=0.0d0   ! integral is zero
   goto 1
end if

if(M.ne.size(partition,dim=1)) then
   write(*,*)'SIZE(PARTITION,DIM=1) .NE. SIZE(PARAM_NO)'
   write(*,*)'IN ROUTINE  INTEGRATE_OVER_PARTITION'
   write(*,*)'PROGRAM STOPPED'
   stop
end if

if((any(param_no(:).gt.size(as,dim=2))).or.(any(     &
   param_no(:).lt.1))) then
      write(*,*)'invalid parameter number in routine'
      write(*,*)'INTEGRATE_OVER_PARTITION'
      write(*,*)'PROGRAM STOPPED'
     stop
   end if

point(0)=0
do k=1,N
   point(k)=k*no_as
   as1(point(k-1)+1:point(k),:)=as
   do i=1,M
      as1(point(k-1)+1:point(k),param_no(i)) =          &
         partition(i,k)
   end do
end do

call func_name(as1,f)
Integrate_over_partition=f(1:point(1))*(partition(1,2)- &
    partition(1,1)) + f(point(N-1)+1:point(N))*         &
    (partition(1,N) - partition(1,N-1))
do k=2,N-1
   Integrate_over_partition=Integrate_over_partition +  &
       f(point(k-1)+1:point(k))*(partition(1,k+1)-      &
       partition(1,k-1))
end do
Integrate_over_partition=Integrate_over_partition*0.5d0

1  END FUNCTION INTEGRATE_OVER_PARTITION




SUBROUTINE MINIMIZE (spare_params,a,maxnumas,func_name, &
                     HESSINreset,itmax,schlength,       &