Fortran 中的“Varargin”

计算科学 正则
2021-12-19 00:59:25

考虑以下问题：我想在 Fortran 中求解一个微分方程（使用 Runge-Kutta 例程），例如x'=-a*x，其中a是一个常数。然后我想为不同的a. 我现在想编写代码，以便将其应用于其他问题。但总的来说，我不知道我的微分方程可能需要什么额外的输入（以及它们的大小和类型）。所以我的问题是，是否可以在 Fortran 中将所有附加参数传递给我的 Runge-Kutta 方法（最好只有一个属性），以便将其传递给我的微分方程，以便在那里使用？有了一组已知的参数，这很容易，但是未知的集合呢？

到目前为止，我的方法是这样的：

PROGRAM main
  USE RKModule

  IMPLICIT NONE

  INTEGER :: l

  REAL*8, PARAMETER :: num_t = 1e3
  REAL*8, PARAMETER :: min_t = 0.0d0
  REAL*8, PARAMETER :: max_t = 10.0d0

  REAL*8, DIMENSION(1), PARAMETER :: a = 1.8d0

  REAL*8, DIMENSION(num_t) :: dGrid_t, grid_t
  COMPLEX*16, DIMENSION(num_t) :: sol
  COMPLEX*16, DIMENSION(1) :: x_in, x_out

  DO l = 1, num_t
    dGrid_t(l) = (max_t - min_t) / DBLE(num_t)
    grid_t(l) = min_t + DBLE(l-0.5d0)*dGrid_t(l)
  END DO

  x_in(1) = dcmplx(1.0d0, 0.0d0)
  x_out(1) = dcmplx(0.0d0, 0.0d0)

  DO l = 1, num_t
    sol(l) = x_in(1)
    CALL rk45(x_out, DGL, x_in, grid_t(l), dGrid_t(l), a)
    x_in = x_out
  END DO

  OPEN(UNIT = 1000, FILE = 'sol.dat', ACTION = 'WRITE')
  WRITE(1000, FMT = '(e14.7e3, 1x, e14.7e3)', ADVANCE = "NO") sol
  CLOSE(1000)

  CONTAINS

  SUBROUTINE DGL(t,f,df,vararg_in)
    REAL*8, INTENT(IN) :: t
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in

    REAL*8, DIMENSION(1) :: b

    b = 1.0d0
    IF (PRESENT(vararg_in)) THEN
      b = vararg_in
    END IF

    df = -b*f

  END SUBROUTINE DGL

END PROGRAM main

还有 Runge Kutta 套路：

MODULE RKModule
  IMPLICIT NONE

  CONTAINS

  SUBROUTINE rk45(x_out, RHS, x_in, t, dt, vararg_in)
    INTERFACE
      SUBROUTINE RHS(t,f,df,a)
    REAL*8, INTENT(IN) :: t
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: a
      END SUBROUTINE
    END INTERFACE
    REAL*8, INTENT(IN) :: t, dt 
    COMPLEX*16, DIMENSION(:), INTENT(IN) :: x_in
    COMPLEX*16, DIMENSION(:), INTENT(OUT) :: x_out
    REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in

    REAL*8 :: tHlp
    COMPLEX*16, DIMENSION(SIZE(x_in)) :: k1 ,k2, k3, k4, dx, x_hlp
    REAL*8, DIMENSION(:), ALLOCATABLE :: vararg_in_copy

    IF (PRESENT(vararg_in)) THEN
      ALLOCATE(vararg_in_copy(SIZE(vararg_in)))
      vararg_in_copy = vararg_in

      tHlp = t
      x_hlp = x_in
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k1 = dt * dx

      tHlp = t + 0.5d0 * dt
      x_hlp = x_in + 0.5d0 * k1
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k2 = dt * dx

      x_hlp = x_in + 0.5d0 * k2
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k3 = dt * dx

      tHlp = t + dt
      x_hlp = x_in + k3
      CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
      k4 = dt * dx

      x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0

      DEALLOCATE(vararg_in_copy)
    ELSE

      tHlp = t
      x_hlp = x_in
      CALL RHS(tHlp, x_hlp, dx)
      k1 = dt * dx

      tHlp = t + 0.5d0 * dt
      x_hlp = x_in + 0.5d0 * k1
      CALL RHS(tHlp, x_hlp, dx)
      k2 = dt * dx

      x_hlp = x_in + 0.5d0 * k2
      CALL RHS(tHlp, x_hlp, dx)
      k3 = dt * dx

      tHlp = t + dt
      x_hlp = x_in + k3
      CALL RHS(tHlp, x_hlp, dx)
      k4 = dt * dx

      x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0

    END IF

  END SUBROUTINE rk45

END MODULE RKModule
1个回答

这个问题更适合 StackOverflow。

您可以使用class(*)或type(c_ptr)。两者都来自 Fortran 2003。后者得到广泛支持，前者仅从 gfortran 4.8 和其他编译器的其他最新版本开始。Notetype(c_ptr)即使在 Fortran 内部也非常有用，无需调用任何 C。

第一种方法：

subroutine RK4(normal arguments, vararg_in)
  class(*), intent(in) :: vararg_in(:)  !for scalar delete (:)

  !if you need to preserve a copy (why?)
  allocate(vararg_copy, source=varrarg_in)

  call RHS(...,vararg_copy) !or directly vararg_in
end subroutine


subroutine DGL(...,vararg_in)
  class(*), intent(in) :: vararg_in(:)  !for scalar delete (:)

  select type (vararg_in)
    select type (expected_type)
     !do something
    class default
     !some exception
  end select
end subroutine

如果您需要不同类型的变量，您可以将任何内容传递给class(*)，可能是派生类型。

另一种方法：

subroutine RK4(normal arguments, vararg_in)
  use iso_c_binding
  type(c_ptr), intent(in) :: vararg_in_ptr

  !Doing a copy of the pointer doesn't make much sense

  call RHS(...,vararg_in_ptr)
end subroutine


subroutine DGL(...,vararg_in_ptr)
  use iso_c_binding
  type(c_ptr), intent(in) :: vararg_in_ptr
  type(expected_type), pointer :: vararg_in(:) !for scalar delete (:)

  !for an array
  call c_f_pointer(vararg_in_ptr, vararg_in, expected_shape)
  !for a scalar just
  call c_f_pointer(vararg_in_ptr, vararg_in)

  !do something
end subroutine
其它你可能感兴趣的问题
上一篇如何减少对称矩阵的计算时间？ 下一篇均衡各种多边形随机镶嵌面积的算法