routines.f90

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

                     deriv_type,eps_vector,num_derivs,  &
                     conv1,conv2,conv3,FRET,nerr)
implicit none
real(8), dimension(:)       ::a,eps_vector
integer, intent(in)         ::HESSINreset,itmax,num_derivs
integer, dimension(:)       ::spare_params
integer  maxnumas,deriv_type
real(8), dimension(maxnumas,size(a))     ::as
real(8), dimension(maxnumas)             ::f
real(8), dimension(size(a))              ::XI,G,DG,HDG
real(8), dimension(size(a),size(a))      ::HESSIN
real(8), dimension(0:maxnumas)           ::x
real(8), intent(in)                      ::schlength
real(8)  conv1,conv2,conv3,FRET,fp,norm
real(8)  x0,x1,x2,xx,f1,alpha,FAC,FAD,FAE,sch,mag
real(8)  line_min_tol,eps,small
integer  nerr,HESrescount,i
integer  iterations,ii,sizea,d_type,n_derivs
logical  reset_d_type

include 'blas.inc'   

interface
   SUBROUTINE func_name(spare_params,as,f,nerr)
     integer, dimension(:)      ::spare_params
     real(8), dimension(:)      ::f
     integer  nerr
     real(8), dimension(:,:)    ::as
   END SUBROUTINE func_name
end interface

!  check for valid derivative type request
   if(deriv_type.lt.0 .or. deriv_type.gt.2) then
      write(*,*) 'invalid derivative type specified'
      write(*,*) 'program stopped'
      stop
   end if

   if(size(eps_vector).ne.size(a)) then
      write(*,*)'incorrect dimensioning between parameter'
      write(*,*)'vector and eps_vector:  PROGRAM STOPPED'
      STOP
   end if

!  set value for 'a small number' relative to 1.0d0
   eps=epsilon(0.0d0)
   small=sqrt(eps)

!  set line minimisation tolerance
   line_min_tol=2.0d0*conv2+eps

!  set first linear partitioning factor
   sch=max((schlength-conv2)/real(maxnumas-1,           &
   kind(0.0d0)),conv2)

!  stretch factor for search length
   mag=real(max(min(HESSINreset,10),2),kind(0.0d0))

!  reset error flag, hessian reset counter
   nerr=0
   HESrescount=0
   iterations=0
   d_type=0
   reset_d_type=.true.
   n_derivs=1
   sizea=size(a)
   
   as(1,:)=a
   CALL func_name(spare_params,as(1:1,:),f(1:1),nerr)
   fret=f(1)
   fp=fret      ! initialise fp
   write(*,*)'initial f = ',fret
   write(*,*)'a = ',a
5  HESSIN=identity(sizea) ! set hessin to identity
6  CALL Dfunc(a,fret,G)
   xi=-G  ! set first search direction (steepest decent)
   do     ! start loop for minimisation
      iterations=iterations+1
      HESrescount=HESrescount+1 !inc Hessian reset counter
      norm=max(norm2(xi),eps)
      xi=xi/norm  ! normalise search direction
      x(0)=0.0d0
      x(1)=conv2                       ! set partition
      as(1,:)=a+x(1)*xi                ! along search
      do i=2,maxnumas                  ! direction
         x(i)=x(i-1)+sch
         as(i,:)=a+x(i)*xi
      end do
      CALL func_name(spare_params,as,f,nerr)
      if (f(1).ge.fret) then           ! see if  at a
         alpha = 0.0d0                 ! true minimum
         fp = fret
         GOTO 99
      end if

