1425阅读
0回复

大家看看这个程序，FORTRAN要转化为C

UID ：41884

注册：2009-09-17
登录：2010-01-13
发帖：2
等级：旁观者

0楼发表于: 2010-01-05 18:47:35

c    ****************************************************************
c     DIFFERENTIAL EVOLUTION STRATEGY(DES) for minization problem
c     Last updated on September 28, 2002
c    ================================================================
c     This program is workable under FORTRAN PowerStation and FORTRAN
c     90/95 on Planet workstations at Temasek Laboratories. To run
c     the program under FORTRAN 90/95 on Planet workstations, the
c     program should be compiled using the following command
c     f95 -YEXT_NAMES=LCS -YEXT_SFX=_ -L/opt/absoft/lib -limsl -limslblas
c     -lV77 -o des des.f
c    ================================================================
      PROGRAM MAIN
c     for FORTRAN PowerStation
      USE PORTLIB
      implicit none
      real*4 current_time
c     for FORTRAN 90/95 on Planet workstations at Temasek Laboratories
c      implicit none
c      real*4 secnds,current_time
      real*8 eps
      integer max_generation
      real*8 crossover_probability,mutation_intensity
      integer length_chromosome,population_size
      real*8, allocatable :: genes_min(:),genes_max(:)
      character*80 input_file,objf_file,result_file,note
      integer j

c     ask for input and output file names
      write(*,*) 'Parameter file name, please.'
      read(*,*) input_file
      write(*,*) 'File name for object function value, please.'
      read(*,*) objf_file
      write(*,*) 'File name for variables, please.'
      read(*,*) result_file
      write(*,*)
c     read control parameters
      open(1,file=input_file)
      read(1,*) note,eps
      read(1,*) note,max_generation
      read(1,*) note,crossover_probability
      read(1,*) note,mutation_intensity
      read(1,*) note,length_chromosome
      read(1,*) note,population_size
      allocate(genes_min(length_chromosome))
      allocate(genes_max(length_chromosome))
      do j=1,length_chromosome
        read(1,*) note,genes_min(j),genes_max(j)
      end do
      close(1)
c     Search the optimal result
      write(*,*) 'Search in progress, please wait.'
      current_time=secnds(0.0)
      call DES(objf_file,result_file,eps,max_generation,
     +         crossover_probability,mutation_intensity,
     +         length_chromosome,population_size,
     +         genes_min,genes_max)
      write(*,*) 'Search over'
      write(*,*) 'Search takes ',secnds(current_time),' s'
      open(3,file=result_file,status='unknown',access='append')
      write(3,*) 'Search takes ',secnds(current_time),' s'
      close(3)
c     remind the output file names
      write(*,*) 'Please view the brief results in ',objf_file
      write(*,*) 'Please view the detailed results in ',result_file
      deallocate(genes_min,genes_max)
      end
c    ****************************************************************
c     Differential Evolution Strategy
c    ================================================================
      SUBROUTINE DES(objf_file,result_file,eps,max_generation,
     +               crossover_probability,mutation_intensity,
     +               length_chromosome,population_size,
     +               genes_min,genes_max)
c     for FORTRAN PowerStation
      USE PORTLIB
      implicit none
      real*4 current_time
c     for FORTRAN 90/95 on Planet workstations at Temasek Laboratories
c      implicit none
c      real*4 secnds,current_time
      character*80 objf_file,result_file
      real*8 eps
      integer max_generation
      real*8 crossover_probability,mutation_intensity
      integer length_chromosome,population_size
      real*8 genes_min(length_chromosome)
      real*8 genes_max(length_chromosome)
c    ================================================================
      integer minpp,maxpp,i,j,gen
      real*8, allocatable :: new_chromosome(:,:),new_objf(:)
      real*8, allocatable :: old_chromosome(:,:),old_objf(:)
