高斯-赛德尔迭代解线性方程组的问题

鹰击长空

两个高斯-赛德尔迭代解线性方程组子程序，算法一样，只是单精度和双精度的区别。麻烦帮忙看看是哪里的问题

单精度

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

    subroutine GS_solver (A,b,n,x)
        implicit none
        
        integer,parameter :: step = 5000  !最大迭代步数
        real   ,parameter :: err  = 1e-7 !误差
        
        integer :: i,j,k,n
        real    :: A(n,n),b(n),x0(n),x(n)
        real    :: x1(n)
        real    :: s,dx
        
        x0 = 0
        x1 = x0
        x  = x0

        
        do k = 1,step
            dx = 0
            
            do i = 1,n
               x1 = x0
               x1(i) = 0
               s = 0 
               s = dot_product( A(i,:), x1 )
               
               x(i) = ( b(i) - s ) / A(i,i) 
               
               dx = dx + ( x0(i)-x(i) )**2
               x0(i) = x(i) 
            end do
            
            dx = sqrt( dx )
            if ( dx < err ) exit
            
            
            x0 = x

    end do
    
 end subroutine GS_solver

双精度

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

    subroutine GSeidelSolver (A,b,n,x)
        implicit none
        
        integer   ,parameter :: MaxStep = 9000   !最大迭代步数
        real(8)   ,parameter :: cov  = 1d-8      !误差
        integer  :: i, j, k, n, step
        real(8)  :: A(n,n), b(n), x(n)
        real(8)  :: x0(n), x1(n)
        real(8)  :: dx, err
        
        x0 = 0d0
        x1 = 0d0
        x  = 0d0
        
        LOOP: do step = 1, MaxStep
            
            dx = 0d0
            do i = 1,n

                x1 = x0
                
                x1(i) = 0
                
                k = 0d0
                k = dot_product( A(i,:), x1 )
                
                x(i) = ( b(i) - k ) / A(i,i)
                dx = dx + ( x(i)-x0(i) )**2
                
                x0(i) = x(i) 
         
           end do
         
        err = dsqrt( dx )
                        
       if ( err < cov )  then
                exit LOOP
        end if

        x0 = x          

  end do LOOP         

 end subroutine GSeidelSolver