fortran 中是否有任何方法可以从外部定义实数 4 或实数 8？

Question

我们知道我们可以从外部优化我们的代码。
我们知道在 fortran 编程中，我们首先在程序中定义变量。然后我们可以从外部（通过读取语句）获取输入，但我们可以编写一个从外部获取某种变量的小代码。

即如果我们在终端中放入 4 种 4（即真实（种类 = 4））变量，即如果我们放入 8 种 8（即真实（种类 = 8）变量，则引入。有什么办法吗。
我知道我们可以做三个单独的 if 循环来定义变量（即 kind 4 、 kind 8 、 kind 16 并重复该程序三次）。
我写的代码是用欧拉法求 y 的值。 我想推广到任何类型并计算所花费的时间。我希望这可以以不那么繁琐的方式完成。

我写的代码：

       program euler
       implicit none
       real(kind=4)::t,h,y,s,e,r
       real(kind=8)::t8,h8,y8,s8,e8,r8
       real(kind=16)::t16,h16,y16,s16,e16,r16
       integer::k,i
       t=0
       t8=0
       t16=0
       y=10
       y8=10
       y16=10
       print *,"enter the kind you want to work with"
       read(*,*) k
       !so if user writes 4 kind 4 variables would do the work
       if(k==4) then
       print *,"enter the grid step"
       read(*,*) h
       r=10*exp(-5.0)
       call cpu_time(s)
       do i=1,999999999
       if(0.le.t.and.t.le.25) then
       y=y-h*y/5.0
       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 =",e-s
      else if(k==8) then
      print *,"enter the grid step"
      read(*,*) h8
     r8=10*exp(-5.0D0)
      call cpu_time(s8)
      do i=1,999999999
      if(0.le.t8.and.t8.le.25) then
      y8=y8-h8*y8/5.0
      t8=t8+h8
      end if
      end do
      call cpu_time(e8)
      print *,"solution is",y8
      print *,"the error is",(r8-y8)/r8
      print *,"time taken in seconds for kind 8 =",e8-s8
      else if(k==16) then
      print *,"enter the grid step"
      read(*,*) h16
      r16=10*exp(-5.0D0)
      call cpu_time(s16)
      do i=1,999999999
      if(0.le.t16.and.t16.le.25) then
      y16=y16-h16*y16/5.0
      t16=t16+h16
      end if
      end do
      call cpu_time(e16)
      print *,"solution is",y16
      print *,"the error is",(r16-y16)/r16
      print *,"time taken in seconds for kind 16 =",e16-s16
      end if
      end program euler
 

但我看起来更聪明，更不麻烦。

Answer 1

我认为这不是你想要的 100%，因为仍然有三个单独的子例程，每个子例程都用于不同的类型，但它们是使用 interface euler 调用的，它根据参数确定使用哪一个

program SO_Euler
use iso_fortran_env, only : sp=>real32, dp=>real64, qp=>real128, i4=>int32, i8=>int64
implicit none

interface euler
    procedure euler_driver_sp, euler_driver_dp, euler_driver_qp
end interface

real(sp), parameter :: r4 = 10*exp(-5.0)
real(dp), parameter :: r8 = 10*exp(-5d0)
real(qp), parameter :: r16 = 10*exp(-5q0)

real(sp) :: h
print *,"enter the grid step"
read(*,*) h    
print *, ""

call euler(r4, h)
call euler(r8, h)
call euler(r16, h)

contains 

subroutine euler_driver_sp(r,h_in)
real(sp), intent(in) :: r
real(sp), intent(in) :: h_in
real(sp) :: h, y, t
integer(i8) :: s, e, rate
integer :: i
    print '(a15,1x,g0)', "kind is ", kind(r)
    h = h_in
    t = 0
    y = 10
    call SYSTEM_CLOCK(s,rate)
    do i=1,999999999
        if(0<=t .and.  t<=25) then
            y=y-h*y/5
            t=t+h
        else
            exit
        end if            
    end do
    call SYSTEM_CLOCK(e,rate)
    print '(a15,1x,g0.15)',"solution is", y 
    print '(a15,1x,g0.15)',"the error is", (r-y)/r
    print '(a15,1x,g0.4,1x,a)',"time taken is", real(e-s)/rate,"seconds"
    print *, ""
end subroutine

subroutine euler_driver_dp(r, h_in)
real(dp), intent(in) :: r
real(sp), intent(in) :: h_in
real(dp) :: h, y, t
integer(i8) :: s, e, rate
integer :: i
    print '(a15,1x,g0)', "kind is ", kind(r)
    h = h_in    !! convert sp=>dp
    t = 0
    y = 10
    call SYSTEM_CLOCK(s,rate)
    do i=1,999999999
        if(0<=t .and.  t<=25) then
            y=y-h*y/5
            t=t+h
        else
            exit
        end if            
    end do
    call SYSTEM_CLOCK(e,rate)
    print '(a15,1x,g0.15)',"solution is", y 
    print '(a15,1x,g0.15)',"the error is", (r-y)/r
    print '(a15,1x,g0.4,1x,a)',"time taken is", real(e-s)/rate,"seconds"
    print *, ""
end subroutine

subroutine euler_driver_qp(r, h_in)
real(qp), intent(in) :: r
real(sp), intent(in) :: h_in
real(qp) :: h, y, t
integer(i8) :: s, e, rate
integer :: i
    print '(a15,1x,g0)', "kind is ", kind(r)
    h = h_in ! convert sp=>qp
    t = 0
    y = 10
    call SYSTEM_CLOCK(s,rate)
    do i=1,999999999
        if(0<=t .and.  t<=25) then
            y=y-h*y/5
            t=t+h
        else
            exit
        end if            
    end do
    call SYSTEM_CLOCK(e,rate)
    print '(a15,1x,g0.15)',"solution is", y 
    print '(a15,1x,g0.15)',"the error is", (r-y)/r
    print '(a15,1x,g0.4,1x,a)',"time taken is", real(e-s)/rate,"seconds"
    print *, ""
end subroutine


end program

这是该过程的一些示例输出

 enter the grid step
0.000002

       kind is  4
    solution is .547848604619503E-01
   the error is .186920538544655
  time taken is .1020 seconds

       kind is  8
    solution is .673793765102040E-01
   the error is .138737586862949E-05
  time taken is .7200E-01 seconds

       kind is  16
    solution is .673793765102226E-01
   the error is .138737559174033E-05
  time taken is 1.535 seconds

请注意，我在 64 位发布模式下编译，浮点模型 不是 快，但 严格 以及扩展实常数精度的选项.