主题:[讨论]Fortran数组越界问题
求助:有个关于最小二乘法拟合趋势的程序,编译没错,但运行是会出现数组越界问题,请教下高手是哪里的问题?
! CC*CHANG 为拟合趋势程序,其中AA-是输入变量,151为资料长度,
! CC 3,4表示拟合的阶数,3=2阶+1,4等于拟合阶数+2
! CC X, AT, ST均为工作数组
! CC 调过此子程序后,AA的时间序列已去除了趋势
PROGRAM MAIN
PARAMETER(N=73,m=36*13)
DIMENSION X(2,N),AT(3),ST(3,4),AA(N),vwnd(m,n),vwnd2(m,n)
!读入数据
open(1,file='vwnd.1996.txt',form='formatted')
do it=1,n
read(1,*)(vwnd(i,it),i=1,m)
enddo
close(1)
write(*,*)'input done'
!线性拟合
do i=1,m
do it=1,n
aa(it)=vwnd(i,it)
! write(*,*) aa(it),it
enddo
CALL CHANG(AA,n,3,4,X,AT,ST)
do j=1,n
vwnd2(i,j)=aa(j) !拟合后赋值输出
enddo
enddo
write(*,*) 'call subpoutine done'
open(1,file='vwndaa.1996.txt',form='formatted')
write(1,*)((vwnd2(i,it),i=1,m),it=1,n)
close(1)
write(*,*) 'output done'
! open(2,file='vwndaa.grd',form='formatted')
! do i=1,m
! write(2,rec=1)(vwnd2(i,it),it=1,n)
! enddo
! close(2)
end
SUBROUTINE CHANG(UU,N,M1,M2,X,A,S)
DIMENSION X(2,N),A(M1),S(M1,M2),UU(N)
DO 100 NT=1,N
100 X(1,NT)=FLOAT(NT)
DO 300 NT=1,N
300 X(2,NT)=UU(Nt)
M=M1-1
DO 5 I=1,M1
DO 5 J=1,M2
5 S(I,J)=0.0
S(1,1)=N
DO 6 J=1,N
6 S(1,M2)=S(1,M2)+X(2,J)
DO 7 I=2,M1
I1=I-1
MI=M1+I-2
DO 7 J=1,N
S(I,1)=S(I,1)+X(1,J)**I1
S(M1,I)=S(M1,I)+X(1,J)**MI
7 S(I,M2)=S(I,M2)+X(1,J)**I1*X(2,J)
DO 8 J=2,M
DO 8 I=J,M
8 S(I,J)=S(I+1,J-1)
DO 9 I=1,M
I1=I+1
DO 9 J=I1,M1
9 S(I,J)=S(J,I)
CALL GS(M1,S,A,M2,FAIL)
DO 500 NT=1,N
500 UU(NT)=UU(NT)-A(1)-A(2)*FLOAT(NT)-A(3)*FLOAT(NT*NT)
200 CONTINUE
RETURN
END
SUBROUTINE GS(N,A,X,N1,FAIL)
DIMENSION A(N,N1),X(N)
DO 11 K=1,N
DO 6 I=K,N
IF(A(I,K).EQ.0.0) GO TO 6
I0=I
GO TO 7
6 CONTINUE
FAIL=-1.0
RETURN
7 IF(I0.EQ.K) GO TO 9
DO 8 J=K,N1
T=A(K,J)
A(K,J)=A(I0,J)
8 A(I0,J)=T
9 K1=K+1
DO 10 J=K1,N1
10 A(K,J)=A(K,J)/A(K,K)
DO 11 I=K1,N
DO 11 J=K1,N1
11 A(I,J)=A(I,J)-A(I,K)*A(K,J)
X(N)=A(N,N1)
NI=N-1
DO 12 K=1,N1
K1=NI-K+1
X(K1)=A(K1,N1)
K2=K1+1
DO 12 J=K2,N
12 X(K1)=X(K1)-A(K1,J)*X(J)
FAIL=0.0
RETURN
END