solve.f90

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

program solve

    ! this tells FORTRAN compiler that gauss_elim
    ! is defined elsewhere, and linker will find it
    external gauss_elim

    integer :: n = 4
    integer i,j

    double precision A(4,4)
    double precision L(4,4)
    double precision U(4,4)

    double precision b(4)
    double precision x(4)
    double precision y(4)

    double precision :: h=1.256637061435917

    do i=1,n
        do j=1,n
            A(i,j)=0
        end do
        A(i,i)=-2/(h*h)
        if ( i .gt. 1 ) then
            A(i,i-1) = 1/(h*h)
            A(i-1,i) = 1/(h*h)
        end if
    end do

    b(1)=0.951056516295154
    b(2)=0.587785252292473
    b(3)=-b(2)
    b(4)=-b(1)

    ! note that there is no underscore at the end
    ! of the function name, unlike the C++ code.
    ! also, the C++ function expects a pointer for 
    ! the last argument, but here it is passed
    ! normally, and by reference
    call gauss_elim( A, L, U, n )

    call printmatrix( 'A', A, n )
    call printmatrix( 'U', U, n )
    call printmatrix( 'L', L, n )
    
    do i=1,n
        y(i)=b(i)
        do j=1,i-1
            y(i)=y(i)-L(i,j)*y(j)
        end do
    end do

    do i=n,1,-1
        x(i)=y(i)
        do j=i+1,n
            x(i)=x(i)-U(i,j)*x(j)
        end do
        x(i)=x(i)/U(i,i)
    end do

    write (*,*) 'x = '
    do i=1,n
        write (*,'(f10.5)') x(i)
    end do

end

subroutine printmatrix( name, A, n )

    double precision A(n,n)
    character name
    integer i, n
    write (*,*) name, ' = '
    do i=1,n
        write (*,'(4f10.5)') A(i,1:n)
    end do

end