gauss_elim.f90

斯坦福Energy211/CME211课《c++编程——地球科学科学家和工程师》的课件

subroutine gauss_elim(A, L, U, n)

    integer n
    double precision A(n,n)
    double precision L(n,n)
    double precision U(n,n)
    double precision m
	
    integer i, j, k

    do i=1,n
        do j=1,n
            L(i,j)=0
            U(i,j)=A(i,j)
        end do
        L(i,i)=1
    end do
	
    ! matrices in C++ are stored in row order,
    ! but FORTRAN uses column order, so here
    ! we work with the transpose
    do j = 1,n
        do i = j+1,n
            m = U(j,i) / U(j,j)
            do k = j,n
                U(k,i) = U(k,i) - m * U(k,j)
            end do
            L(j,i) = m
        end do	
    end do	
end