Numerik, Parallelrechnen und Hochleistung

program numerik
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: n = 1000
  real(dp) :: h, summe, x, pi_naehrung
  integer :: i

  ! Trapezregel: Integral von sin(x) von 0 bis pi
  h = acos(-1.0_dp) / n
  summe = 0.5_dp * (sin(0.0_dp) + sin(n * h))
  do i = 1, n - 1
    summe = summe + sin(i * h)
  end do
  print "(a, f10.6)", "Integral sin(x) = ", summe * h

  ! Leibniz-Reihe: pi/4 = 1 - 1/3 + 1/5 - ...
  summe = 0.0_dp
  do i = 0, 999999
    summe = summe + (-1.0_dp)**i / (2 * i + 1)
  end do
  pi_naehrung = 4.0_dp * summe
  print "(a, f12.8)", "pi ~ ", pi_naehrung

  ! Newton-Verfahren: sqrt(2)
  x = 1.0_dp
  do i = 1, 6
    x = x - (x * x - 2.0_dp) / (2.0_dp * x)
  end do
  print "(a, f12.10, a, es9.2)", "sqrt(2) ~ ", x, "   Fehler ", abs(x - sqrt(2.0_dp))
end program numerik

Integral sin(x) =   1.999998
pi ~   3.14159165
sqrt(2) ~ 1.4142135624   Fehler  2.22E-16

program waerme
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: n = 11, schritte = 200
  real(dp) :: u(n), neu(n), alpha
  integer :: t

  alpha = 0.25_dp                   ! stabil für alpha <= 0.5
  u = 0.0_dp
  u(1) = 100.0_dp                    ! linkes Ende heiß
  u(n) = 0.0_dp                      ! rechtes Ende kalt

  do t = 1, schritte
    neu = u
    neu(2:n-1) = u(2:n-1) + alpha * (u(1:n-2) - 2.0_dp * u(2:n-1) + u(3:n))
    u = neu
  end do
  print "(11f6.1)", u                 ! nahezu lineares Temperaturprofil
end program waerme

100.0  89.9  79.7  69.6  59.6  49.6  39.6  29.6  19.7   9.9   0.0

program parallel
  implicit none
  integer, parameter :: n = 1000000
  real, allocatable :: a(:), b(:), c(:)
  integer :: i

  allocate(a(n), b(n), c(n))
  a = 1.5
  b = 2.5

  do concurrent (i = 1:n)               ! Iterationen sind voneinander unabhängig
    c(i) = a(i) * b(i) + 1.0
  end do
  print *, c(1), c(n), sum(c) / n

  !$omp parallel do
  do i = 1, n
    c(i) = sqrt(a(i) + b(i))
  end do
  !$omp end parallel do
  print *, c(1)
end program parallel

4.75000000       4.75000000       4.77924681
   2.00000000

! Coarray-Beispiel (gfortran mit OpenCoarrays: caf / cafrun)
program caf
  implicit none
  integer :: summe[*]
  summe = this_image() * 10
  sync all
  if (this_image() == 1) print *, "Summe:", sum([(summe[i], i = 1, num_images())])
end program caf

program messen
  implicit none
  integer, parameter :: n = 300
  real :: a(n, n), b(n, n)
  real :: t0, t1
  integer :: i, j

  call random_number(a)
  call cpu_time(t0)
  do j = 1, n
    do i = 1, n                ! innerer Index zuerst: cachefreundlich
      b(i, j) = a(i, j) * 2.0
    end do
  end do
  call cpu_time(t1)
  print *, "Ergebnis stimmt:", all(abs(b - 2.0 * a) < 1.0e-6), " Zeit nicht negativ:", t1 >= t0
end program messen

Ergebnis stimmt: T  Zeit nicht negativ: T