Неожиданное поведение с основным гауссовым устранением в Fortran

Question

Я пытаюсь решить систему линейных уравнений в матричной форме, используя одну из многих широко доступных подпрограмм, в частности эту функцию gauss (наряду с образцовыми данными с известным результатом), но независимо от того, какую стратегию устранения я использую, я постоянно получаю неправильные результаты и не знаю, куда обратиться в этот момент. Мой код:

MODULE gaussMod 
    CONTAINS 
    function gauss(a,b) result(x) 
     implicit none 
     real*8 :: b(:), a(size(b), size(b)) 
     real*8 :: x(size(b)) 

     real*8 :: r(size(b)) 
     integer i,j, neq 

     neq = size(b) 

     do i =1, neq 
      r = a(:,i)/a(i,i) 

      do j = i+1, neq 
       a(j,:) = a(j,:) - r(j)*a(i,:) 
       b(j) = b(j) - r(j)*b(j) 
      enddo 
     enddo 

     do i= neq, 1, -1 
      x(i) = (b(i) - sum(a(i, i+1:) * x(i+1:)))/a(i,i) 
     enddo 
    END function 
END MODULE 

SUBROUTINE outFile(n,x) 
    IMPLICIT none 
    INTEGER n 
    INTEGER i 
    REAL*8, DIMENSION(n) :: x 


    OPEN(UNIT=20,FILE="solution.csv",action="write",status="replace") 

    DO i=1, n 
     WRITE (20,"(1(f0.30,',',:))") x(i) 
    END DO 
    CLOSE(20) 
END SUBROUTINE 

PROGRAM elimtest 
     USE gaussMod 
     IMPLICIT none 
     INTEGER, PARAMETER :: coeff_kind = selected_real_kind(p=30, r=99) 
     INTEGER n 
     REAL*8, DIMENSION(:,:), ALLOCATABLE :: a 
     REAL*8, DIMENSION(:), ALLOCATABLE :: b 
     REAL*8, DIMENSION(:), ALLOCATABLE :: x 
     n = 4 

     ALLOCATE(a(n,n)) 
     ALLOCATE(b(n)) 
     ALLOCATE(x(n)) 

      a(1,1) = 18. 
      a(1,2) = -6. 
      a(1,3) = -6. 
      a(1,4) = 0. 
      a(2,1) = -6. 
      a(2,2) = 12. 
      a(2,3) = 0. 
      a(2,4) = -6. 
      a(3,1) = -6. 
      a(3,2) = 0. 
      a(3,3) = 12. 
      a(3,4) = -6. 
      a(4,1) = 0. 
      a(4,2) = -6. 
      a(4,3) = -6. 
      a(4,4) = 18. 

      b(1) = 60. 
      b(2) = 0. 
      b(3) = 20. 
      b(4) = 0. 


     x = gauss(a,b) 
     CALL outFile(n,x) 


END PROGRAM 

Любая помощь была бы принята с благодарностью!

Answer 1

Похоже, в следующей строке

b(j) = b(j) - r(j)*b(j) 

есть опечатка

b(j) = b(j) - r(j)*b(i) 

С этой модификацией ваш код дает правильный результат (возможно!):

x(:) = 8.3333333333333339 6.6666666666666670 8.3333333333333339 5.0000000000000000