运行时出现 Program Exception - access violation，不知道什么原因？ - Fortran Coder

[Fortran] 纯文本查看 复制代码
module gauss

  !高斯消元模块：主要为牛顿法中解线性方程所调用

  contains

  subroutine lineq(A,b,x,N)

    implicit none

    integer(kind=2) :: id_max,i,k,N!主元素标号

    real(kind=8) :: A(N,N),b(N),x(N)

    real(kind=8) :: Aup(N,N),bup(N)

    real(kind=8) :: Ab(N,N+1)  !增广矩阵

    real(kind=8) :: vtemp1(N+1),vtemp2(N+1)

    real(kind=8) :: elmax,temp

    Ab(1:N,1:N) = A

    Ab(:,N+1) = b

    do k=1,N-1

      elmax = dabs(Ab(k,k))

      id_max = k

      do i=k+1,N

        if(abs(Ab(i,k))>elmax) then

          elmax = Ab(i,k)

          id_max = i

        endif

      enddo

      vtemp1 = Ab(k,:)

      vtemp2 = Ab(id_max,:)

      Ab(k,:) = vtemp2

      Ab(id_max,:) = vtemp1

      do i=k+1,N

        temp = Ab(i,k)/Ab(k,k)

        Ab(i,:) = Ab(i,:) - temp*Ab(k,:)

      enddo

    enddo

    Aup(:,:) = Ab(1:N,1:N)

    bup(:) = Ab(:,N+1)

    call uptri(Aup,bup,x,N)

    return

  end subroutine lineq

  

  !上三角线性方程组

  subroutine uptri(A,b,x,N)

    implicit none

    integer(kind=2) :: i,j,N

    real(kind=8) :: A(N,N),b(N),x(N)

    x(N) = b(N)/A(N,N)

    do i=n-1,1,-1

      x(i) = b(i)

      do j=i+1,N

        x(i) = x(i) - A(i,j)*x(j)

      enddo

      x(i) = x(i)/A(i,i)

    enddo

    return 

   end subroutine uptri

end module gauss

  

module newton

  contains

  subroutine solve

    use gauss

    implicit  none

    integer(kind=2) :: i,itmax = 500!最大允许迭代次数

    integer(kind=2) :: p,q,N=2

    real(kind=8) :: x(2),f(2),dx(2)

    real(kind=8) :: x0(2),l

    real(kind=8),allocatable :: s(:),t(:),s0(:),t0(:)

    real(kind=8) :: df(2,2) !偏导数矩阵

    real(kind=8) :: tol,dx2

    q=(p+1)*(p+2)/2

    allocate(s(q))

    allocate(t(q))

    allocate(s0(q))

    allocate(t0(q))

    !随机给定初值

    call RANDOM_SEED()

    call RANDOM_NUMBER(x)

    tol = 1d-4  !计算精度控制

    open(unit=11,file = 'parm.txt')

    read(11,*) p,x0,l

    write(*,*) p,x0,l

    do i=1,q

      read(11,100) s0(i),t0(i) 

      write(*,100) s0(i),t0(i) 

      pause '=============================='

    enddo

    !x = x0 

    s = s0/100

    t = t0/100

100 format(E15.8E2,1x,E15.8E2)

    open(unit=11,file = 'result.txt')

    write(11,101)

101 format(/,t6,'牛顿法计算非线性方程组迭代序列',/) 

    do i=1,itmax 

      call func(f,p,x,s,t,x0,l)

      call jac(df,p,x,s,t,x0,l)

      call lineq(df,-f,dx,N)

      x = x + dx

      write(*,102) i , x*l+x0