c     record control parameters
      open(2,file=objf_file)
      open(3,file=result_file)
      write(3,*) 'operation parameters'
      write(3,*) 'Convergence criteria ',real(eps)
      write(3,*) 'Maximum number of evolution generations ',
     +            max_generation
      write(3,*) 'crossover probability ',real(crossover_probability)
      write(3,*) 'mutation intensity ',real(mutation_intensity)
      write(3,*) 'Number of optimizaztion parameters ',length_chromosome
      write(3,*) 'population size ',population_size
      write(3,*)
      write(3,*) 'search range'
      do i=1,length_chromosome
        write(3,70) i,genes_min(i),genes_max(i)
      end do
      write(3,*)
c     allocate memory for population
      allocate(new_chromosome(population_size,length_chromosome))
      allocate(new_objf(population_size))
      allocate(old_chromosome(population_size,length_chromosome))
      allocate(old_objf(population_size))
c     initialize the random number generator
      call RANDOMIZE()

c     generating the initial population
      gen=0
      current_time=secnds(0.0)
      write(*,*) 'Generating initial population'
      call INITIALIZE(population_size,length_chromosome,
     +                genes_min,genes_max,new_chromosome,new_objf)
      call STATISTICS(population_size,new_objf,minpp,maxpp)
      write(*,*) 'Minimum objective function value ',new_objf(minpp)
      write(*,*) 'Maximum objective function value ',new_objf(maxpp)
      write(2,*) gen,new_objf(minpp)
      write(3,80) gen,new_objf(minpp)
      do 10 j=1,length_chromosome
        write(3,90) j,new_chromosome(minpp,j)
10    end do
      write(3,*)
      write(*,*) 'Initial population takes ',secnds(current_time)
      write(*,*)
c     searching for optimal solution by evolution
20    current_time=secnds(0.0)
      gen=gen+1
      write(*,*) 'Evolution to generation ',gen
c     update the old population
      do 40 i=1,population_size
        do 30 j=1,length_chromosome
          old_chromosome(i,j)=new_chromosome(i,j)
30      end do
        old_objf(i)=new_objf(i)
40    end do
c     generate new population
      call EVOLUTION(crossover_probability,mutation_intensity,
     +               length_chromosome,population_size,
     +               genes_min,genes_max,
     +               old_chromosome,old_objf,
     +               new_chromosome,new_objf)
      call STATISTICS(population_size,new_objf,minpp,maxpp)
      write(*,*) 'Minimum objective function value ',new_objf(minpp)
      write(*,*) 'Maximum objective function value ',new_objf(maxpp)
      write(2,*) gen,new_objf(minpp)
      write(3,80) gen,new_objf(minpp)
      do 50 j=1,length_chromosome
        write(3,90) j,new_chromosome(minpp,j)
50    end do
      write(3,*)
      write(*,*) 'Evolution to generation ',gen,' takes ',
     +            secnds(current_time)
      write(*,*)
c     convergence examination
      if(new_objf(minpp).lt.eps) then
        write(*,*) 'Minimum in capture, congraduation!'
        go to 60
      else
        if(gen.ge.max_generation) then
          write(*,*) 'Maximum iterations use up.'
          go to 60
        end if
        go to 20
      end if
60    deallocate(new_chromosome,new_objf)
      deallocate(old_chromosome,old_objf)
      close(2)
      close(3)
70    format(1x,'parameter ',i3,' from ',f8.2,' to ',f8.2)
80    format(1x,'generation ',i3,' minimum objective function value ',
     +       g10.3)
90    format('x(',i2,')=',g10.3)
      return
      end
c    ****************************************************************
c     generate the initial population.
c    ================================================================
      SUBROUTINE INITIALIZE(population_size,length_chromosome,
     +                      genes_min,genes_max,
     +                      new_chromosome,new_objf)
      implicit none
      integer population_size,length_chromosome
      real*8 genes_min(length_chromosome)
      real*8 genes_max(length_chromosome)
      real*8 new_chromosome(population_size,length_chromosome)
      real*8 new_objf(population_size)
c    ================================================================
      integer i,j
      real*8 RANDOM_REAL,OBJF
      real*8, allocatable :: chromosome(:)
