Why julia / transonic!

kuramoto_sivashinksy/codes/count_heads.jl

@@ -0,0 +1,7 @@
+function count_heads(n)
+    c::Int = 0
+    for i=1:n
+        c += rand(Bool)
+    end
+    c
+end

kuramoto_sivashinksy/codes/f.cpp

@@ -0,0 +1,4 @@
+double f(double x) {
+  return 4*x*(1-x);
+}
+

kuramoto_sivashinksy/codes/fN.cpp

@@ -0,0 +1,26 @@
+#include <stdlib.h>
+#include <iostream>
+#include <iomanip>
+#include <ctime>
+
+using namespace std;
+
+inline double f(double x) {
+  return 4*x*(1-x);
+}
+
+int main(int argc, char* argv[]) {
+  double x = argc > 1 ? atof(argv[1]) : 0.0;
+
+  double t0 = clock();
+  for (int n=0; n<1000000; ++n)
+    x = f(x);
+  double t1 = clock();
+
+  cout << "t = " << (t1-t0)/CLOCKS_PER_SEC << " seconds" << endl;
+  cout << setprecision(17);
+  cout << "x = " << x << endl;
+
+  return 0;
+}
+

kuramoto_sivashinksy/codes/ks-r2c.f90

@@ -0,0 +1,108 @@
+ module KS
+   implicit none
+   save
+   integer,          parameter :: Nx = 8192
+   double precision, parameter :: dt = 1d0/16d0
+   double precision, parameter :: T  = 200d0
+   double precision, parameter :: pi = 3.14159265358979323846d0
+   double precision, parameter :: Lx = (pi/16d0)*Nx
+ end module KS
+
+!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ program main
+   use KS
+   implicit none
+   integer :: i,Nt,r, Nruns=5, skip=1
+   double precision :: x(0:Nx-1), u(0:Nx-1)
+   real :: tstart, tend, avgtime
+
+   Nt = floor(T/dt)
+   x = Lx * (/(i,i=0,Nx-1)/) / dble(Nx)
+   u = cos(x) + 0.2d0*sin(x/8d0) + 0.01d0*cos(x/16d0)
+
+   avgtime = 0.
+   do r = 1, Nruns
+      call cpu_time(tstart)
+      call ksintegrate(Nt, u)
+      call cpu_time(tend)
+      if(r>skip)  avgtime = avgtime + tend-tstart
+   end do
+   avgtime = avgtime/(Nruns-skip)
+   print*, 'avgtime (seconds) = ', avgtime
+
+   print*, 'norm(u) = ', norm(u)
+
+   open(10, status='unknown', file='u.dat')
+   write(10,'(2e20.12)') (x(i),u(i), i=0,Nx-1)
+   close(10)
+
+ end program main
+
+!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+ double precision function norm(us,Nx)
+    double complex, intent(in) :: us(0:Nx-1)
+    norm = sqrt(sum(abs(us))/Nx)
+ end function norm
+
+!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ subroutine ksintegrate(Nt, u)
+   use KS
+   implicit none
+   integer,          intent(in)    :: Nt
+   double precision, intent(inout) :: u(0:Nx-1)
+   double precision :: dt2, dt32, Nx_inv
+   double precision :: A(0:Nx-1), B(0:Nx-1), L(0:Nx-1), kx(0:Nx-1)
+   double complex   :: alpha(0:Nx-1), G(0:Nx-1), Nn(0:Nx-1), Nn1(0:Nx-1)
+   double precision, save :: u_(0:Nx-1), uu(0:Nx-1) 
+   double complex,   save :: us(0:Nx-1), uus(0:Nx-1)
+   logical, save :: planned=.false.
+   integer*8 :: plan_u2us, plan_us2u, plan_uu2uus, plan_uus2uu
+   integer :: n, fftw_patient=32, fftw_estimate=64
+
+   if(.not.planned) then
+      call dfftw_plan_dft_r2c_1d(plan_u2us, Nx, u_, us, fftw_estimate)
+      call dfftw_plan_dft_c2r_1d(plan_us2u, Nx, us, u_, fftw_estimate)
+      call dfftw_plan_dft_r2c_1d(plan_uu2uus, Nx, uu, uus, fftw_estimate)
+      call dfftw_plan_dft_c2r_1d(plan_uus2uu, Nx, uus, uu, fftw_estimate)
+      planned = .true.
+   end if
+
+   do n = 0, Nx/2-1
+      kx(n) = n
+   end do
+   kx(Nx/2) = 0d0
+   do n = Nx/2+1, Nx-1
+      kx(n) = -Nx + n;
+   end do
+   alpha = (2d0*pi/Lx)*kx
+   L = alpha**2 - alpha**4
+   G = -0.5d0*dcmplx(0d0,1d0)*alpha
+
+   Nx_inv = 1d0/Nx
+   dt2  = dt/2d0
+   dt32 = 3d0*dt/2d0
+   A = 1d0 + dt2*L
+   B = 1d0/(1d0 - dt2*L)
+
+   uu = u*u
+   call dfftw_execute(plan_uu2uus)
+   Nn  = G*uus
+   Nn1 = Nn
+
+   u_ = u
+   call dfftw_execute(plan_u2us)
+
+   do n = 1, Nt
+      Nn1 = Nn
+      call dfftw_execute(plan_us2u)
+      u_ = u_ * Nx_inv
+      uu = u_*u_
+      call dfftw_execute(plan_uu2uus)
+      Nn = G*uus
+      us = B * (A*us + dt32*Nn - dt2*Nn1)
+   end do
+
+   call dfftw_execute(plan_us2u)
+   u = u_ * Nx_inv
+
+ end subroutine ksintegrate

