ex2: fix benchmark

- work around bugs in savebmp by changing YSIZE
- optimize calculate by changing YSIZE
- write to total buffer in reverse order for correct bmp

.gitignore

@@ -54,4 +54,5 @@ dkms.conf
 
 *.zip
 *.pdf
-*.png
+*.bmp
+!main*

exercise0/Makefile

@@ -1,13 +1,17 @@
-build:
-	gcc main.c bitmap.c -o main -O3
+build: main.c bitmap.c
+	gcc main.c bitmap.c -o main -O2
 
-run:
-	make build
+run: build
 	./main
 
-show:
-	make run
+show: run
 	feh after.bmp
 
-convert:
+convert: run
 	magick after.bmp after.png # for showing image in pdf
+
+zip: before.bmp main.pdf Makefile main.c bitmap.h bitmap.c
+	zip handin.zip before.bmp main.pdf Makefile main.c bitmap.h bitmap.c
+
+unzip: handin.zip
+	unzip handin.zip -d handin

exercise1/.gitignore

@@ -0,0 +1,3 @@
+data
+images
+*.mp4

exercise1/Makefile

@@ -0,0 +1,19 @@
+CFLAGS+= -std=c99 -O2 -Wall -Wextra
+LDLIBS+= -lm
+TARGETS=wave_1d
+IMAGES=$(shell find data -type f | sed s/\\.dat/.png/g | sed s/data/images/g )
+.PHONY: all clean dirs plot movie
+all: dirs ${TARGETS}
+dirs:
+	mkdir -p data images
+plot: ${IMAGES}
+images/%.png: data/%.dat
+	./plot_image.sh $<
+movie: ${IMAGES}
+	ffmpeg -y -an -i images/%5d.png -vcodec libx264 -pix_fmt yuv420p -profile:v baseline -level 3 -r 12 wave.mp4
+clean:
+	-rm -fr ${TARGETS} data images wave.mp4
+zip: theory.pdf Makefile wave_1d.c
+	zip handin.zip theory.pdf Makefile wave_1d.c
+unzip: handin.zip
+	unzip handin.zip -d handin

exercise1/README.md

@@ -0,0 +1,4 @@
+* make       : builds the executable 'wave_1d'
+* ./wave\_1d : fills the 'data/' directory with stored time steps
+* make plot  : converts saved time steps to png files under 'images/', using gnuplot. Runs faster if launched with e.g. 4 threads (make -j4 plot).
+* make movie : converts collection of png files under 'images' into an mp4 movie file

exercise1/plot_image.sh

@@ -0,0 +1,10 @@
+#! /usr/bin/env bash
+SIZE=1024
+DATAFILE=$1
+IMAGEFILE=`echo $1 | sed s/dat$/png/ | sed s/data/images/`
+cat <<END_OF_SCRIPT | gnuplot -
+set term png
+set output "$IMAGEFILE"
+set yrange[-1:1]
+plot "$DATAFILE" binary array=${SIZE} format='%double' with lines
+END_OF_SCRIPT

exercise1/theory.typ

