Как компилятор gfortran обрабатывает числовые коэффициенты в выражениях? - PullRequest
0 голосов
/ 05 июня 2018

В настоящее время я работаю над библиотекой временных интеграторов.Мне интересно, что лучше для вычислительной эффективности: вычисление коэффициентов для каждого временного шага или использование их в качестве параметров из внешнего репозитория.Первый случай будет таким:

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

И с использованием предварительно вычисленных коэффициентов:

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

Я провел некоторый тест, и разница во времени между двумя схемами не имеет значения,Компилятор предварительно обрабатывает эти коэффициенты умножения?

...