c     initialize the random number generator
      call RANDOMIZE()
      allocate(chromosome(length_chromosome))
      do 20 i=1,population_size
        do 10 j=1,length_chromosome
          chromosome(j)=genes_min(j)
     +                 +RANDOM_REAL()*(genes_max(j)-genes_min(j))
          new_chromosome(i,j)=chromosome(j)
10      end do
        new_objf(i)=OBJF(length_chromosome,chromosome)
20    end do
      deallocate(chromosome)
      return
      end
c    ****************************************************************
c     generate a new population by performing the genetic operations,
c     i.e., crossover, mutation and selection.
c    ================================================================
      SUBROUTINE EVOLUTION(crossover_probability,mutation_intensity,
     +                     length_chromosome,population_size,
     +                     genes_min,genes_max,
     +                     old_chromosome,old_objf,
     +                     new_chromosome,new_objf)
      implicit none
      real*8 crossover_probability,mutation_intensity
      integer length_chromosome,population_size
      real*8 genes_min(population_size),genes_max(population_size)
      real*8 old_chromosome(population_size,length_chromosome)
      real*8 old_objf(population_size)
      real*8 new_chromosome(population_size,length_chromosome)
      real*8 new_objf(population_size)
c    ================================================================
      integer minpp,maxpp,sister1,sister2,RANDOM_INTEGER,counter
      integer i,j,jcross
      real*8 OBJF,baby_objf
      real*8, allocatable :: baby_chromosome(:)
      logical FLIP

c     initialize the random number generator
      call RANDOMIZE()
      allocate(baby_chromosome(length_chromosome))
c     find the best individual in the old population. He is the father
c     for all children. He will have a child with each mother with the
c     help of two of his sisters. All the other individuals in the old
c     population are mothers. The child will compete with their mother
c     for a position in the next generation.
      call STATISTICS(population_size,old_objf,minpp,maxpp)
      do 110 i=1,population_size
c     elitist model. The best individual in the old population is
c     kept.
        if(i.eq.minpp) then
          do 10 j=1,length_chromosome
            new_chromosome(i,j)=old_chromosome(i,j)
10        end do
          new_objf(i)=old_objf(i)
          go to 110
        end if
20      sister1=RANDOM_INTEGER(population_size)
        if(sister1.eq.i) go to 20
30      sister2=RANDOM_INTEGER(population_size)
        if(sister2.eq.i.or.sister2.eq.sister1) go to 30
c     generate crossover mate (mutation)
        do 60 j=1,length_chromosome
c     for minization problem
          baby_chromosome(j)=old_chromosome(minpp,j)
     +                      +mutation_intensity
     +                       *(old_chromosome(sister2,j)
     +                         -old_chromosome(sister1,j))
c     for maximization problem
c          baby_chromosome(j)=old_chromosome(maxpp,j)
c     +                      +mutation_intensity
c     +                       *(old_chromosome(sister2,j)
c     +                         -old_chromosome(sister1,j))
c     to avoid the genes of the crossover mate exceed the upper limit
40        if(baby_chromosome(j).gt.genes_max(j)) then
            baby_chromosome(j)=(baby_chromosome(j)+genes_min(j))/2.0
            go to 40
          end if
c     to avoid the genes of the crossover mate exceed the lower limit
50        if(baby_chromosome(j).lt.genes_min(j)) then
            baby_chromosome(j)=(baby_chromosome(j)+genes_max(j))/2.0
            go to 50
          end if
60      end do
        baby_objf=OBJF(length_chromosome,baby_chromosome)
c     the baby competes with his mother
        if(baby_objf.lt.old_objf(i)) then ! for minimization problem
c        if(new_objf(i).gt.old_objf(i)) then ! for maximization problem
c     The baby wins
          do 70 j=1,length_chromosome
            new_chromosome(i,j)=baby_chromosome(j)
70        end do
          new_objf(i)=baby_objf
          go to 110
        end if
c     the mother wins, crossover to deliver another baby
        counter=0
        do 80 j=1,length_chromosome
          if(FLIP(crossover_probability)) then
            counter=counter+1
          else
            baby_chromosome(j)=old_chromosome(i,j)
          end if
