前言：研究生的时候，利用测井曲线分析水淹层，导师让用rbf神经网络进行预测，为了搞懂神经网络，遍学各种神经网络算法，这是我学习中用fortran编写的BP神经网络程序，参考了很多经典的程序。也许对大家有用。

program Bp_Neutral_Net 

character*1 cc

character*15 wfilename

integer m_samplenum

common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum

common/nn/n2

common/ee/ee,es1,es

cc**********************************************************

double precision x(300),x1(300),y(300),dy

double precision w_ih(300,300),w_ho(300,300)

double precision sample_in(20000,300),sample_t(20000,300),

     &sample_out(20000,300),sample_x1(20000,300)

common/weight/w_ih,w_ho

common/sample/sample_in,sample_t,sample_out,sample_x1

integer simu_num

cc ee=1.0/4000.0

ee=1.0/8000.0

cc es1=0.01

es1=1.0/8000.0

cc es=0.01

      es=1.0/8000.0

write(*,*)'Do you want to learn(l) or calculate(c) or derive(d)?'

read(*,*)cc

if(cc.eq.'l'.or.cc.eq.'L') then

call bp_ini

call readsampledata

call Learning

write(*,*)'Enter the name of bp file!'

cc read(*,*)wfilename

open(2,file='bp.dat')

write(2,*)n2

write(2,*)m_inputnum,m_hidenum,m_outnum

write(2,*)((w_ih(i,j),j=1,m_inputnum+1),i=1,m_hidenum)

write(2,*)((w_ho(i,j),j=1,m_hidenum+1),i=1,m_outnum)

close(2)

elseif(cc.eq.'c'.or.cc.eq.'C') then

open(2,file='bp.dat')

read(2,*)n2

read(2,*)m_inputnum,m_hidenum,m_outnum

read(2,*)((w_ih(i,j),j=1,m_inputnum+1),i=1,m_hidenum)

read(2,*)((w_ho(i,j),j=1,m_hidenum+1),i=1,m_outnum)

close(2)

cc write(*,*)((w_ih(i,j),j=1,m_inputnum+1),i=1,m_hidenum)

write(*,*)'Enter the name of simu file!'

read(*,*)wfilename

open(5,file=wfilename)

read(5,*)simu_num

write(*,*) simu_num

cc write(*,*)'enter x:'

      do 500 j=1,simu_num

read(5,*)(x(i),i=1,m_inputnum)

call calculate_output(x,x1,y)

open(6,file='output.txt')

write(6,*)(y(i),i=1,m_outnum)

500   continue

      close(6)

close(5)

elseif(cc.eq.'d'.or.cc.eq.'D') then

open(2,file='bp.dat')

read(2,*)m_inputnum,m_hidenum,m_outnum

read(2,*)((w_ih(i,j),j=1,m_inputnum+1),i=1,m_hidenum)

read(2,*)((w_ho(i,j),j=1,m_hidenum+1),i=1,m_outnum)

close(2)

read(*,*)(x(i),i=1,m_inputnum)

read(*,*)i,j

call derivevalue(x,i,j,dy)

endif

end

subroutine calculate_output(x,x1,y)

double precision x(300),x1(300),y(300)

common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum

common/ee/ee,es1,es

common/nn/n2

do 230 i=1,m_inputnum

x(i)=x(i)*ee

230 continue

write(*,*)ee,es1,es,n2

call outputgenerate(x,x1,y)

write(*,*)(y(i),i=1,m_outnum)

if(n2.eq.0) then

do 235 i=1,m_outnum

y(i)=y(i)/es1

235 continue

elseif(n2.gt.0) then

do 240 i=1,m_outnum

y(i)=y(i)/es-2.0

240 continue

endif

end

subroutine bp_ini

double precision w_ih(300,300),w_ho(300,300)

double precision rn

common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum

common/weight/w_ih,w_ho

cc5 write(*,*)'The element number of input layer:'

cc read(*,*)m_inputnum

cc write(*,*)'The element number of hide layer:'

cc read(*,*)m_hidenum

cc write(*,*)'The element number of output layer:'

cc read(*,*)m_outnum

cc if(m_inputnum.gt.30.or.m_hidenum.gt.30.or.m_outnum.gt.30) then

cc write(*,*)'The net is too large!'

cc goto 5

cc endif

      m_inputnum=3

m_hidenum=3

m_outnum=1

r1=10

r2=100

do 20 i=1,m_hidenum

do 20 j=1,m_inputnum+1

20 w_ih(i,j)=rn(r1)

do 40 i=1,m_outnum

do 40 j=1,m_hidenum+1

40 w_ho(i,j)=rn(r2)

end

function rn(R)

      double precision rn

s=65536

u=2053

v=13849

r=u*r+v

m=r/s

r=r-m*s

rn=r/s

return 

end

subroutine readsampledata

character*15 filename

common/nn/n2

common/ee/ee,es1,es

common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum

double precision sample_in(20000,300),sample_t(20000,300),

     &sample_out(20000,300),sample_x1(20000,300)

common/sample/sample_in,sample_t,sample_out,sample_x1

cc     write(*,*)m_inputnum,m_hidenum,m_outnum

write(*,*)'Enter the name of sample file:'

read(*,*)filename

open(1,file=filename)

read(1,*)m_samplenum

CC WRITE(*,*)M_SAMPLENUM

do 50 i=1,m_samplenum

read(1,*)sample_in(i,1),sample_in(i,2),sample_in(i,3),

     &sample_out(i,1)

cc read(1,*)(sample_in(i,j),j=1,m_inputnum)

cc read(1,*)(sample_out(i,j),j=1,m_outnum)

cc WRITE(*,*)(sample_in(i,j),j=1,m_inputnum)