kuramoto_sivashinksy/codes/ksbenchmark.c

@@ -0,0 +1,174 @@
+#include <math.h>
+#include <time.h>
+#include <complex.h>
+#include <fftw3.h>
+#include <stdlib.h>
+#include <stdio.h>
+//#include <stdbool.h>
+
+typedef complex double Complex;
+
+void ksintegrate(Complex* u, int Nx, double Lx, double dt, int Nt);
+double ksnorm(Complex* u, int Nx);
+
+const double pi = 3.14159265358979323846; /* why doesn't M_PI from math.j work? :-( */
+
+int main(int argc, char* argv[]) {
+
+  /* to be set as command-line args */
+  int Nx = 0;         /* number of x gridpoints */
+  int Nruns = 5;      /* number of benchmarking runs */
+  int printnorms = 0; /* can't figure out C bool type :-( */
+
+  if (argc < 2) {
+    printf("please provide one integer argument Nx\n");
+    exit(1);
+  }
+  Nx = atoi(argv[1]);  
+
+  if (argc >= 3)
+    Nruns = atoi(argv[2]);
+
+  if (argc >= 4)
+    printnorms = atoi(argv[3]);
+
+  printf("Nx == %d\n", Nx);
+
+  double dt = 1.0/16.0;
+  double T  = 200.0;
+  double Lx = (pi/16.0)*Nx;
+  double dx = Lx/Nx;
+  int Nt = (int)(T/dt);
+
+  double* x   = malloc(Nx*sizeof(double));
+  Complex* u0 = malloc(Nx*sizeof(Complex));
+  Complex* u  = malloc(Nx*sizeof(Complex));
+
+  for (int n=0; n<Nx; ++n) {
+    x[n] = n*dx;
+    u0[n] = cos(x[n]) + 0.1*sin(x[n]/8) + 0.01*cos((2*pi/Lx)*x[n]);
+  }
+
+  int skip = 1;
+  double avgtime = 0.0;
+
+  for (int run=0; run<Nruns; ++run) {
+
+    for (int n=0; n<Nx; ++n)
+      u[n] = u0[n];
+
+    clock_t tic = clock();
+    ksintegrate(u, Nx, Lx, dt, Nt);
+    clock_t toc = clock();
+
+    double cputime =((double)(toc - tic))/CLOCKS_PER_SEC;
+    if (run >= skip) 
+      avgtime += cputime;
+
+    printf("cputtime == %f\n", cputime);
+  }
+  printf("norm(u(0)) == %f\n", ksnorm(u0, Nx)); 
+  printf("norm(u(T)) == %f\n", ksnorm(u, Nx)); 
+
+
+  avgtime /= (Nruns-skip);
+  printf("avgtime == %f\n", avgtime);
+
+  free(u0);
+  free(u);
+  free(x);
+}
+
+void ksintegrate(Complex* u, int Nx, double Lx, double dt, int Nt) {
+  int* kx        = malloc(Nx*sizeof(int));
+  double*  alpha = malloc(Nx*sizeof(double));  
+  Complex* D     = malloc(Nx*sizeof(Complex)); 
+  Complex* G     = malloc(Nx*sizeof(Complex)); 
+  double*  L     = malloc(Nx*sizeof(double));  
+
+  for (int n=0; n<Nx/2; ++n)      
+    kx[n] = n;
+  kx[Nx/2] = 0.0;                 
+  for (int n=Nx/2+1; n<Nx; ++n)   
+    kx[n] = -Nx + n;
+
+  const double a = 2*pi/Lx;
+  for (int n=0; n<Nx; ++n) {   
+    alpha[n] = a*kx[n];   
+    double alpha2 = alpha[n]*alpha[n];  /* alpha[n]^2 */
+    D[n] = alpha[n]*I;                     
+    L[n] = alpha2*(1-alpha2);
+    G[n] = -0.5*alpha[n]*I;
+  }
+  double dt2  = dt/2;
+  double dt32 = 3*dt/2;
+  double Nx_inv = 1.0/Nx;
+  double* A = malloc(Nx*sizeof(double));
+  double* B = malloc(Nx*sizeof(double));
+  for (int n=0; n<Nx; ++n) {
+    A[n] = 1.0 + dt2*L[n];
+    B[n] = 1.0/(1.0 - dt2*L[n]);
+  }
+  Complex* uu = malloc(Nx*sizeof(Complex));  
+  for (int n=0; n<Nx; ++n) 
+    uu[n] = u[n]*u[n];
+
+  fftw_complex* u_fftw  =  u;
+  fftw_complex* uu_fftw =  uu;
+  fftw_plan u_fftw_plan   = fftw_plan_dft_1d(Nx, u_fftw,  u_fftw,  FFTW_FORWARD, FFTW_ESTIMATE);
+  fftw_plan u_ifftw_plan  = fftw_plan_dft_1d(Nx, u_fftw,  u_fftw,  FFTW_BACKWARD, FFTW_ESTIMATE);
+  fftw_plan uu_fftw_plan  = fftw_plan_dft_1d(Nx, uu_fftw, uu_fftw, FFTW_FORWARD,  FFTW_ESTIMATE);
+  fftw_plan uu_ifftw_plan = fftw_plan_dft_1d(Nx, uu_fftw, uu_fftw, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+  fftw_execute(uu_fftw_plan); 
+  Complex* Nn = malloc(Nx*sizeof(Complex));
+  for (int n=0; n<Nx; ++n) 
+    Nn[n] = G[n]*uu[n];
+  Complex* Nn1 = malloc(Nx*sizeof(Complex));
+
+  fftw_execute(u_fftw_plan);
+
+  for (int s=0; s<Nt; ++s) {
+    for (int n=0; n<Nx; ++n) {
+      Nn1[n] = Nn[n];
+      uu[n] = u[n];
+    }      
+    fftw_execute(uu_ifftw_plan);
+    for (int n=0; n<Nx; ++n) 
+      uu[n] *= Nx_inv; 
+    for (int n=0; n<Nx; ++n) 
+      uu[n] = uu[n]*uu[n];       /* compute uu = u^2 physical */
+    fftw_execute(uu_fftw_plan);  /* transform uu physical -> spectral */
+    for (int n=0; n<Nx; ++n) 
+      Nn[n] = G[n]*uu[n];        /* compute Nn == -u u_x, spectral */
+    for (int n=0; n<Nx; ++n) 
+      u[n] = B[n] * (A[n] * u[n] + dt32*Nn[n] - dt2*Nn1[n]);
+  }
+  fftw_execute(u_ifftw_plan);
+  for (int n=0; n<Nx; ++n) 
+      u[n] *= Nx_inv; 
+
+  fftw_destroy_plan(u_fftw_plan);
+  fftw_destroy_plan(uu_fftw_plan);
+  fftw_destroy_plan(u_ifftw_plan);
+  fftw_destroy_plan(uu_ifftw_plan);
+  free(kx);
+  free(alpha);
+  free(D);
+  free(G);
+  free(L);
+  free(A);
+  free(B);
+  free(uu);
+  free(Nn);
+  free(Nn1);
+}
+
+inline double square(double x) {return x*x;}
+
+double ksnorm(Complex* u, int Nx) {
+  double s = 0.0;
+  for (int n=0; n<Nx; ++n)
+    s += square(cabs(u[n]));
+  return sqrt(s/Nx);
+}