80      end do
c     Full crossover (counter=population_size) is prohibited
        if(counter.eq.length_chromosome) then
          jcross=RANDOM_INTEGER(length_chromosome)
          baby_chromosome(jcross)=old_chromosome(i,jcross)
        end if
        baby_objf=OBJF(length_chromosome,baby_chromosome)
c     the baby competes with his mother
        if(baby_objf.lt.old_objf(i)) then ! for minimization problem
c        if(new_objf(i).gt.old_objf(i)) then ! for maximization problem
c     The baby wins
          do 90 j=1,length_chromosome
            new_chromosome(i,j)=baby_chromosome(j)
90        end do
          new_objf(i)=baby_objf
        else
c     The mother wins
          do 100 j=1,length_chromosome
            new_chromosome(i,j)=old_chromosome(i,j)
100       end do
          new_objf(i)=old_objf(i)
        end if
110   end do
      deallocate(baby_chromosome)
      return
      end
c    ****************************************************************
c     Benouli Experiment
c    ================================================================
      logical FUNCTION FLIP(probability)
      implicit none
      integer RANDOM_INTEGER
      real*8 probability,p
      p=RANDOM_INTEGER(20000)/20000.0
      FLIP=.false.
      if(p.le.probability) FLIP=.true.
      return
      end
c    ****************************************************************
c     initialize the real random generator
c    ================================================================
      SUBROUTINE RANDOMIZE()
      integer, static :: IFF
c     common /randomnumber/IFF,ix1,ix2,ix3,R
c     static IFF,ix1,ix2,ix3,R
      IFF=0
      return
      end
c    ****************************************************************
c     integer random number generator. It generates an integer random
c     number in [1,N]
c    ================================================================
      integer FUNCTION RANDOM_INTEGER(N)
      implicit none
      integer N
c    ================================================================
      real*8 RANDOM_REAL
      RANDOM_INTEGER=N*RANDOM_REAL()+1
      if(RANDOM_INTEGER.gt.N) RANDOM_INTEGER=N
      if(RANDOM_INTEGER.lt.1) RANDOM_INTEGER=1
      return
      end
c    ****************************************************************
c     real random number generator. It generates a real random number
c     in [0,1]. It also initialize the random generator if IFF=0
c    ================================================================
      real*8 FUNCTION RANDOM_REAL()
c     for FORTRAN Powerstation
      USE Portlib
c     for FORTRAN 90 on Planet Workstations at Temasek Laboratories,
c     program should be compiled as f90 -lV77
c      real*4 secnds
      integer SEED
      integer, static :: IFF,ix1,ix2,ix3
      real*8, static :: R(97)
      real*8 RM1,RM2
c     common /randomnumber/IFF,ix1,ix2,ix3,R
c     static IFF,ix1,ix2,ix3,R
      parameter(M1=259200,IA1=7141,IC1=54773,RM1=1.0/M1)
      parameter(M2=134456,IA2=8121,IC2=28411,RM2=1.0/M2)
      parameter(M3=243000,IA3=4561,IC3=51349)

      if(IFF.eq.0) then
        SEED=int(secnds(0.0)*100.0)
        SEED=mod(SEED,100000)
        IFF=1
        ix1=mod(IA1*SEEd+IC1,M1)
        ix2=mod(Ix1,M2)
        ix1=mod(IA1*ix1+IC1,M1)
        ix3=mod(Ix1,M3)
        do 10 j=1,97
          ix1=mod(IA1*ix1+IC1,M1)
          ix2=mod(IA2*ix2+IC2,M2)
          R(j)=(dble(Ix1)+dble(ix2)*RM2)*RM1
10      end do
      end if
20    ix1=mod(IA1*ix1+IC1,M1)
      ix2=mod(IA2*ix2+IC2,M2)
      ix3=mod(IA3*ix3+IC3,M3)
      j=1+(97*ix3)/M3
      if(j.gt.97.or.j.lt.1) go to 20
      RANDOM_REAL=R(j)
      R(j)=(dble(Ix1)+dble(ix2)*RM2)*RM1
      return
      end
