最近用fortran编写最小二乘曲线模拟程序，一直出现问题 - Fortran Coder

[Fortran] 纯文本查看 复制代码
program Console1

implicit none

integer::i,j,n,status,m,f,k

real*8 Xmin,Ymin,px,px1,average,yy,dt1,dt2,dt3,z,d1,p,c,d2,g,q,dt,err

real*8,allocatable::X(:),Y(:),d(:),h(:),a(:)

real s(20),t(20),b(20)

open(10,file="E:\Fortran\task\2017.8.3\Console1\Console1\11.dat",status="old")

k=0

do while(.true.)

  read(10,*,iostat=status)

  if(status/=0)exit

  k=k+1

end do

n=k

allocate(X(n),Y(n))

rewind(10)

do j=1,n

  read(10,*)X(j),Y(j)

end do

Xmin=minval(X)

Ymin=minval(Y)

close(10)

if(Ymin<=0.0)then

  Y=log10(Y-Ymin+1.0)

else

    Y=log10(Y)

end if

if(Xmin<=0.0)then

  X=log10(X-Xmin+1.0)

else

  X=log10(X)

end if

read (*,*)m

px=(X(1)+X(n))/2.0

X=X-px

f=m-1

px1=((real(n)+1.0)/sum(x**((2.0)*f)))**((1.0)/(2.0*f))

x=px1*x

d=y

loop1:do

    call hpir1(x,y,a,n,m,dt1,dt2,dt3)

    write(9,*)(i,a(i),i=1,m)

    yy=0

    average=sum(x)/real(n)

    outter:do j=1,n

        h(0)=0

        inner:do i=1,m

            h(j)=h(j-1)+a(i)*((x(j)-average)**(i-1))

        end do inner

        yy=max(yy,abs(y(j)-h(j)))

    end do outter

    if(y(j)>h(j)+yy/real(2.0))then

            y(j)=h(j)+yy/real(2.0)

    elseif(y(j)<h(j)-yy/real(2.0))then

            y(j)=h(j)-yy/real(2.0)

    end if

    call RMS(i,y,d,n,err)

    if(err<1.0e-6.or.j==999)then

        exit loop1

    else

        j=j+1

    end if

end do loop1

open(9,file="result01")

write(9,400)j,m,n,err,dt1,yy

400 format(3(i3.3,2x),es12.4,2x,es12.4,2x,f12.4)

end program



subroutine RMS(i,y,d,n,err)

implicit none

integer,intent(in)::n

real*8,intent(in)::d(i),y(i)

real*8 err

integer::i

err=0

do i=1,n

    err=err+(d(i)-y(i))**2

end do

err=sqrt(err/real(n))

return 

end



subroutine hpir1(x,y,a,n,m,dt1,dt2,dt3)

real*8 x(n),y(n),a(m),s(20),t(20),b(20)

double precision dt1,dt2,dt3,z,d1,p,c,d2,g,q,dt

do 5 i=1,m

5 a(i)=0.0

if (m>n) m=n

if (m>20) m=20

z=0.0

do 10 i=1,n

10 z=z+x(i)/n

b(1)=1.0

d1=n

p=0.0

c=0.0

do 20 i=1,n 

    p=p+(x(i)-z)

    c=c+y(i)

20 continue

c=c/d1

p=p/d1

a(1)=c*b(1)

if (m>1) then

t(2)=1.0

t(1)=-p

d2=0.0

c=0.0

g=0.0

do 30 i=1,n

    q=x(i)-z-p

    d2=d2+q*q

    c=y(i)*q+c

    g=(x(i)-z)*q*q+g

30 continue

c=c/d2

p=g/d2

q=d2/d1

d1=d2

a(2)=c*t(2)

a(1)=c*t(1)+a(1)

endif

do 100 j=3,m

    s(j)=t(j-1)

    s(j-1)=-p*t(j-1)+t(j-2)

if (j>=4) then

do 40  k=j-2,2,-1

40 s(k)=-p*t(k)+t(k-1)-q*b(k)

end if

s(1)=-p*T(1)-q*b(1)

d2=0.0

c=0.0

g=0.0

do 70 i=1,n

    q=s(j)

    do 60 k=j-1,1,-1

60  q=q*(x(i)-z)+s(k)

    d2=d2+q*q

    c=y(i)*q+c

    g=(x(i)-z)*q*q+g

70 continue

c=c/d2

p=g/d2

q=d2/d1

d1=d2

a(j)=c*s(j)

t(j)=s(j)

do 80 k=j-1,1,-1

    a(k)=c*s(k)+a(k)

    b(k)=t(k)

    t(k)=s(k)

80 continue

100 continue

dt1=0.0

dt2=0.0

dt3=0.0

do 120 i=1,n

    q=a(m)

    do 110 k=m-1,1,-1

110 q=q*(x(i)-z)+a(k)

dt=q-y(i)

if (abs(dt)>dt3) dt3=abs(dt)

dt1=dt1+dt*dt

dt2=dt2+abs(dt)

120 continue

return

end