2

I'm currently working on a library of temporal integrators. I'm wondering what's best for computing efficiency: Computing the coefficients for every time step, or using them as parameters from an external repository. The first case would be this way:

Subroutine getRKF45(this, dt, y, z)
  Implicit none
  Class(diffEc), Intent(InOut) :: this
  Real(8), Intent(In) :: dt
  Real(8), dimension(:), Intent(InOut) :: y, z
  Real(8), dimension(6,this%n) :: k
  Call field(this%deriv, this%x, this%t)
  k(1,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*(1.d0/5.d0), this%t+dt*(1.d0/5.d0))
  k(2,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*(3.d0/40.d0)+k(2,:)*dt*(9.d0/40.d0) &
     , this%t + dt*(3.d0/10.d0))
  k(3,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*(3.d0/10.d0)-k(2,:)*dt*(9.d0/10.d0) &
     +k(3,:)*dt*(6.d0/5.d0), this%t + dt*(3.d0/5.d0))
  k(4,:) = this%deriv(:)
  Call field(this%deriv, this%x-k(1,:)*dt*(11.d0/54.d0)+k(2,:)*dt*(5.d0/2.d0) &
     -k(3,:)*dt*(70.d0/27.d0)+k(4,:)*dt*(35.d0/27.d0), this%t + dt)
  k(5,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*(1631.d0/55296.d0)+k(2,:)*dt*(175.d0/512.d0) &
     +k(3,:)*dt*(575.d0/13824.d0)+k(4,:)*dt*(44275.d0/110592.d0) &
     +k(5,:)*dt*(253.d0/4096.d0), this%t + dt*(7.d0/8.d0))
  k(6,:) = this%deriv(:)
  y(:) = this%x(:) + dt*(k(1,:)*(37.d0/378.d0)+k(3,:)*(250.d0/621.d0)+k(4,:)*(125.d0/594.d0) &
     +k(6,:)*(512.d0/1771.d0))
  z(:) = this%x(:) + dt*(k(1,:)*(2825.d0/27648.d0)+k(3,:)*(18575.d0/48384.d0) &
     +k(4,:)*(13525.d0/55296.d0)+k(5,:)*(277.d0/14336.d0)+k(6,:)*(1.d0/4.d0))
End Subroutine getRKF45

And using coefficients pre-computed:

Subroutine getRKF45(this, dt, y, z)
  Use IntegratorUtilities, only: rk45Coef
  Implicit none
  Class(diffEc), Intent(InOut) :: this
  Real(8), Intent(In) :: dt
  Real(8), dimension(:), Intent(InOut) :: y, z
  Real(8), dimension(6,this%n) :: k
  Call field(this%deriv, this%x, this%t)
  k(1,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(1), this%t+dt*rk45Coef(1))
  k(2,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(2)+k(2,:)*dt*rk45Coef(3) &
     , this%t + dt*rk45Coef(4))
  k(3,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(4)-k(2,:)*dt*rk45Coef(5) &
     +k(3,:)*dt*rk45Coef(6), this%t + dt*rk45Coef(7))
  k(4,:) = this%deriv(:)
  Call field(this%deriv, this%x-k(1,:)*dt*rk45Coef(8)+k(2,:)*dt*rk45Coef(9) &
     -k(3,:)*dt*rk45Coef(10)+k(4,:)*dt*rk45Coef(11), this%t + dt)
  k(5,:) = this%deriv(:)
  Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(12)+k(2,:)*dt*rk45Coef(13) &
     +k(3,:)*dt*rk45Coef(14)+k(4,:)*dt*rk45Coef(15) &
     +k(5,:)*dt*rk45Coef(16), this%t + dt*rk45Coef(17))
  k(6,:) = this%deriv(:)
  y(:) = this%x(:) + dt*(k(1,:)*rk45Coef(18)+k(3,:)*rk45Coef(19)+k(4,:)*rk45Coef(20) &
     +k(6,:)*rk45Coef(21))
  z(:) = this%x(:) + dt*(k(1,:)*rk45Coef(22)+k(3,:)*rk45Coef(23) &
     +k(4,:)*rk45Coef(24)+k(5,:)*rk45Coef(25)+k(6,:)*rk45Coef(26))
End Subroutine getRKF45

I've run some test and the time difference between the two schemes is not relevant. Does the compiler pre-process these coefficients multiplications?

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
Facundo Nicolas Airaudo
  • 21
  • 1
  • Welcome, remember to take the Welcome [tour]. Maybe just the part of code is not relevant in general. Did you test whether those lines are really important for the running time? Maybe it spends more time in `call field(..)` anyway. – Vladimir F Героям слава Jun 05 '18 at 14:34
  • We have to know how is `rk45Coef` defined. Is it a constant array or a variable array? – Vladimir F Героям слава Jun 05 '18 at 14:34
  • BTW using `real(8)` is ugly and not portable and does not have to be compatible with the double precision you use for your constants like `1.d0`. More at https://stackoverflow.com/questions/838310/fortran-90-kind-parameter – Vladimir F Героям слава Jun 05 '18 at 14:38
  • 3
    gfortran will constant-fold the numerical constants. For example, (1.d0/4.d0) will be computed at compile to be 0.25d0. You can inspect an intermediate representation of your code with -fdump-tree-original. – Steve Jun 05 '18 at 18:07

0 Answers0