102   format(i5,2f16.10)      

      dx2 = dsqrt(dx(1)**2+dx(2)**2)

      if (dx2 < tol) exit

    enddo

    return

  end subroutine solve

  

  subroutine func(f,p,x,s,t,x0,l)

    implicit none

    real(kind=8) :: x(2),f(2)

    real(kind=8) :: x0(2),l

    integer(kind=2) :: p

    integer :: i,j,k=0

    real(kind=8) :: s((p+1)*(p+2)/2),t((p+1)*(p+2)/2)

    f = 0

    !f(1) = s(1)+x0(1)-10+(s(2)+l)*x(1)+s(3)*x(2)+s(4)*x(1)**2+s(5)*x(1)*x(2)+s(6)*x(2)**2+&

    !&s(7)*x(1)**3+s(8)*x(1)**2*x(2)+s(9)*x(1)*x(2)**2+s(10)*x(2)**3!+&

    !&s(11)*x(1)**4+s(12)*x(1)**3*x(2)+s(13)*x(1)**2*x(2)**2+s(14)*x(1)*x(2)**3+s(15)*x(2)**4

    !f(2) = t(1)+x0(2)-4.8+t(2)*x(1)+(t(3)+l)*x(2)+t(4)*x(1)**2+t(5)*x(1)*x(2)+t(6)*x(2)**2+&

    !&t(7)*x(1)**3+t(8)*x(1)**2*x(2)+t(9)*x(1)*x(2)**2+t(10)*x(2)**3!+&

    !&t(11)*x(1)**4+t(12)*x(1)**3*x(2)+t(13)*x(1)**2*x(2)**2+t(14)*x(1)*x(2)**3+t(15)*x(2)**4

    !s(2)=s(2)+l

    !t(3)=t(3)+l

    !s(1)=s(1)+x0(1)-10

    !t(1)=t(1)+x0(1)-4.8

    do i=0,p

      do j=0,i

        k = k+1

        f(1)=f(1)+s(k)*x(1)**(i-j)*x(2)**j

      enddo

    enddo

    k=0    

    do i=0,p

      do j=0,i

        k = k+1

        f(2)=f(2)+t(k)*x(1)**(i-j)*x(2)**j

      enddo

    enddo

  end subroutine func

  subroutine jac(df,p,x,s,t,x0,l)

    implicit none

    real(kind=8) :: x(2),df(2,2)

    integer(kind=2) :: p,i,j,k,m

    real(kind=8) :: s((p+1)*(p+2)/2),t((p+1)*(p+2)/2)

    real(kind=8) :: x0(2),l

    !df(1,1) = s(2)+l+2*s(4)*x(1)+s(5)*x(2)+3*s(7)*x(1)**2+2*s(8)*x(1)*x(2)+s(9)*x(2)**2!+&

              !&4*s(11)*x(1)**3+3*s(12)*x(1)**2*x(2)+2*s(13)*x(1)*x(2)**2+s(14)*x(2)**3

    !df(2,1) = s(3)+s(5)*x(1)+2*s(6)*x(2)+s(8)*x(1)**2+2*s(9)*x(1)*x(2)+3*s(10)*x(2)**3!+&

              !&s(12)*x(1)**3+2*s(13)*x(1)**2*x(2)+3*s(14)*x(1)*x(2)**2+4*s(15)*x(2)**3

    !df(1,2) = t(2)+2*t(4)*x(1)+t(5)*x(2)+3*t(7)*x(1)**2+2*t(8)*x(1)*x(2)+t(9)*x(2)**2!+&

              !&4*t(11)*x(1)**3+3*t(12)*x(1)**2*x(2)+2*t(13)*x(1)*x(2)**2+t(14)*x(2)**3

    !df(2,2) = t(3)+l+t(5)*x(1)+2*t(6)*x(2)+t(8)*x(1)**2+2*t(9)*x(1)*x(2)+3*t(10)*x(2)**3!+&

              !&s(12)*x(1)**3+2*s(13)*x(1)**2*x(2)+3*s(14)*x(1)*x(2)**2+4*s(15)*x(2)**3

    !s(2)=s(2)+l

    !t(3)=t(3)+l           

    df=0         

    !df(1,1)

    m=0

    do i=0,p-1

      k = i+1

      m = m+1

      do j=0,i

        m=m+1

        df(1,1)=df(1,1)+k*s(m)*x(1)**(i-j)*x(2)*j

        k=k-1 

      enddo

    enddo

    !df(1,2)

    m=0

    do i=0,p-1

      k = i+1

      m = m+1

      do j=0,i

        m=m+1

        df(1,2)=df(1,2)+k*t(m)*x(1)**(i-j)*x(2)*j

        k=k-1

      enddo

    enddo

    !df(2,1)

    m=1

    do i=0,p-1

      k = 1

      m = m+1

      do j=0,i

        m=m+1        

        df(2,1)=df(2,1)+k*s(m)*x(1)**(i-j)*x(2)*j

        k=k+1

      enddo

    enddo

    !df(2,2)

    m=1

    do i=0,p-1

      k = 1

      m = m+1

      do j=0,i

        m=m+1        

        df(2,2)=df(2,2)+k*t(m)*x(1)**(i-j)*x(2)*j

        k=k+1

      enddo

    enddo

  end subroutine jac

end module newton

  

  

program main

  use newton

  call solve

end program main