*23456789012345678901234567890123456789012345678901234567890123456789012 function funmin(p,n) implicit double precision (a-h,o-z) DIMENSION P(N) save include 'para' * parameter (llw=5,llm=1,lls=1,nnd=1,nn=llw+lls+2+nnd*(llm+1)) * mm is the number of the inpout data points * parameter (mm=50) dimension polm(nnd),x(nnd),gfw(nn) dimension xset(mm,nnd),y(mm) data itric/0/ * if (itric.eq.0) then itric=1 if (n.ne.nn) then write(*,*)' SUBSUM: # of Variables N not equal to parameter nn' write(*,*)' N = ',N,' nn = ', nn stop endif * * Fetch Data: ******* USER PROVIDED SUBROUTINE: FETDAT call fetdat(xset,mm,nnd,y) ******* endif * funmin=0. do i=1,mm do j=1,nnd x(j)=xset(i,j) enddo call fwann(x,nnd,p,nn,llw,llm,lls,polw,polm,pols,fwnn,gfw) RES= fwnn-y(i) funmin=funmin+res**2 enddo call getfla(1,xf1) if (xf1.eq.1.) then call setfla(1,0.) itap = nunit() open(itap,file='res') rewind itap do i=1,mm do j=1,nnd x(j)=xset(i,j) enddo call fwann(x,nnd,p,nn,llw,llm,lls,polw,polm,pols,fwnn,gfw) write(itap,*) (x(j),j=1,nnd), fwnn, y(i) enddo close(itap) endif end *======================================================================= subroutine fetdat(xset,mm,nnd,y) implicit double precision (a-h,o-z) dimension xset(mm,nnd),y(mm) nt=nunit() open(nt,file='datal') rewind(nt) do m=1,mm read(nt,*)(xset(m,j),j=1,nnd) enddo close(nt) ntt=nunit() open(ntt,file='datat') rewind(ntt) do m=1,mm read(ntt,*)y(m) enddo close(ntt) end *======================================================================= subroutine granal(n,p,grad) implicit double precision (a-h,o-z) DIMENSION P(N), GRAD(N) save include 'para' * parameter (llw=5,llm=1,lls=1,nnd=1,nn=llw+lls+2+nnd*(llm+1)) * mm is the number of the inpout data points * parameter (mm=50) dimension polm(nnd),x(nnd),gfw(nn) dimension xset(mm,nnd),y(mm) data itrc/0/ * if (itrc.eq.0) then itrc=1 if(n.ne.nn) STOP 'GRANAL:' * * Fetch Data: USER PROVIDED SUBROUTINE: FETDAT ******* call fetdat(xset,mm,nnd,y) ******* endif * do ix = 1,n grad(ix)= 0.d0 enddo do i=1,mm do j=1,nnd x(j)=xset(i,j) enddo call fwann(x,nnd,p,n,llw,llm,lls,polw,polm,pols,fwnn,gfw) do ix = 1,n grad(ix)= grad(ix) + 2*(fwnn-y(i))*gfw(ix) enddo enddo end *======================================================================= subroutine polynom(p,n,lw,lm,ls,nd,s,polw,polm,pols) implicit double precision (a-h,o-z) save dimension p(n),polm(nd) * W(0:LW) coef: p(1:LW+1) * MU(0:LM,ND) coef: P(LW+2:LW+1+ND*(LM+1)) * SIG(i) coef: P(LW+2+ND*(LM+1):LW+1+ND*(LM+1)+LS+1) data itr/0/ if(itr.eq.0) then if(N.ne.LW+LS+2+ND*(LM+1)) then write(*,*)'Parameter count inconsistent' stop endif itr=1 endif * polw = p(lw+1) do i=1,lw polw=polw*s+p(lw+1-i) enddo * pols=p(n) do i=1,ls pols = pols*s+p(n-i) enddo * do j=1,nd polm(j) = p(lw+1+j*(lm+1)) do i=1,lm polm(j)=polm(j)*s+p(lw+1+j*(lm+1)-i) enddo enddo end *======================================================================= *23456789012345678901234567890123456789012345678901234567890123456789012 *======================================================================= subroutine fwann(x,nd,p,n,lw,lm,ls,polw,polm,pols,fwnn,gfw) implicit double precision (a-h,o-z) parameter (mi=100) save dimension x(nd),p(n),polm(nd),gfw(n),tse(mi) data itrr/0/ if (itrr.eq.0) then itrr=1 pi=acos(-1.d0) do i=1,mi tse(i) = cos((2.d0*i-1.d0)/(2.d0*mi)*pi) enddo endif * fwnn = 0. do i=1,n gfw(i)=0. enddo do mx=1,mi s=tse(mx) call polynom(p,n,lw,lm,ls,nd,s,polw,polm,pols) anor=0.d0 do i=1,nd anor=anor+(x(i)-polm(i))**2 enddo * expn=-anor/pols**2/2.d0 exfac = exp(expn) fs=polw*exfac fwnn = fwnn + fs do i=0,lw k=i+1 gfw(k)=gfw(k)+s**i*exfac enddo do j=1,nd do ia=0,lm k=lw+1+ia+(lm+1)*j-lm gfw(k)=gfw(k)+fs*s**ia*(x(j)-polm(j))/pols**2 enddo enddo do i = 0,ls k=lw+1+nd*(lm+1)+1+i gfw(k) = gfw(k)+fs*s**i*anor/pols**3 enddo enddo fwnn = fwnn*pi/mi do i=1,n gfw(i)=gfw(i)*pi/mi enddo end *======================================================================= *======================================================================= subroutine janal(m,n,p,fj,ld) implicit double precision (a-h,o-z) DIMENSION P(N),FJ(LD,n) save include 'para' * parameter (llw=5,llm=1,lls=1,nnd=1,nn=llw+lls+2+nnd*(llm+1)) * mm is the number of the inpout data points * parameter (mm=50) dimension polm(nnd),x(nnd),gfw(nn) dimension xset(mm,nnd),y(mm) data itrc/0/ * if (itrc.eq.0) then itrc=1 * * Fetch Data: USER PROVIDED SUBROUTINE: FETDAT ******* call fetdat(xset,mm,nnd,y) ******* endif * do i=1,mm do j=1,nnd x(j)=xset(i,j) enddo call fwann(x,nnd,p,n,llw,llm,lls,polw,polm,pols,fwnn,gfw) do ka = 1,n fj(i,ka)=gfw(ka) enddo enddo end