@@ -0,0 +1,70 @@
+#import "@preview/cetz:0.3.2";
+#import "@preview/cetz-plot:0.1.1": plot
+#import "@preview/physica:0.9.4": *
+#import "@preview/plotsy-3d:0.1.0": plot-3d-parametric-surface
+#import "@preview/fletcher:0.5.4" as fletcher: diagram, edge, node
+
+#set page(paper: "a4", margin: (x: 2.6cm, y: 2.8cm), numbering: "1 : 1")
+#set par(justify: true, leading: 0.52em)
+
+#let FONT_SIZE = 18pt;
+#set text(font: "FreeSerif", size: FONT_SIZE, lang: "us")
+#show math.equation: set text(font: "Euler Math", size: (FONT_SIZE * 1.0), lang: "en")
+
+#set heading(numbering: none)
+#show heading.where(level: 1): it => {
+  rect(inset: FONT_SIZE / 2)[#it]
+}
+
+// Dracula palette
+#let bg = rgb("#282a36")
+#let fg = rgb("#f8f8f2")
+#let sec = rgb("#44475a")
+
+
+#align(center)[
+  #text(size: FONT_SIZE * 2, weight: "bold")[#underline[exercise 0]]
+]
+
+these are my solutions to exercise set 1 of TDT4200.
+
+this document was created using
+#link("https://typst.app/")[#text(blue.darken(5%))[typst]].
+
+#v(42pt)
+
+#outline(title: none)
+
+#v(42pt)
+
+= macros
+
+macros define programmer-friendly syntax that is preprocessed at compile time,
+thus incurring no performance penality. they constitute a fundamental part in
+meta programming (see lisp). excessive use of macros may obfuscate semantics,
+but can often be used to make the code clearer to read. in this case, it helps
+simplify the indexing of the buffers for easier to read calculations.
+
+= boundary conditions
+
+the boundary condition could be changed from reflection to wrapping around. this
+can be achieved by setting the left boundary to be the right-most point, and the
+right boundary to be the left-most point.
+
+= no memory
+
+if you don't allocate memory in T1, the buffers will be unallocated and you will
+get indexing errors as you try to access the buffers using the predefined
+macros.
+
+basically, you segfault.
+
+= `float const *a` vs `float *const b`
+
+- `float const *a` is a constant pointer to some memory address, assigned to the
+  constant name `a`. since `a` is constant, it cannot be reassigned.
+- `float *const b` is a pointer to some constant value, assigned to the
+  variable `b`. since `b` points to a constant, the dereferenced value cannot be
+  changed.
+
+these are not the same.

exercise1/wave_1d.c

@@ -0,0 +1,120 @@
+#define _XOPEN_SOURCE 600
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+// Option to change numerical precision.
+typedef int64_t int_t;
+typedef double real_t;
+
+// Simulation parameters: size, step count, and how often to save the state.
+const int_t
+    N = 1024,
+    max_iteration = 4000,
+    snapshot_freq = 10;
+
+// Wave equation parameters, time step is derived from the space step.
+const real_t
+    c = 1.0,
+    dx = 1.0;
+real_t
+    dt;
+
+// Buffers for three time steps, indexed with 2 ghost points for the boundary.
+real_t
+    *buffers[3] = { NULL, NULL, NULL };
+
+#define U_prv(i) buffers[0][(i) + 1]
+#define U(i) buffers[1][(i) + 1]
+#define U_nxt(i) buffers[2][(i) + 1]
+
+// Save the present time step in a numbered file under 'data/'.
+void domain_save(int_t step) {
+  char filename[256];
+  sprintf(filename, "data/%.5ld.dat", step);
+  FILE *out = fopen(filename, "wb");
+  fwrite(&U(0), sizeof(real_t), N, out);
+  fclose(out);
+}
+
+// TASK: T1
+// Set up our three buffers, fill two with an initial cosine wave,
+// and set the time step.
+void domain_initialize(void) {
+  // allocate memory
+  for (int i = 0; i < 3; i++)
+    buffers[i] = malloc((N + 2) * sizeof(real_t));
+
+  // apply initial condition
+  real_t x = 0.0f;
+  for (int i = 0; i < N; i++) {
+    U_prv(i) = cos(2 * M_PI * x);
+    U(i) = U_prv(i);
+    U_nxt(i) = 0.0f;
+    x += (real_t)dx / (real_t)N;
+  }
+
+  // set the time-step dt according to CFL
+  dt = dx / c;
+}
+
+// TASK T2:
+// Return the memory to the OS.
+void domain_finalize(void) {
+  for (int i = 0; i < 3; i++) {
+    free(buffers[i]);
+    buffers[i] = NULL;
+  }
+}
+
+// TASK: T3
+// Rotate the time step buffers.
+void rotate_buffers(void) {
+  real_t *temp = buffers[0];
+  buffers[0] = buffers[1];
+  buffers[1] = buffers[2];
+  buffers[2] = temp;
+}
+
+// TASK: T4
+// Derive step t+1 from steps t and t-1.
+void time_step(void) {
+  for (int i = 0; i < N; i++) {
+    U_nxt(i) = -U_prv(i) + 2 * U(i) + dt * dt * c * c / (dx * dx) * (U(i - 1) + U(i + 1) - 2 * U(i));
+  }
+}
+
+// TASK: T5
+// Neumann (reflective) boundary condition.
+void boundary(void) {
+  U_prv(-1) = U_prv(1);
+  U(-1) = U(1);
+  U_nxt(-1) = U_nxt(1);
+  U_prv(N) = U_prv(N - 2);
+  U(N) = U(N - 2);
+  U_nxt(N) = U_nxt(N - 2);
+}
+
+// TASK: T6
+// Main time integration.
+void simulate(void) {
+  int_t iteration = 0;
+  while (iteration < max_iteration) {
+    boundary();
+    time_step();
+    rotate_buffers();
+    domain_save(iteration / snapshot_freq);
+    iteration++;
+  }
+}
+
+int main(void) {
+  domain_initialize();
+
+  simulate();
+
+  domain_finalize();
+  exit(EXIT_SUCCESS);
+}