10    do i=2,maxnumas
         if(f(i).ge.f(i-1)) then  ! know f(1)<fret
            x2=x(i)    ! Now we have found a
            x1=x(i-1)  ! point with lower function
            f1=f(i-1)  ! value. Now find a point further
             x0=x(i-2) ! on in the search direction which
             goto 20   ! gives a function value higher
         end if        ! than the previous point, you
      end do           ! have now bracketed a minimum.

      f(1)=f(maxnumas)              ! if have not found
      x(0:1)=x(maxnumas-1:maxnumas) ! 'a' then look
      as(1,:)=as(maxnumas,:)        ! further along search
      do i=2,maxnumas               ! direction
         x(i)=x(i-1)+sch
         as(i,:)=a+x(i)*xi
      end do
      CALL func_name(spare_params,as(2:,:),f(2:),nerr)
      go to 10

20    if(abs(x2-x0).le.line_min_tol) goto 30  ! corner
      xx=(x2-x0)/real(maxnumas-1,kind(0.0d0)) ! minimum
      x(0:1)=x0                               ! like a
      as(1,:)=a+x(1)*xi                       ! scared
      do i=2,maxnumas                         ! animal
         x(i)=x(i-1)+xx
         as(i,:)=a+x(i)*xi
      end do
      CALL func_name(spare_params,as,f,nerr)

      do i=2,maxnumas-1
         if(f(i).ge.f(i-1)) then
            x2=x(i)
            x1=x(i-1)
            f1=f(i-1)
            x0=x(i-2)
            goto 20
         end if
      end do

      x2=x(maxnumas)   ! if there is round off error
      x1=x(maxnumas-1) ! on iteration MAXNUMAS
      f1=f(maxnumas-1) ! then it will exit the above
      x0=x(maxnumas-2) ! do loop anyway, Therefore only
      goto 20          ! count to maxnumas-1 and go here.

30    fp=fret
      fret=f1
      alpha=x1
      xi=alpha*xi
      a=a+xi

!     Assuming that iterations get closer to minimum lets
!     decrease search length so we get a finer partition.
!     Use mag*last search length as an estimate.
      sch=max((alpha*mag+conv2)/real(maxnumas-1,        &
               kind(0.0d0)),conv2)

40    if (alpha.lt.eps) goto 99        ! see if converged
      if (abs(fret - fp).lt.eps) goto 99 
      if ((abs((fret-fp)/fp).le.conv1).and.             &
          (alpha.le.conv2).and.(norm2(G).le.conv3))     &
          goto 99 

      DG=G
      CALL Dfunc(a,fret,G)
      n_derivs=1
      reset_d_type=.true.
      write(*,*)'iteration',iterations,'fret = ',fret
      write(*,*)'a = ',a
      open(1,file='a.vec',status='old')
      rewind(1)
      do ii=1,sizea
         write(1,*) a(ii)
      end do
      close(1)

50    if (HESrescount.ge.HESSINreset) then   ! test for
         HESrescount = 0                     ! clearing
         HESSIN=identity(sizea)              ! hessian
         xi = -G
         goto 60
      endif 

!     calculate next search direction
      xi=xi*norm
      DG=G-DG
      HDG=matmul(HESSIN,DG)
      fac=dot_product(DG,xi)
      fae=dot_product(DG,HDG)
      if((abs(fac).le.small).or.(abs(fae).le.small)) then
         HESrescount=HESSINreset
         goto 50
      end if
      fac = 1.0d0/fac 
      fad = 1.0d0/fae
      DG=fac*xi-fad*HDG
      HESSIN=HESSIN+fac*outer_product(xi,xi)-fad*       &
             outer_product(HDG,HDG)+fae*                &
             outer_product(DG,DG)
      xi=-matmul(HESSIN,G)

60    if(iterations.eq.itmax) exit
   end do                ! end of iteration loop
   nerr = - 1              ! flag too many iterations
   return

99 if(HESrescount.ge.2) then
      HESrescount=0
      goto 5
   end if

   HESrescount=0
   if(d_type.lt.deriv_type) then
      d_type=d_type+1
      goto 6
   end if

   if(reset_d_type) then
      d_type=0
      reset_d_type=.false.
      goto 6
   end if

   if(n_derivs.lt.num_derivs) then
      n_derivs=n_derivs+1
      d_type=0
      reset_d_type=.true.
      goto 6
   end if

   conv1 = abs (fret - fp) 
   conv2 = alpha 
   d_type=deriv_type
   n_derivs=1
   call dfunc(a,fret,G)
   conv3 = norm2(G)