cc WRITE(*,*)(sample_out(i,j),j=1,m_outnum)

50 continue

n2=0

do 55 i=1,m_samplenum

do 54 j=1,m_outnum

if(sample_out(i,j).lt.0) n2=n2+1

54 continue

55 continue

do 49 i=1,m_samplenum 

do 56 j=1,m_inputnum

sample_in(i,j)=sample_in(i,j)*ee

56 continue

          if(n2.eq.0) then

do 57 j=1,m_outnum

sample_out(i,j)=sample_out(i,j)*es1

57 continue

elseif(n2.gt.0) then

do 53 j=1,m_outnum

sample_out(i,j)=(2+sample_out(i,j))*es

53 continue

endif

49 continue

close(1)

cc WRITE(*,*)(sample_in(i,j),j=1,m_inputnum)

cc WRITE(*,*)(sample_out(i,j),j=1,m_outnum)

end

subroutine outputgenerate(x,x1,y)

double precision x(300),x1(300),y(300),s,sig

double precision w_ih(300,300),w_ho(300,300)

common/bp/m_inputnum,m_hidenum,m_outnum

common/weight/w_ih,w_ho

x1(m_hidenum+1)=-1.

x(m_inputnum+1)=-1.

do 80 i=1,m_hidenum

s=0.

do 90 j=1,m_inputnum+1

s=s+w_ih(i,j)*x(j)

90 continue

x1(i)=sig(s)

80 continue

do 100 i=1,m_outnum

s=0.

do 110 j=1,m_hidenum+1

s=s+w_ho(i,j)*x1(j)

110 continue

y(i)=sig(s)

100 continue

end

subroutine Learning

double precision w_ih(300,300),w_ho(300,300)

double precision sample_in(20000,300),sample_t(20000,300),

     &sample_out(20000,300),sample_x1(20000,300)

common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum

common/weight/w_ih,w_ho

common/sample/sample_in,sample_t,sample_out,sample_x1

integer training_num

double precision  x(300),x1(300),y(300)

double precision w1(300,300),w2(300,300),w11(300,300),w22(300,300)

double precision err,err_whole,err_whole1

write(*,*)'Enter the error:'

read(*,*)err_xy

training_num=0

do 105 i=1,m_outnum

do 105 j=1,m_hidenum+1

105 w1(i,j)=w_ho(i,j)

do 110 i=1,m_hidenum

do 110 j=1,m_inputnum+1

110 w2(i,j)=w_ih(i,j)

err_whole1=0.

120 training_num=training_num+1

err_whole=0.

do 130 i=1,m_samplenum

do 135 k=1,m_inputnum

x(k)=sample_in(i,k)

135 continue

call outputgenerate(x,x1,y)

do 136 k=1,m_hidenum+1

sample_x1(i,k)=x1(k)

136 continue

do 137 k=1,m_outnum

sample_t(i,k)=y(k)

137 continue

err=0.

do 140 j=1,m_outnum

err=err+(sample_t(i,j)-sample_out(i,j))**2

140 continue

err=err*0.5

err_whole=err_whole+err

130 continue

uu=0.5

cc write(*,*)'err_whole=',err_whole,'err_whole1=',err_whole1

if(err_whole-err_whole1.lt.0) then

uu=1.20*uu

alfa=0.9

else

uu=0.60*uu

alfa=0.

endif

do 150 j=1,m_outnum

do 160 k=1,m_hidenum+1

s=0.

do 170 i=1,m_samplenum

s=s+(sample_t(i,j)-sample_out(i,j))*(1-sample_t(i,j))

     &*sample_t(i,j)*sample_x1(i,k)

170 continue

w11(j,k)=w_ho(j,k)-uu*s+alfa*(w_ho(j,k)-w1(j,k))

w1(j,k)=w_ho(j,k)

w_ho(j,k)=w11(j,k)

160 continue

150 continue

do 180 j=1,m_hidenum

do 190 k=1,m_inputnum+1

        ss=0.

do 200 i=1,m_samplenum

s=0.

do 210 l=1,m_outnum

s=s+(sample_t(i,l)-sample_out(i,l))*(1-sample_t(i,l))*

     &sample_t(i,l)*w_ho(l,j)

210 continue

ss=ss+s*(1-sample_x1(i,j))*sample_x1(i,j)*sample_in(i,k)

200 continue

w22(j,k)=w_ih(j,k)-uu*ss+alfa*(w_ih(j,k)-w2(j,k))

w2(j,k)=w_ih(j,k)

w_ih(j,k)=w22(j,k)

190 continue

180 continue

cc write(*,*)((w_ho(i,j),j=1,hide_num+1),i=1,out_num)

cc write(*,*)'err_whole=',err_whole

err_whole1=err_whole

if(training_num.gt.1000000) goto 111

if(err_whole.gt.err_xy) goto 120

111 write(*,*)'Training is over!'

cc write(*,*)'err_whole=',err_whole

write(*,*)'The training number is',training_num

end

subroutine derivevalue(x,iy,jx,dy)

double precision x(300),dy

integer iy,jx

double precision y(300),x1(300)

double precision w_ih(300,300),w_ho(300,300)

common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum

common/weight/w_ih,w_ho

call outputgenerate(x,x1,y)

s=0.

do 220 k=1,m_hidenum+1

s=s+(1-y(iy))*y(iy)*w_ho(iy,k)*(1-x1(k))*x1(k)*w_ih(k,jx)

220 continue      

dy=s

end  

function sig_d(x)

double precision sig_d,sig

double precision x

sig_d=sig(x)*(1.-sig(x))

end

function sig(x)

double precision sig

double precision x

sig=1./(1.+exp(-x))

end