c    ****************************************************************
c     statistics of the population
c    ================================================================
      SUBROUTINE STATISTICS(N,value,minpp,maxpp)
      implicit none
      integer N,minpp,maxpp
      real*8 value(N)
c    ================================================================
      integer i
      real*8 minvalue,maxvalue
      minpp=1
      maxpp=1
      minvalue=value(1)
      maxvalue=value(1)
      do 10 i=2,N
        if(value(i).lt.minvalue) then
          minpp=i
          minvalue=value(i)
        end if
        if(value(i).gt.maxvalue) then
          maxpp=i
          maxvalue=value(i)
        end if
10    end do
      return
      end
c    ****************************************************************
c     Objective function
c    ================================================================
      real*8 FUNCTION OBJF(N,x)
implicit none
c====defination sentences============================================

      integer N
real*8 x(N)
real*8, parameter :: pi=3.1415926
      real*8, parameter :: eps=1d-3
c   real*8 error,error1,error2
real*8,allocatable :: TTEDB(:)
real*8,allocatable :: freq(:)
integer Nfreq
      real*8 lowest_freq,freq_step
       real*8 p,d1,d2,w1,w2,g1,g2,Er,psi
c     ===============================================================
       real*8 lambda,TTETE
real*8 tmp,error1
       real*8 Fw,Fg
       real*8 X1p,X2p,B1p,B2p
       real*8 X1,X2,B1,B2
       real*8 Y,ft,fr1,fr2,error3
      integer i,count
N=5
      p=x(1)
w1=x(2)
w2=x(3)
g1=x(4)
g2=x(5)
Nfreq=181
lowest_freq=1.
freq_step=.1
ER=1.35
psi=0

c error=0.
error1=0.
c error2=0.
allocate(TTEDB(Nfreq))
allocate(freq(Nfreq))
d1=p-g1
d2=d1-2*w1-2*g2
      psi=2.*pi*psi/360 ! degree to radian

DO 11 i=1,Nfreq
        freq(i)=lowest_freq+freq_step*(i-1)
        lambda=300./freq(i) !light_speed=3.d8;the unit of lambda is mm
        X1p=Fw(p,w1,psi,lambda)                             ! (3.9a)
        X2p=Fw(p,w2,psi,lambda)                             ! (3.9b)
        X1=2*X1p*X2p/(X1p+X2p)*d1/p                         ! (3.8)
        X2=Fw(p,2*w2,psi,lambda)*d2/p                       ! (3.12)
        B1p=4*Er*Fg(p,g1,psi,lambda)                        ! (3.11)
        B2p=4*Er*Fg(p,g2,psi,lambda)                        ! (3.14)
        B1=0.75*B1p*d1/p                                    ! (3.10)
        B2=d2/p*B1p*B2p/(B1p+B2p)                           ! (3.13)
        Y=1/(X1-1/B1)+1/(X2-1/B2)       !
        TTETE=4./(4+Y*Y)  ! (2.12)
        IF(TTETE.LT.eps) TTETE=eps
          TTEDB(i)=10.*DLOG10(TTETE)
11 end do
      count=0
do i=2,Nfreq-1
     if(TTEDB(i).le.TTEDB(i-1).and.TTEDB(i).lt.TTEDB(i+1)) then
        if(count.eq.0) then
        fr1=freq(i)
           count=count+1
         else
          fr2=freq(i)
         end if
     end if
    if(TTEDB(i).ge.TTEDB(i-1).and.TTEDB(i).gt.TTEDB(i+1)) then
              ft=freq(i)
    end if
      end do
if (abs(fr1-9.0).lt.0.1.and.abs(fr2-15.0).lt.0.1)then
    error3=0
else
    error3=0.1
end if
do i=1,10
    if (TTEDB(i).ge.-0.5) then
          tmp=0
         else
          tmp=abs(TTEDB(i)+0.5)
         end if
     ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

共1条评分

tensor

rf币 +5

积极参与论坛交流，加分！

2010-01-05

发帖回复