exercise2/.gitignore

@@ -0,0 +1 @@
+out/

exercise2/Makefile

@@ -0,0 +1,33 @@
+CC      := gcc
+CFLAGS  := -Wall -Wextra -std=c17
+SRC     := mandel_mpi.c
+TARGET  := $(SRC:.c=)
+OUTDIR  := out
+OUT     := $(OUTDIR)/$(TARGET)
+ARGS    := 1
+NBENCH  := 1
+NPROC   := 3
+
+.PHONY: all clean run parallel
+
+all: clean parallel show
+
+$(OUT): $(SRC)
+	$(CC) $(CFLAGS) -o $(OUT) $<
+
+run: $(OUT)
+	cd $(OUTDIR) && ./$(TARGET) $(ARGS)
+
+clean:
+	rm -rf $(OUTDIR)/*
+
+time: $(OUT)
+	python3 bench.py './$(OUT) 0' $(NBENCH)
+
+parallel: $(SRC)
+	mpicc -o $(OUT) $(SRC)
+	cd $(OUTDIR) && mpirun -np $(NPROC) $(TARGET) $(ARGS)
+
+show: $(OUTDIR)/mandel.bmp
+	feh $<
+

exercise2/bench.py

@@ -0,0 +1,24 @@
+#!/usr/bin/env python3
+
+from statistics import mean
+from subprocess import DEVNULL, run
+from sys import argv, exit
+from time import perf_counter
+
+
+def benchmark(cmd, N=10):
+    times = []
+    for _ in range(N):
+        start = perf_counter()
+        run(cmd, shell=True, stdout=DEVNULL, stderr=DEVNULL)
+        times.append(perf_counter() - start)
+    return mean(times)
+
+
+if __name__ == "__main__":
+    if len(argv) < 2:
+        print("usage: python3 bench.py 'cmd' [N]")
+        exit(1)
+    cmd = argv[1]
+    N = int(argv[2]) if len(argv) > 2 else 10
+    print(f"{benchmark(cmd, N):.4f}")

exercise2/mandel_cpu.c

@@ -0,0 +1,100 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define XSIZE 2560
+#define YSIZE 2048
+
+#define MAXITER 255
+
+double xleft=-2.01;
+double xright=1;
+double yupper,ylower;
+double ycenter=1e-6;
+double step;
+
+int pixel[XSIZE*YSIZE];
+
+#define PIXEL(i,j) ((i)+(j)*XSIZE)
+
+typedef struct {
+	double real,imag;
+} complex_t;
+
+void calculate() {
+	for(int i=0;i<XSIZE;i++) {
+		for(int j=0;j<YSIZE;j++) {
+			/* Calculate the number of iterations until divergence for each pixel.
+			   If divergence never happens, return MAXITER */
+			complex_t c,z,temp;
+			int iter=0;
+			c.real = (xleft + step*i);
+			c.imag = (ylower + step*j);
+			z = c;
+			while(z.real*z.real + z.imag*z.imag < 4) {
+				temp.real = z.real*z.real - z.imag*z.imag + c.real;
+				temp.imag = 2*z.real*z.imag + c.imag;
+				z = temp;
+				if(++iter==MAXITER) break;
+			}
+			pixel[PIXEL(i,j)]=iter;
+		}
+	}
+}
+
+typedef unsigned char uchar;
+
+/* save 24-bits bmp file, buffer must be in bmp format: upside-down */
+void savebmp(char *name,uchar *buffer,int x,int y) {
+	FILE *f=fopen(name,"wb");
+	if(!f) {
+		printf("Error writing image to disk.\n");
+		return;
+	}
+	unsigned int size=x*y*3+54;
+	uchar header[54]={'B','M',size&255,(size>>8)&255,(size>>16)&255,size>>24,0,
+		0,0,0,54,0,0,0,40,0,0,0,x&255,x>>8,0,0,y&255,y>>8,0,0,1,0,24,0,0,0,0,0,0,
+		0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+	fwrite(header,1,54,f);
+	fwrite(buffer,1,XSIZE*YSIZE*3,f);
+	fclose(f);
+}
+
+/* given iteration number, set a colour */
+void fancycolour(uchar *p,int iter) {
+	if(iter==MAXITER);
+	else if(iter<8) { p[0]=128+iter*16; p[1]=p[2]=0; }
+	else if(iter<24) { p[0]=255; p[1]=p[2]=(iter-8)*16; }
+	else if(iter<160) { p[0]=p[1]=255-(iter-24)*2; p[2]=255; }
+	else { p[0]=p[1]=(iter-160)*2; p[2]=255-(iter-160)*2; }
+}
+
+int main(int argc,char **argv) {
+	if(argc==1) {
+		puts("Usage: MANDEL n");
+		puts("n decides whether image should be written to disk (1=yes, 0=no)");
+		return 0;
+	}
+	/* Calculate the range in the y-axis such that we preserve the
+	   aspect ratio */
+	step=(xright-xleft)/XSIZE;
+	yupper=ycenter+(step*YSIZE)/2;
+	ylower=ycenter-(step*YSIZE)/2;
+
+	calculate();
+
+	if(strtol(argv[1],NULL,10)!=0) {
+		/* create nice image from iteration counts. take care to create it upside
+		   down (bmp format) */
+		unsigned char *buffer=calloc(XSIZE*YSIZE*3,1);
+		for(int i=0;i<XSIZE;i++) {
+			for(int j=0;j<YSIZE;j++) {
+				int p=((YSIZE-j-1)*XSIZE+i)*3;
+				fancycolour(buffer+p,pixel[PIXEL(i,j)]);
+			}
+		}
+		/* write image to disk */
+		savebmp("mandel2.bmp",buffer,XSIZE,YSIZE);
+	}
+	return 0;
+}

exercise2/mandel_mpi.c

@@ -0,0 +1,158 @@
+#include <math.h>
+#include <mpi.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#define XSIZE 2560
+#define YSIZE 2048
+
+#define MAXITER 255
+
+double xleft = -2.01;
+double xright = 1;
+double yupper, ylower;
+double ycenter = 1e-6;
+double step;
+
+int pixel[XSIZE * YSIZE];
+
+#define PIXEL(i, j) ((i) + (j) * XSIZE)
+
+typedef struct {
+  double real, imag;
+} complex_t;
+
+int rank, size;
+bool save = false;
+
+void calculate() {
+  for (int i = 0; i < XSIZE; i++) {
+    for (int j = 0; j < YSIZE; j++) {
+      /* Calculate the number of iterations until divergence for each pixel.
+         If divergence never happens, return MAXITER */
+      complex_t c, z, temp;
+      int iter = 0;
+      c.real = (xleft + step * i);
+      c.imag = (ylower + step * j);
+      z = c;
+      while (z.real * z.real + z.imag * z.imag < 4) {
+        temp.real = z.real * z.real - z.imag * z.imag + c.real;
+        temp.imag = 2 * z.real * z.imag + c.imag;
+        z = temp;
+        if (++iter == MAXITER)
+          break;
+      }
+      pixel[PIXEL(i, j)] = iter;
+    }
+  }
+}
+
+typedef unsigned char uchar;
+
+/* save 24-bits bmp file, buffer must be in bmp format: upside-down */
+void savebmp(char *name, uchar *buffer, int x, int y) {
+  FILE *f = fopen(name, "wb");
+  if (!f) {
+    printf("Error writing image to disk.\n");
+    return;
+  }
+  unsigned int size = x * y * 3 + 54;
+  uchar header[54] = { 'B', 'M', size & 255, (size >> 8) & 255, (size >> 16) & 255, size >> 24, 0,
+                       0, 0, 0, 54, 0, 0, 0, 40, 0, 0, 0, x & 255, x >> 8, 0, 0, y & 255, y >> 8, 0, 0, 1, 0, 24, 0, 0, 0, 0, 0, 0,
+                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
+  fwrite(header, 1, 54, f);
+  fwrite(buffer, 1, XSIZE * YSIZE * 3, f);
+  fclose(f);
+}
+
+/* given iteration number, set a colour */
+void fancycolour(uchar *p, int iter) {
+  if (iter == MAXITER)
+    ;
+  else if (iter < 8) {
+    p[0] = 128 + iter * 16;
+    p[1] = p[2] = 0;
+  } else if (iter < 24) {
+    p[0] = 255;
+    p[1] = p[2] = (iter - 8) * 16;
+  } else if (iter < 160) {
+    p[0] = p[1] = 255 - (iter - 24) * 2;
+    p[2] = 255;
+  } else {
+    p[0] = p[1] = (iter - 160) * 2;
+    p[2] = 255 - (iter - 160) * 2;
+  }
+}
+
+int main(int argc, char **argv) {
+  if (argc <= 1) {
+    printf("usage: mandel_mpi [y/n]\n");
+    printf("y/n denotes if you want to save the resulting\n");
+    printf("calculations to mandel.bmp or not.\n");
+    return 1;
+  } else {
+    char arg = argv[1][0];
+    save = arg == 'y' || arg == '1';
+  }
+  MPI_Init(&argc, &argv);
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &size);
+
+#undef YSIZE
+#define YSIZE (2048 / size)
+
+  step = (xright - xleft) / XSIZE;
+  double old_yupper = ycenter + (step * YSIZE * size) / 2;
+  double old_ylower = ycenter - (step * YSIZE * size) / 2;
+  double delta = (old_yupper - old_ylower) / size;
+
+  ylower = old_ylower + rank * delta;
+  yupper = old_ylower + (rank + 1) * delta;
+  ycenter = (ylower + yupper) / 2;
+
+  calculate();
+
+  printf("after calculate %d\n", rank);
+
+  const int bufsize = XSIZE * YSIZE * 3;
+  unsigned char *buffer = calloc(bufsize, 1);
+  for (int i = 0; i < XSIZE; i++) {
+    for (int j = 0; j < YSIZE; j++) {
+      int p = ((YSIZE - j - 1) * XSIZE + i) * 3;
+      fancycolour(buffer + p, pixel[PIXEL(i, j)]);
+    }
+  }
+
+  printf("after buffer %d\n", rank);
+
+  if (rank != 0)
+    MPI_Send(buffer, bufsize, MPI_UNSIGNED_CHAR, 0, 0, MPI_COMM_WORLD);
+
+  printf("hello again from %d\n", rank);
+
+  // use process 0 as central process: gather calulations and output image
+  if (rank == 0) {
+    unsigned char *total = calloc(bufsize, size);
+    memcpy(total + (size - 1) * bufsize, buffer, bufsize);
+    printf("after memcpy 0\n");
+
+    for (int i = 1; i < size; i++) {
+      MPI_Recv(total + (size - 1 - i) * bufsize, bufsize, // write reverse order
+               MPI_UNSIGNED_CHAR, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      printf("after direct receive %d\n", i);
+    }
+    printf("after loop\n");
+
+#undef YSIZE
+#define YSIZE 2048
+
+    if (save)
+      savebmp("mandel.bmp", total, XSIZE, YSIZE);
+    printf("after save %d\n", save);
+  }
+
+  MPI_Finalize();
+  return 0;
+}

flake.nix

@@ -16,8 +16,17 @@
         buildInputs = with pkgs; [
           gcc
           gnumake
-          mpi
+          imagemagick_light
+          typst
+          typstyle
+          zip
+          unzip
           feh
+          # exercise 1: 1D wave equation plotting utils
+          gnuplot
+          ffmpeg
+          mpi
+          python3
         ];
         shellHook = ''
           echo welcome!