PROGRAM gauss_syoukyohou
    implicit none
    ! example
    double precision :: A(4,4) = reshape( (/ &
        10., 1., 4., 0., &
        1., 10., 5., -1., &
        4., 5., 10., 7., &
        0., -1., 7., 9. &
        /), shape(A) )
    double precision :: b(4) = (/ 15., 15., 26., 15. /)
    double precision :: x(4) = 0

    double precision temp_r
    integer i, j, k, dim, temp_i
    
    dim = 4

    do i = 1, dim
        temp_r = A(i,1)
        do j = 2, dim
            if (temp_r < A(i,j)) then
                temp_r = A(i,j)
            end if
        end do
        do j = 1, dim
            A(i,j) = A(i,j) / temp_r
        end do
        b(i) = b(i) / temp_r
    end do

    do i = 1, dim - 1
        temp_i = i
        do k = i, dim
            if (A(i,i) < A(k,i)) then
                temp_i = k
            end if
        end do
        do k = 1, dim
            temp_r = A(i,k)
            A(i,k) = A(temp_i,k)
            A(temp_i,k) = temp_r
        end do
        temp_r = b(i)
        b(i) = b(temp_i)
        b(temp_i) = temp_r    
        do k = i + 1, dim
            temp_r = A(k,i) / A(i,i)
            do j = i, dim
                A(k,j) = A(k,j) - A(i,j) * temp_r
            end do
            b(k) = b(k) - b(i) * temp_r
        end do
    end do
    
    do i = dim, 1, -1
        x(i) = b(i) / A(i,i)
        do j = 1, dim
            b(j) = b(j) - A(j,i) * x(i)
        end do
    end do

    print *, '################'
    print *, 'answer'
    do i = 1, dim
        print *, x(i)
    end do

END PROGRAM