Fortran能否从外部指定real变量kind值简化多精度欧拉法代码

Fortran属于静态编译型语言，变量类型必须在编译阶段确定，无法在运行时动态修改变量的kind参数，我们可以通过泛型接口封装重复逻辑的方式消除冗余代码：将重复的欧拉法计算、耗时统计、结果输出逻辑封装为对应不同kind的子例程，主程序仅保留kind选择的分支逻辑即可，无需重复定义多组变量、编写多份重复的运算代码。

优化后完整代码 #

module euler_mod
    implicit none
    ! 泛型接口定义，可统一调用不同精度的计算子例程
    generic :: euler_calc => euler_calc_k4, euler_calc_k8, euler_calc_k16
contains
    ! kind=4 精度的计算子例程
    subroutine euler_calc_k4()
        implicit none
        real(kind=4) :: t, h, y, s, e, r
        integer :: i
        t = 0.0_4
        y = 10.0_4
        print *, "enter the grid step"
        read(*,*) h
        r = 10.0_4 * exp(-5.0_4)
        call cpu_time(s)
        do i = 1, 999999999
            if (0.0_4 <= t .and. t <= 25.0_4) then
                y = y - h * y / 5.0_4
                t = t + h
            end if
        end do
        call cpu_time(e)
        print *, "solution is", y
        print *, "the error is", (r - y)/r
        print *, "time taken in seconds for kind 4 =", e - s
    end subroutine euler_calc_k4

    ! kind=8 精度的计算子例程
    subroutine euler_calc_k8()
        implicit none
        real(kind=8) :: t, h, y, s, e, r
        integer :: i
        t = 0.0_8
        y = 10.0_8
        print *, "enter the grid step"
        read(*,*) h
        r = 10.0_8 * exp(-5.0_8)
        call cpu_time(s)
        do i = 1, 999999999
            if (0.0_8 <= t .and. t <= 25.0_8) then
                y = y - h * y / 5.0_8
                t = t + h
            end if
        end do
        call cpu_time(e)
        print *, "solution is", y
        print *, "the error is", (r - y)/r
        print *, "time taken in seconds for kind 8 =", e - s
    end subroutine euler_calc_k8

    ! kind=16 精度的计算子例程
    subroutine euler_calc_k16()
        implicit none
        real(kind=16) :: t, h, y, s, e, r
        integer :: i
        t = 0.0_16
        y = 10.0_16
        print *, "enter the grid step"
        read(*,*) h
        r = 10.0_16 * exp(-5.0_16)
        call cpu_time(s)
        do i = 1, 999999999
            if (0.0_16 <= t .and. t <= 25.0_16) then
                y = y - h * y / 5.0_16
                t = t + h
            end if
        end do
        call cpu_time(e)
        print *, "solution is", y
        print *, "the error is", (r - y)/r
        print *, "time taken in seconds for kind 16 =", e - s
    end subroutine euler_calc_k16
end module euler_mod

program euler
    use euler_mod
    implicit none
    integer :: k
    print *, "enter the kind you want to work with"
    read(*,*) k
    select case(k)
        case(4)
            call euler_calc_k4()
        case(8)
            call euler_calc_k8()
        case(16)
            call euler_calc_k16()
        case default
            print *, "invalid kind value, only 4/8/16 are supported"
    end select
end program euler

进一步优化建议 #

如果使用的编译器支持Fortran 2018标准，可以使用template模板特性，仅编写一份模板子例程即可实例化出不同kind的版本，进一步降低代码冗余度。

内容的提问来源于stack exchange，提问作者user187604