RETURN 

CONTAINS
   SUBROUTINE dfunc(a,fret,G)

!  derivative of func_name for minimisation routine
   IMPLICIT none
   real(8), dimension(:)                 :: a
   real(8), dimension(4*size(a),size(a)) :: as
   real(8), dimension(:)                 :: G
   real(8), dimension(size(a))           :: delta
   real(8), dimension(4*size(a))         :: term
   integer  i, sizea
   real(8)  fret
  
   sizea=size(a)
!  set up as
   as=spread(a,dim=1,ncopies=sizea*4)
   where(eps_vector>0.0d0) delta=max(eps_vector*abs(a), &
                                     small)
   where(eps_vector<0.0d0) delta=max(-eps_vector,small)
   where(eps_vector.eq.0.0d0) delta=max(small*abs(a),   &
                                        small)
   delta=delta*real(n_derivs,kind(0.0d0))
   do i=1,sizea
      as(i,i)=a(i)-delta(i)
      as(sizea+i,i)=a(i)+delta(i)
      as(2*sizea+i,i)=a(i)-2.0d0*delta(i)
      as(3*sizea+i,i)=a(i)+2.0d0*delta(i)
   end do
  
   if (d_type.eq.0) then        ! forward approximation
      CALL func_name(spare_params,as(sizea+1:2*sizea,   &
                     :),term(sizea+1:2*sizea),nerr)
      g=(term(sizea+1:sizea+sizea)-fret)/delta
   end if
     
   if(d_type.eq.1) then       ! symmetric approximation
      CALL func_name(spare_params,as(1:2*sizea,:),      &
                     term(1:2*sizea),nerr)
     g=(term(sizea+1:sizea+sizea)-term(1:sizea))/       &
       (2.0d0*delta)
   end if
     
   if(d_type.eq.2) then ! stable symmetric approximation
      CALL func_name(spare_params,as(1:4*sizea,:),      &
                     term(1:4*sizea),nerr)
      g=(term(2*sizea+1:3*sizea)-term(3*sizea+1:        &
        4*sizea) + 8.0d0*( term(sizea+1:2*sizea)-       &
        term(1:sizea)))/(12.0d0*delta)
   end if
   RETURN 
   END SUBROUTINE dfunc

END SUBROUTINE MINIMIZE




SUBROUTINE N01RNG(seed,random)
implicit none
real(8), dimension(:)          ::random
integer, dimension(:)          ::seed
integer, dimension(size(seed)) ::seed1
real(8), dimension(min(size(random),1000))          ::W
real(8), dimension(2*size(W))                       ::V
logical, dimension(size(W))                         ::flag
integer N,i,M,number,dim1,dim2,dim3

dim1=size(W)
dim2=dim1+1
dim3=dim1+dim1
number=size(random)
call random_seed(put=seed)
i=0
do
   call random_number(V)
   V=2.0d0*V-1.0d0
   W=V(1:dim1)*V(1:dim1) + V(dim2:dim3)*V(dim2:dim3)
   flag=.false.
   where (W <= 1.0d0 .and. W > 0.0d0) 
      W=log(W)/W
      W=sqrt(-W-W)
      V(1:dim1)=V(1:dim1)*W
      V(dim2:dim3)=V(dim2:dim3)*W
      flag=.true.  !true means it is part of N(0,1) series
   endwhere
   N=count(flag)
   V(1:N)=pack(V(1:dim1),flag)
   V(N+1:N+N)=pack(V(dim2:dim3),flag)
   M=min(N+N,number-i)
   random(i+1:i+M)=V(1:M)
   i=i+M
   if(i.ge.number) exit
end do
call random_seed(get=seed1)
seed=seed1
return

END SUBROUTINE N01RNG