tasks/c/pngwriter.c

+#include <png.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <limits.h>
+#include "pngwriter.h"
+
+/* Datatype for RGB pixel */
+typedef struct {
+    uint8_t red;
+    uint8_t green;
+    uint8_t blue;
+} pixel_t;
+
+void cmap(const int value, const int max_val, pixel_t *pix);
+
+/* Heat colormap from black to white */
+// *INDENT-OFF*
+static int heat_colormap[256][3] = {
+    {  0,   0,   0}, { 35,   0,   0}, { 52,   0,   0}, { 60,   0,   0},
+    { 63,   1,   0}, { 64,   2,   0}, { 68,   5,   0}, { 69,   6,   0},
+    { 72,   8,   0}, { 74,  10,   0}, { 77,  12,   0}, { 78,  14,   0},
+    { 81,  16,   0}, { 83,  17,   0}, { 85,  19,   0}, { 86,  20,   0},
+    { 89,  22,   0}, { 91,  24,   0}, { 92,  25,   0}, { 94,  26,   0},
+    { 95,  28,   0}, { 98,  30,   0}, {100,  31,   0}, {102,  33,   0},
+    {103,  34,   0}, {105,  35,   0}, {106,  36,   0}, {108,  38,   0},
+    {109,  39,   0}, {111,  40,   0}, {112,  42,   0}, {114,  43,   0},
+    {115,  44,   0}, {117,  45,   0}, {119,  47,   0}, {119,  47,   0},
+    {120,  48,   0}, {122,  49,   0}, {123,  51,   0}, {125,  52,   0},
+    {125,  52,   0}, {126,  53,   0}, {128,  54,   0}, {129,  56,   0},
+    {129,  56,   0}, {131,  57,   0}, {132,  58,   0}, {134,  59,   0},
+    {134,  59,   0}, {136,  61,   0}, {137,  62,   0}, {137,  62,   0},
+    {139,  63,   0}, {139,  63,   0}, {140,  65,   0}, {142,  66,   0},  
+    {148,  71,   0}, {149,  72,   0}, {149,  72,   0}, {151,  73,   0},
+    {151,  73,   0}, {153,  75,   0}, {153,  75,   0}, {154,  76,   0},
+    {154,  76,   0}, {154,  76,   0}, {156,  77,   0}, {156,  77,   0},
+    {157,  79,   0}, {157,  79,   0}, {159,  80,   0}, {159,  80,   0},
+    {159,  80,   0}, {160,  81,   0}, {160,  81,   0}, {162,  82,   0},
+    {162,  82,   0}, {163,  84,   0}, {163,  84,   0}, {165,  85,   0},
+    {165,  85,   0}, {166,  86,   0}, {166,  86,   0}, {166,  86,   0},
+    {168,  87,   0}, {168,  87,   0}, {170,  89,   0}, {170,  89,   0},
+    {171,  90,   0}, {171,  90,   0}, {173,  91,   0}, {173,  91,   0},
+    {174,  93,   0}, {174,  93,   0}, {176,  94,   0}, {176,  94,   0},
+    {177,  95,   0}, {177,  95,   0}, {179,  96,   0}, {179,  96,   0},
+    {180,  98,   0}, {182,  99,   0}, {182,  99,   0}, {183, 100,   0},
+    {183, 100,   0}, {185, 102,   0}, {185, 102,   0}, {187, 103,   0},
+    {187, 103,   0}, {188, 104,   0}, {188, 104,   0}, {190, 105,   0},
+    {191, 107,   0}, {191, 107,   0}, {193, 108,   0}, {193, 108,   0},
+    {194, 109,   0}, {196, 110,   0}, {196, 110,   0}, {197, 112,   0},
+    {197, 112,   0}, {199, 113,   0}, {200, 114,   0}, {200, 114,   0},
+    {202, 116,   0}, {202, 116,   0}, {204, 117,   0}, {205, 118,   0},
+    {205, 118,   0}, {207, 119,   0}, {208, 121,   0}, {208, 121,   0},
+    {210, 122,   0}, {211, 123,   0}, {211, 123,   0}, {213, 124,   0},
+    {214, 126,   0}, {214, 126,   0}, {216, 127,   0}, {217, 128,   0},
+    {217, 128,   0}, {219, 130,   0}, {221, 131,   0}, {221, 131,   0},
+    {222, 132,   0}, {224, 133,   0}, {224, 133,   0}, {225, 135,   0},
+    {227, 136,   0}, {227, 136,   0}, {228, 137,   0}, {230, 138,   0},
+    {230, 138,   0}, {231, 140,   0}, {233, 141,   0}, {233, 141,   0},
+    {234, 142,   0}, {236, 144,   0}, {236, 144,   0}, {238, 145,   0},
+    {239, 146,   0}, {241, 147,   0}, {241, 147,   0}, {242, 149,   0},
+    {244, 150,   0}, {244, 150,   0}, {245, 151,   0}, {247, 153,   0},
+    {247, 153,   0}, {248, 154,   0}, {250, 155,   0}, {251, 156,   0},
+    {251, 156,   0}, {253, 158,   0}, {255, 159,   0}, {255, 159,   0},
+    {255, 160,   0}, {255, 161,   0}, {255, 163,   0}, {255, 163,   0},
+    {255, 164,   0}, {255, 165,   0}, {255, 167,   0}, {255, 167,   0},
+    {255, 168,   0}, {255, 169,   0}, {255, 169,   0}, {255, 170,   0},
+    {255, 172,   0}, {255, 173,   0}, {255, 173,   0}, {255, 174,   0},
+    {255, 175,   0}, {255, 177,   0}, {255, 178,   0}, {255, 179,   0},
+    {255, 181,   0}, {255, 181,   0}, {255, 182,   0}, {255, 183,   0},
+    {255, 184,   0}, {255, 187,   7}, {255, 188,  10}, {255, 189,  14},
+    {255, 191,  18}, {255, 192,  21}, {255, 193,  25}, {255, 195,  29},
+    {255, 197,  36}, {255, 198,  40}, {255, 200,  43}, {255, 202,  51},
+    {255, 204,  54}, {255, 206,  61}, {255, 207,  65}, {255, 210,  72},
+    {255, 211,  76}, {255, 214,  83}, {255, 216,  91}, {255, 219,  98},
+    {255, 221, 105}, {255, 223, 109}, {255, 225, 116}, {255, 228, 123},
+    {255, 232, 134}, {255, 234, 142}, {255, 237, 149}, {255, 239, 156},
+    {255, 240, 160}, {255, 243, 167}, {255, 246, 174}, {255, 248, 182},
+    {255, 249, 185}, {255, 252, 193}, {255, 253, 196}, {255, 255, 204},
+    {255, 255, 207}, {255, 255, 211}, {255, 255, 218}, {255, 255, 222},
+    {255, 255, 225}, {255, 255, 229}, {255, 255, 233}, {255, 255, 236},
+    {255, 255, 240}, {255, 255, 244}, {255, 255, 247}, {255, 255, 255}
+};
+// *INDENT-ON*
+/* Save the two dimensional array as a png image
+ * Arguments:
+ *   double *data is a pointer to an array of nx*ny values
+ *   int nx is the number of COLUMNS to be written
+ *   int ny is the number of ROWS to be written
+ *   char *fname is the name of the picture
+ *   char lang is either 'c' or 'f' denoting the memory
+ *             layout. That is, if 'f' is given, then rows
+ *             and columns are swapped.
+ */
+int save_png(int *data, const int height, const int width, const char *fname)
+{
+    FILE *fp;
+    png_structp pngstruct_ptr = NULL;
+    png_infop pnginfo_ptr = NULL;
+    png_byte **row_pointers = NULL;
+    int i, j, max_value;
+
+    /* Determine the maximum value */
+    max_value = INT_MIN;
+    for (i = 0; i < height * width; i++) {
+        max_value = data[i] > max_value ? data[i] : max_value;
+    }
+
+    /* Default return status is failure */
+    int status = -1;
+
+    int pixel_size = 3;
+    int depth = 8;
+
+    fp = fopen(fname, "wb");
+    if (fp == NULL) {
+        goto fopen_failed;
+    }
+
+    pngstruct_ptr =
+        png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+
+    if (pngstruct_ptr == NULL) {
+        goto pngstruct_create_failed;
+    }
+
+    pnginfo_ptr = png_create_info_struct(pngstruct_ptr);
+
+    if (pnginfo_ptr == NULL) {
+        goto pnginfo_create_failed;
+    }
+
+    if (setjmp(png_jmpbuf(pngstruct_ptr))) {
+        goto setjmp_failed;
+    }
+
+    png_set_IHDR(pngstruct_ptr, pnginfo_ptr, (size_t) width,
+                 (size_t) height, depth, PNG_COLOR_TYPE_RGB,
+                 PNG_INTERLACE_NONE,
+                 PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+
+    row_pointers = png_malloc(pngstruct_ptr, height * sizeof(png_byte *));
+
+    for (i = 0; i < height; i++) {
+        png_byte *row = png_malloc(pngstruct_ptr,
+                                   sizeof(uint8_t) * width * pixel_size);
+        row_pointers[i] = row;
+
+        // Branch according to the memory layout
+        for (j = 0; j < width; j++) {
+            pixel_t pixel;
+            // Scale the values so that values between
+            // 0 and 100 degrees are mapped to values
+            // between 0 and 255
+            cmap(data[j + i * width], max_value, &pixel);
+            *row++ = pixel.red;
+            *row++ = pixel.green;
+            *row++ = pixel.blue;
+        }
+    }
+
+    png_init_io(pngstruct_ptr, fp);
+    png_set_rows(pngstruct_ptr, pnginfo_ptr, row_pointers);
+    png_write_png(pngstruct_ptr, pnginfo_ptr,
+                  PNG_TRANSFORM_IDENTITY, NULL);
+
+    status = 0;
+
+    for (i = 0; i < height; i++) {
+        png_free(pngstruct_ptr, row_pointers[i]);
+    }
+    png_free(pngstruct_ptr, row_pointers);
+
+setjmp_failed:
+pnginfo_create_failed:
+    png_destroy_write_struct(&pngstruct_ptr, &pnginfo_ptr);
+pngstruct_create_failed:
+    fclose(fp);
+fopen_failed:
+    return status;
+}
+
+/* This routine sets the RGB values for the pixel_t structure using
+ * the colormap data heat_colormap. If the value is outside the
+ * acceptable png values 0,255 blue or red color is used instead. */
+void cmap(const int value, const int max_value, pixel_t *pix)
+{
+    int ival;
+
+    ival = ((double) value / (double) max_value) * 255;
+    /* Check for wrong values */
+    if (ival > 255) {
+        ival = 255;
+    }
+    if (ival < 0) {
+        ival = 0;
+    }
+
+    /* Pick color from heat_colormap */
+    pix->red = heat_colormap[ival][0];
+    pix->green = heat_colormap[ival][1];
+    pix->blue = heat_colormap[ival][2];
+}

tasks/c/pngwriter.h

+#ifndef PNGWRITER_H_
+#define PNGWRITER_H_
+
+int save_png(int *data, const int nx, const int ny, const char *fname);
+
+#endif

tasks/c/solution/Makefile

+ifeq ($(COMP),)
+COMP=cray
+endif
+
+COMMONDIR=../../../heat/common
+
+ifeq ($(COMP),cray)
+FC=ftn
+CC=cc
+FCFLAGS=-O3
+CCFLAGS=-O3 -I/appl/opt/libpng/include -I$(COMMONDIR)
+LDFLAGS=-L/appl/opt/libpng/lib
+LIBS=-lpng -lz
+endif
+
+ifeq ($(COMP),gnu)
+FC=gfortran
+CC=gcc
+FCFLAGS=-fopenmp -O3 -Wall
+CCFLAGS=-fopenmp -O3 -Wall -I$(COMMONDIR)
+LDFLAGS=
+LIBS=-lpng
+endif
+
+EXE=mandelbrot
+OBJS=mandelbrot.o
+OBJS_PNG=pngwriter.o
+CORRECT_OBJS_PNG=$(COMMONDIR)/pngwriter.o
+
+all: $(EXE)
+
+$(COMMONDIR)/pngwriter.o: $(COMMONDIR)/pngwriter.c $(COMMONDIR)/pngwriter.h
+mandelbrot.o: mandelbrot.c
+
+pngwriter.o: pngwriter.c pngwriter.h
+	$(CC) $(CCFLAGS) -c -o $@ $<
+
+$(EXE): $(OBJS) $(OBJS_PNG)
+	$(CC) $(CCFLAGS) $(OBJS) $(OBJS_PNG) -o $@ $(LDFLAGS) $(LIBS)
+
+%.o: %.F90
+	$(FC) $(FCFLAGS) -c $< -o $@
+
+%.o: %.c
+	$(CC) $(CCFLAGS) -c $< -o $@
+
+.PHONY: clean
+clean:
+	-/bin/rm -f $(EXE) a.out *.o *.mod *.png *~ mandelbrot.png

tasks/c/solution/mandelbrot.c

+/*
+ * This example is based on the code of Andrew V. Adinetz
+ * https://github.com/canonizer/mandelbrot-dyn
+ * Licensed under The MIT License
+ */
+
+#include <omp.h>
+#include <complex.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "pngwriter.h"
+
+// Maximum number of iterations
+const int MAX_ITER_COUNT = 512;
+// Marker for different iteration counts
+const int DIFF_ITER_COUNT = -1;
+// Maximum recursion depth
+const int MAX_DEPTH = 6;
+// Region size below which do per-pixel
+const int MIN_SIZE = 32;
+// Subdivision factor along each axis
+const int SUBDIV = 4;
+
+
+// |z|^2 of a complex number z
+float abs2(complex v)
+{
+    return creal(v) * creal(v) + cimag(v) * cimag(v);
+}
+
+// The kernel to count per-pixel values of the portion of the Mandelbrot set
+// Does not need to be edited
+int kernel(int w, int h, complex cmin, complex cmax,
+           int x, int y)
+{
+    complex dc = cmax - cmin;
+    float fx = (float)x / w;
+    float fy = (float)y / h;
+    complex c = cmin + fx * creal(dc) + fy * cimag(dc) * I;
+    int iteration = 0;
+    complex z = c;
+    while (iteration < MAX_ITER_COUNT && abs2(z) < 2 * 2) {
+        z = z * z + c;
+        iteration++;
+    }
+    return iteration;
+}
+
+
+/* Computes the Mandelbrot image recursively
+ * At each call, the image is divided into smaller blocks (by a factor of
+ * subdiv), and the function is called recursively with arguments corresponding
+ * to subblock. When maximum recursion depth is reached or size of block
+ * is smaller than predefined minimum, one starts to calculate actual pixel
+ * values
+ *
+ * - - - - - - - -           -----       -----
+ * |             |           |   |       |   |
+ * |             |           -----       -----
+ * |             |  -->                         -->   ...
+ * |             |           -----       -----
+ * |             |           |   |       |   |
+ * |             |           -----       -----
+ * ---------------
+ */
+void mandelbrot_block(int *iter_counts, int w, int h, complex cmin,
+                      complex cmax, int x0, int y0, int d, int depth)
+{
+
+// TODO Parallelize the recursive function call
+// with OpenMP tasks
+
+    int block_size = d / SUBDIV;
+    if (depth + 1 < MAX_DEPTH && block_size > MIN_SIZE) {
+        // Subdivide recursively
+        for (int i = 0; i < SUBDIV; i++) {
+            for (int j = 0; j < SUBDIV; j++) {
+                #pragma omp task
+                mandelbrot_block(iter_counts, w, h, cmin, cmax,
+                                 x0 + i * block_size, y0 + j * block_size,
+                                 d / SUBDIV, depth + 1);
+            }
+        }
+    } else {
+        // Last recursion level reached, calculate the values
+        for (int i = x0; i < x0 + d; i++) {
+            for (int j = y0; j < y0 + d; j++) {
+                iter_counts[j * w + i] = kernel(w, h, cmin, cmax, i, j);
+            }
+        }
+    }
+}
+
+
+int main(int argc, char **argv)
+{
+    // Picture size, should be power of two
+    const int w = 512;
+    const int h = w;
+    int *iter_counts;
+
+    complex cmin, cmax;
+
+    int pic_bytes = w * h * sizeof(int);
+    iter_counts = (int *)malloc(pic_bytes);
+
+    cmin = -1.5 + -1.0 * I;
+    cmax = 0.5 + 1.0 * I;
+
+    double t1 = omp_get_wtime();
+
+// TODO create parallel region. How many threads should be calling
+// mandelbrot_block in this uppermost level?
+    #pragma omp parallel
+    #pragma omp single
+    {
+        mandelbrot_block(iter_counts, w, h, cmin, cmax,
+                         0, 0, w, 1);
+    }
+    double t2 = omp_get_wtime();
+
+    // Save the image to a PNG file
+    save_png(iter_counts, w, h, "mandelbrot.png");
+
+    double walltime = t2 - t1;
+    // Print the timings
+    printf("Mandelbrot set computed in %.3lf s, at %.3lf Mpix/s\n",
+           walltime, h * w * 1e-6 / walltime);
+
+    free(iter_counts);
+    return 0;
+}
+

tasks/c/solution/pngwriter.c

+#include <png.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <limits.h>
+#include "pngwriter.h"
+
+/* Datatype for RGB pixel */
+typedef struct {
+    uint8_t red;
+    uint8_t green;
+    uint8_t blue;
+} pixel_t;
+
+void cmap(const int value, const int max_val, pixel_t *pix);
+
+/* Heat colormap from black to white */
+// *INDENT-OFF*
+static int heat_colormap[256][3] = {
+    {  0,   0,   0}, { 35,   0,   0}, { 52,   0,   0}, { 60,   0,   0},
+    { 63,   1,   0}, { 64,   2,   0}, { 68,   5,   0}, { 69,   6,   0},
+    { 72,   8,   0}, { 74,  10,   0}, { 77,  12,   0}, { 78,  14,   0},
+    { 81,  16,   0}, { 83,  17,   0}, { 85,  19,   0}, { 86,  20,   0},
+    { 89,  22,   0}, { 91,  24,   0}, { 92,  25,   0}, { 94,  26,   0},
+    { 95,  28,   0}, { 98,  30,   0}, {100,  31,   0}, {102,  33,   0},
+    {103,  34,   0}, {105,  35,   0}, {106,  36,   0}, {108,  38,   0},
+    {109,  39,   0}, {111,  40,   0}, {112,  42,   0}, {114,  43,   0},
+    {115,  44,   0}, {117,  45,   0}, {119,  47,   0}, {119,  47,   0},
+    {120,  48,   0}, {122,  49,   0}, {123,  51,   0}, {125,  52,   0},
+    {125,  52,   0}, {126,  53,   0}, {128,  54,   0}, {129,  56,   0},
+    {129,  56,   0}, {131,  57,   0}, {132,  58,   0}, {134,  59,   0},
+    {134,  59,   0}, {136,  61,   0}, {137,  62,   0}, {137,  62,   0},
+    {139,  63,   0}, {139,  63,   0}, {140,  65,   0}, {142,  66,   0},  
+    {148,  71,   0}, {149,  72,   0}, {149,  72,   0}, {151,  73,   0},
+    {151,  73,   0}, {153,  75,   0}, {153,  75,   0}, {154,  76,   0},
+    {154,  76,   0}, {154,  76,   0}, {156,  77,   0}, {156,  77,   0},
+    {157,  79,   0}, {157,  79,   0}, {159,  80,   0}, {159,  80,   0},
+    {159,  80,   0}, {160,  81,   0}, {160,  81,   0}, {162,  82,   0},
+    {162,  82,   0}, {163,  84,   0}, {163,  84,   0}, {165,  85,   0},
+    {165,  85,   0}, {166,  86,   0}, {166,  86,   0}, {166,  86,   0},
+    {168,  87,   0}, {168,  87,   0}, {170,  89,   0}, {170,  89,   0},
+    {171,  90,   0}, {171,  90,   0}, {173,  91,   0}, {173,  91,   0},
+    {174,  93,   0}, {174,  93,   0}, {176,  94,   0}, {176,  94,   0},
+    {177,  95,   0}, {177,  95,   0}, {179,  96,   0}, {179,  96,   0},
+    {180,  98,   0}, {182,  99,   0}, {182,  99,   0}, {183, 100,   0},
+    {183, 100,   0}, {185, 102,   0}, {185, 102,   0}, {187, 103,   0},
+    {187, 103,   0}, {188, 104,   0}, {188, 104,   0}, {190, 105,   0},
+    {191, 107,   0}, {191, 107,   0}, {193, 108,   0}, {193, 108,   0},
+    {194, 109,   0}, {196, 110,   0}, {196, 110,   0}, {197, 112,   0},
+    {197, 112,   0}, {199, 113,   0}, {200, 114,   0}, {200, 114,   0},
+    {202, 116,   0}, {202, 116,   0}, {204, 117,   0}, {205, 118,   0},
+    {205, 118,   0}, {207, 119,   0}, {208, 121,   0}, {208, 121,   0},
+    {210, 122,   0}, {211, 123,   0}, {211, 123,   0}, {213, 124,   0},
+    {214, 126,   0}, {214, 126,   0}, {216, 127,   0}, {217, 128,   0},
+    {217, 128,   0}, {219, 130,   0}, {221, 131,   0}, {221, 131,   0},
+    {222, 132,   0}, {224, 133,   0}, {224, 133,   0}, {225, 135,   0},
+    {227, 136,   0}, {227, 136,   0}, {228, 137,   0}, {230, 138,   0},
+    {230, 138,   0}, {231, 140,   0}, {233, 141,   0}, {233, 141,   0},
+    {234, 142,   0}, {236, 144,   0}, {236, 144,   0}, {238, 145,   0},
+    {239, 146,   0}, {241, 147,   0}, {241, 147,   0}, {242, 149,   0},
+    {244, 150,   0}, {244, 150,   0}, {245, 151,   0}, {247, 153,   0},
+    {247, 153,   0}, {248, 154,   0}, {250, 155,   0}, {251, 156,   0},
+    {251, 156,   0}, {253, 158,   0}, {255, 159,   0}, {255, 159,   0},
+    {255, 160,   0}, {255, 161,   0}, {255, 163,   0}, {255, 163,   0},
+    {255, 164,   0}, {255, 165,   0}, {255, 167,   0}, {255, 167,   0},
+    {255, 168,   0}, {255, 169,   0}, {255, 169,   0}, {255, 170,   0},
+    {255, 172,   0}, {255, 173,   0}, {255, 173,   0}, {255, 174,   0},
+    {255, 175,   0}, {255, 177,   0}, {255, 178,   0}, {255, 179,   0},
+    {255, 181,   0}, {255, 181,   0}, {255, 182,   0}, {255, 183,   0},
+    {255, 184,   0}, {255, 187,   7}, {255, 188,  10}, {255, 189,  14},
+    {255, 191,  18}, {255, 192,  21}, {255, 193,  25}, {255, 195,  29},
+    {255, 197,  36}, {255, 198,  40}, {255, 200,  43}, {255, 202,  51},
+    {255, 204,  54}, {255, 206,  61}, {255, 207,  65}, {255, 210,  72},
+    {255, 211,  76}, {255, 214,  83}, {255, 216,  91}, {255, 219,  98},
+    {255, 221, 105}, {255, 223, 109}, {255, 225, 116}, {255, 228, 123},
+    {255, 232, 134}, {255, 234, 142}, {255, 237, 149}, {255, 239, 156},
+    {255, 240, 160}, {255, 243, 167}, {255, 246, 174}, {255, 248, 182},
+    {255, 249, 185}, {255, 252, 193}, {255, 253, 196}, {255, 255, 204},
+    {255, 255, 207}, {255, 255, 211}, {255, 255, 218}, {255, 255, 222},
+    {255, 255, 225}, {255, 255, 229}, {255, 255, 233}, {255, 255, 236},
+    {255, 255, 240}, {255, 255, 244}, {255, 255, 247}, {255, 255, 255}
+};
+// *INDENT-ON*
+/* Save the two dimensional array as a png image
+ * Arguments:
+ *   double *data is a pointer to an array of nx*ny values
+ *   int nx is the number of COLUMNS to be written
+ *   int ny is the number of ROWS to be written
+ *   char *fname is the name of the picture
+ *   char lang is either 'c' or 'f' denoting the memory
+ *             layout. That is, if 'f' is given, then rows
+ *             and columns are swapped.
+ */
+int save_png(int *data, const int height, const int width, const char *fname)
+{
+    FILE *fp;
+    png_structp pngstruct_ptr = NULL;
+    png_infop pnginfo_ptr = NULL;
+    png_byte **row_pointers = NULL;
+    int i, j, max_value;
+
+    /* Determine the maximum value */
+    max_value = INT_MIN;
+    for (i = 0; i < height * width; i++) {
+        max_value = data[i] > max_value ? data[i] : max_value;
+    }
+
+    /* Default return status is failure */
+    int status = -1;
+
+    int pixel_size = 3;
+    int depth = 8;
+
+    fp = fopen(fname, "wb");
+    if (fp == NULL) {
+        goto fopen_failed;
+    }
+
+    pngstruct_ptr =
+        png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+
+    if (pngstruct_ptr == NULL) {
+        goto pngstruct_create_failed;
+    }
+
+    pnginfo_ptr = png_create_info_struct(pngstruct_ptr);
+
+    if (pnginfo_ptr == NULL) {
+        goto pnginfo_create_failed;
+    }
+
+    if (setjmp(png_jmpbuf(pngstruct_ptr))) {
+        goto setjmp_failed;
+    }
+
+    png_set_IHDR(pngstruct_ptr, pnginfo_ptr, (size_t) width,
+                 (size_t) height, depth, PNG_COLOR_TYPE_RGB,
+                 PNG_INTERLACE_NONE,
+                 PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+
+    row_pointers = png_malloc(pngstruct_ptr, height * sizeof(png_byte *));
+
+    for (i = 0; i < height; i++) {
+        png_byte *row = png_malloc(pngstruct_ptr,
+                                   sizeof(uint8_t) * width * pixel_size);
+        row_pointers[i] = row;
+
+        // Branch according to the memory layout
+        for (j = 0; j < width; j++) {
+            pixel_t pixel;
+            // Scale the values so that values between
+            // 0 and 100 degrees are mapped to values
+            // between 0 and 255
+            cmap(data[j + i * width], max_value, &pixel);
+            *row++ = pixel.red;
+            *row++ = pixel.green;
+            *row++ = pixel.blue;
+        }
+    }
+
+    png_init_io(pngstruct_ptr, fp);
+    png_set_rows(pngstruct_ptr, pnginfo_ptr, row_pointers);
+    png_write_png(pngstruct_ptr, pnginfo_ptr,
+                  PNG_TRANSFORM_IDENTITY, NULL);
+
+    status = 0;
+
+    for (i = 0; i < height; i++) {
+        png_free(pngstruct_ptr, row_pointers[i]);
+    }
+    png_free(pngstruct_ptr, row_pointers);
+
+setjmp_failed:
+pnginfo_create_failed:
+    png_destroy_write_struct(&pngstruct_ptr, &pnginfo_ptr);
+pngstruct_create_failed:
+    fclose(fp);
+fopen_failed:
+    return status;
+}
+
+/* This routine sets the RGB values for the pixel_t structure using
+ * the colormap data heat_colormap. If the value is outside the
+ * acceptable png values 0,255 blue or red color is used instead. */
+void cmap(const int value, const int max_value, pixel_t *pix)
+{
+    int ival;
+
+    ival = ((double) value / (double) max_value) * 255;
+    /* Check for wrong values */
+    if (ival > 255) {
+        ival = 255;
+    }
+    if (ival < 0) {
+        ival = 0;
+    }
+
+    /* Pick color from heat_colormap */
+    pix->red = heat_colormap[ival][0];
+    pix->green = heat_colormap[ival][1];
+    pix->blue = heat_colormap[ival][2];
+}

tasks/c/solution/pngwriter.h

+#ifndef PNGWRITER_H_
+#define PNGWRITER_H_
+
+int save_png(int *data, const int nx, const int ny, const char *fname);
+
+#endif

tasks/fortran/Makefile

+ifeq ($(COMP),)
+COMP=cray
+endif
+
+COMMONDIR=../../../heat/common
+
+ifeq ($(COMP),cray)
+FC=ftn
+CC=cc
+FCFLAGS=-O3
+CCFLAGS=-O3 -I/appl/opt/libpng/include -I$(COMMONDIR)
+LDFLAGS=-L/appl/opt/libpng/lib
+LIBS=-lpng -lz
+endif
+
+ifeq ($(COMP),gnu)
+FC=gfortran
+CC=gcc
+FCFLAGS=-fopenmp -O3 -Wall
+CCFLAGS=-fopenmp -O3 -Wall -I$(COMMONDIR)
+LDFLAGS=
+LIBS=-lpng
+endif
+
+EXE=mandelbrot
+OBJS=mandelbrot.o pngwriter_mod.o
+OBJS_PNG=pngwriter.o
+CORRECT_OBJS_PNG=$(COMMONDIR)/pngwriter.o
+
+all: $(EXE)
+
+$(COMMONDIR)/pngwriter.o: $(COMMONDIR)/pngwriter.c $(COMMONDIR)/pngwriter.h
+pngwriter_mod.o: pngwriter_mod.F90
+mandelbrot.o: mandelbrot.F90 pngwriter_mod.o
+
+pngwriter.o: pngwriter.c pngwriter.h
+	$(CC) $(CCFLAGS) -c -o $@ $<
+
+$(EXE): $(OBJS) $(OBJS_PNG)
+	$(FC) $(FCFLAGS) $(OBJS) $(OBJS_PNG) -o $@ $(LDFLAGS) $(LIBS)
+
+%.o: %.F90
+	$(FC) $(FCFLAGS) -c $< -o $@
+
+%.o: %.c
+	$(CC) $(CCFLAGS) -c $< -o $@
+
+.PHONY: clean
+clean:
+	-/bin/rm -f $(EXE) a.out *.o *.mod *.png *~ mandelbrot.png

tasks/fortran/mandelbrot.F90

+program mandelbrot
+
+  use iso_fortran_env, only : REAL64
+  use pngwriter
+  use omp_lib
+
+  implicit none
+
+  integer, parameter :: dp = REAL64
+
+  integer, parameter :: max_iter_count = 512
+  integer, parameter :: max_depth = 6
+  integer, parameter :: min_size = 32
+  integer, parameter :: subdiv = 4
+
+  integer, parameter :: w = 2048
+  integer, parameter :: h = w
+
+  integer, pointer, dimension(:,:) :: iter_counts
+
+  integer :: stat
+
+  complex(dp) :: cmin, cmax
+  real(dp) :: t0, t1
+
+  allocate(iter_counts(w, h))
+
+
+  t0 = omp_get_wtime()
+
+  cmin = (-1.5, -1.0)
+  cmax = (0.5, 1.0)
+
+  ! TODO create parallel region. How many threads should be calling
+  ! mandelbrot_block in this uppermost level?
+
+  call mandelbrot_block(iter_counts, w, h, cmin, cmax, 0, 0, w, 1)
+
+  t1 = omp_get_wtime()
+
+  stat = save_png(iter_counts, h, w, 'mandelbrot.png')
+
+  deallocate(iter_counts)
+
+  write(*,*) 'Mandelbrot set computed in', t1 - t0, 's'
+
+contains
+
+  ! Computes the Mandelbrot image recursively
+  ! At each call, the image is divided into smaller blocks (by a factor of
+  ! subdiv), and the function is called recursively with arguments corresponding
+  ! to subblock. When maximum recursion depth is reached or size of block
+  ! is smaller than predefined minimum, one starts to calculate actual pixel
+  ! values
+  !
+  ! - - - - - - - -           -----       -----
+  ! |             |           |   |       |   |
+  ! |             |           -----       -----
+  ! |             |  -->                         -->   ...
+  ! |             |           -----       -----
+  ! |             |           |   |       |   |
+  ! |             |           -----       -----
+  ! ---------------
+  !
+
+  recursive subroutine mandelbrot_block(iter_counts, w, h, cmin, cmax, &
+       x0, y0, d, depth)
+    implicit none
+
+    integer, pointer, dimension(:,:), intent(inout) :: iter_counts
+    integer, intent(in) :: w, h, x0, y0, d, depth
+    complex(dp), intent(in) :: cmin, cmax
+
+    integer :: block_size, i, j
+
+    ! TODO Parallelize the recursive function call
+    ! with OpenMP tasks
+
+    block_size = d / subdiv
+    if ((depth + 1 < max_depth) .and. (block_size > min_size)) then
+       ! Subdivide recursively
+       do i=0, subdiv - 1
+          do j=0, subdiv - 1
+             call mandelbrot_block(iter_counts, w, h, cmin, cmax, &
+                  x0 + i*block_size, y0 + j*block_size, &
+                  block_size, depth + 1)
+          end do
+       end do
+    else
+       ! Last recursion level reached, calculate the values
+       do j = y0 + 1, y0 + d
+          do i = x0 + 1, x0 + d
+             iter_counts(i, j) = kernel(h, w, cmin, cmax, i-1, j-1)
+          end do
+       end do
+    end if
+
+  end subroutine mandelbrot_block
+
+  ! Calculate iteration count for a pixel.
+  ! This function does not need to be edited
+  integer function kernel(h, w, cmin, cmax, x, y)
+    implicit none
+
+    integer, intent(in) :: h, w, x, y
+    complex(dp), intent(in) :: cmin, cmax
+
+    integer :: iteration
+    complex(dp) :: z, dc, c
+    real(dp) :: fx, fy
+
+    dc = cmax - cmin
+    fx = real(x, dp) / w
+    fy = real(y, dp) / h
+    c = cmin + fx * real(dc) +  fy * imag(dc) * (0.0, 1.0)
+
+    z = c
+    iteration = 0
+
+    do while (iteration < max_iter_count .and. abs(z)**2 < 4)
+       z = z**2 + c
+       iteration = iteration + 1
+    end do
+
+    kernel = iteration
+
+  end function kernel
+
+end program mandelbrot

tasks/fortran/pngwriter.c

+#include <png.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <limits.h>
+#include "pngwriter.h"
+
+/* Datatype for RGB pixel */
+typedef struct {
+    uint8_t red;
+    uint8_t green;
+    uint8_t blue;
+} pixel_t;
+
+void cmap(const int value, const int max_val, pixel_t *pix);
+
+/* Heat colormap from black to white */
+// *INDENT-OFF*
+static int heat_colormap[256][3] = {
+    {  0,   0,   0}, { 35,   0,   0}, { 52,   0,   0}, { 60,   0,   0},
+    { 63,   1,   0}, { 64,   2,   0}, { 68,   5,   0}, { 69,   6,   0},
+    { 72,   8,   0}, { 74,  10,   0}, { 77,  12,   0}, { 78,  14,   0},
+    { 81,  16,   0}, { 83,  17,   0}, { 85,  19,   0}, { 86,  20,   0},
+    { 89,  22,   0}, { 91,  24,   0}, { 92,  25,   0}, { 94,  26,   0},
+    { 95,  28,   0}, { 98,  30,   0}, {100,  31,   0}, {102,  33,   0},
+    {103,  34,   0}, {105,  35,   0}, {106,  36,   0}, {108,  38,   0},
+    {109,  39,   0}, {111,  40,   0}, {112,  42,   0}, {114,  43,   0},
+    {115,  44,   0}, {117,  45,   0}, {119,  47,   0}, {119,  47,   0},
+    {120,  48,   0}, {122,  49,   0}, {123,  51,   0}, {125,  52,   0},
+    {125,  52,   0}, {126,  53,   0}, {128,  54,   0}, {129,  56,   0},
+    {129,  56,   0}, {131,  57,   0}, {132,  58,   0}, {134,  59,   0},
+    {134,  59,   0}, {136,  61,   0}, {137,  62,   0}, {137,  62,   0},
+    {139,  63,   0}, {139,  63,   0}, {140,  65,   0}, {142,  66,   0},  
+    {148,  71,   0}, {149,  72,   0}, {149,  72,   0}, {151,  73,   0},
+    {151,  73,   0}, {153,  75,   0}, {153,  75,   0}, {154,  76,   0},
+    {154,  76,   0}, {154,  76,   0}, {156,  77,   0}, {156,  77,   0},
+    {157,  79,   0}, {157,  79,   0}, {159,  80,   0}, {159,  80,   0},
+    {159,  80,   0}, {160,  81,   0}, {160,  81,   0}, {162,  82,   0},
+    {162,  82,   0}, {163,  84,   0}, {163,  84,   0}, {165,  85,   0},
+    {165,  85,   0}, {166,  86,   0}, {166,  86,   0}, {166,  86,   0},
+    {168,  87,   0}, {168,  87,   0}, {170,  89,   0}, {170,  89,   0},
+    {171,  90,   0}, {171,  90,   0}, {173,  91,   0}, {173,  91,   0},
+    {174,  93,   0}, {174,  93,   0}, {176,  94,   0}, {176,  94,   0},
+    {177,  95,   0}, {177,  95,   0}, {179,  96,   0}, {179,  96,   0},
+    {180,  98,   0}, {182,  99,   0}, {182,  99,   0}, {183, 100,   0},
+    {183, 100,   0}, {185, 102,   0}, {185, 102,   0}, {187, 103,   0},
+    {187, 103,   0}, {188, 104,   0}, {188, 104,   0}, {190, 105,   0},
+    {191, 107,   0}, {191, 107,   0}, {193, 108,   0}, {193, 108,   0},
+    {194, 109,   0}, {196, 110,   0}, {196, 110,   0}, {197, 112,   0},
+    {197, 112,   0}, {199, 113,   0}, {200, 114,   0}, {200, 114,   0},
+    {202, 116,   0}, {202, 116,   0}, {204, 117,   0}, {205, 118,   0},
+    {205, 118,   0}, {207, 119,   0}, {208, 121,   0}, {208, 121,   0},
+    {210, 122,   0}, {211, 123,   0}, {211, 123,   0}, {213, 124,   0},
+    {214, 126,   0}, {214, 126,   0}, {216, 127,   0}, {217, 128,   0},
+    {217, 128,   0}, {219, 130,   0}, {221, 131,   0}, {221, 131,   0},
+    {222, 132,   0}, {224, 133,   0}, {224, 133,   0}, {225, 135,   0},
+    {227, 136,   0}, {227, 136,   0}, {228, 137,   0}, {230, 138,   0},
+    {230, 138,   0}, {231, 140,   0}, {233, 141,   0}, {233, 141,   0},
+    {234, 142,   0}, {236, 144,   0}, {236, 144,   0}, {238, 145,   0},
+    {239, 146,   0}, {241, 147,   0}, {241, 147,   0}, {242, 149,   0},
+    {244, 150,   0}, {244, 150,   0}, {245, 151,   0}, {247, 153,   0},
+    {247, 153,   0}, {248, 154,   0}, {250, 155,   0}, {251, 156,   0},
+    {251, 156,   0}, {253, 158,   0}, {255, 159,   0}, {255, 159,   0},
+    {255, 160,   0}, {255, 161,   0}, {255, 163,   0}, {255, 163,   0},
+    {255, 164,   0}, {255, 165,   0}, {255, 167,   0}, {255, 167,   0},
+    {255, 168,   0}, {255, 169,   0}, {255, 169,   0}, {255, 170,   0},
+    {255, 172,   0}, {255, 173,   0}, {255, 173,   0}, {255, 174,   0},
+    {255, 175,   0}, {255, 177,   0}, {255, 178,   0}, {255, 179,   0},
+    {255, 181,   0}, {255, 181,   0}, {255, 182,   0}, {255, 183,   0},
+    {255, 184,   0}, {255, 187,   7}, {255, 188,  10}, {255, 189,  14},
+    {255, 191,  18}, {255, 192,  21}, {255, 193,  25}, {255, 195,  29},
+    {255, 197,  36}, {255, 198,  40}, {255, 200,  43}, {255, 202,  51},
+    {255, 204,  54}, {255, 206,  61}, {255, 207,  65}, {255, 210,  72},
+    {255, 211,  76}, {255, 214,  83}, {255, 216,  91}, {255, 219,  98},
+    {255, 221, 105}, {255, 223, 109}, {255, 225, 116}, {255, 228, 123},
+    {255, 232, 134}, {255, 234, 142}, {255, 237, 149}, {255, 239, 156},
+    {255, 240, 160}, {255, 243, 167}, {255, 246, 174}, {255, 248, 182},
+    {255, 249, 185}, {255, 252, 193}, {255, 253, 196}, {255, 255, 204},
+    {255, 255, 207}, {255, 255, 211}, {255, 255, 218}, {255, 255, 222},
+    {255, 255, 225}, {255, 255, 229}, {255, 255, 233}, {255, 255, 236},
+    {255, 255, 240}, {255, 255, 244}, {255, 255, 247}, {255, 255, 255}
+};
+// *INDENT-ON*
+/* Save the two dimensional array as a png image
+ * Arguments:
+ *   double *data is a pointer to an array of nx*ny values
+ *   int nx is the number of COLUMNS to be written
+ *   int ny is the number of ROWS to be written
+ *   char *fname is the name of the picture
+ *   char lang is either 'c' or 'f' denoting the memory
+ *             layout. That is, if 'f' is given, then rows
+ *             and columns are swapped.
+ */
+int save_png(int *data, const int height, const int width, const char *fname)
+{
+    FILE *fp;
+    png_structp pngstruct_ptr = NULL;
+    png_infop pnginfo_ptr = NULL;
+    png_byte **row_pointers = NULL;
+    int i, j, max_value;
+
+    /* Determine the maximum value */
+    max_value = INT_MIN;
+    for (i = 0; i < height * width; i++) {
+        max_value = data[i] > max_value ? data[i] : max_value;
+    }
+
+    /* Default return status is failure */
+    int status = -1;
+
+    int pixel_size = 3;
+    int depth = 8;
+
+    fp = fopen(fname, "wb");
+    if (fp == NULL) {
+        goto fopen_failed;
+    }
+
+    pngstruct_ptr =
+        png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+
+    if (pngstruct_ptr == NULL) {
+        goto pngstruct_create_failed;
+    }
+
+    pnginfo_ptr = png_create_info_struct(pngstruct_ptr);
+
+    if (pnginfo_ptr == NULL) {
+        goto pnginfo_create_failed;
+    }
+
+    if (setjmp(png_jmpbuf(pngstruct_ptr))) {
+        goto setjmp_failed;
+    }
+
+    png_set_IHDR(pngstruct_ptr, pnginfo_ptr, (size_t) width,
+                 (size_t) height, depth, PNG_COLOR_TYPE_RGB,
+                 PNG_INTERLACE_NONE,
+                 PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+
+    row_pointers = png_malloc(pngstruct_ptr, height * sizeof(png_byte *));
+
+    for (i = 0; i < height; i++) {
+        png_byte *row = png_malloc(pngstruct_ptr,
+                                   sizeof(uint8_t) * width * pixel_size);
+        row_pointers[i] = row;
+
+        // Branch according to the memory layout
+        for (j = 0; j < width; j++) {
+            pixel_t pixel;
+            // Scale the values so that values between
+            // 0 and 100 degrees are mapped to values
+            // between 0 and 255
+            cmap(data[j + i * width], max_value, &pixel);
+            *row++ = pixel.red;
+            *row++ = pixel.green;
+            *row++ = pixel.blue;
+        }
+    }
+
+    png_init_io(pngstruct_ptr, fp);
+    png_set_rows(pngstruct_ptr, pnginfo_ptr, row_pointers);
+    png_write_png(pngstruct_ptr, pnginfo_ptr,
+                  PNG_TRANSFORM_IDENTITY, NULL);
+
+    status = 0;
+
+    for (i = 0; i < height; i++) {
+        png_free(pngstruct_ptr, row_pointers[i]);
+    }
+    png_free(pngstruct_ptr, row_pointers);
+
+setjmp_failed:
+pnginfo_create_failed:
+    png_destroy_write_struct(&pngstruct_ptr, &pnginfo_ptr);
+pngstruct_create_failed:
+    fclose(fp);
+fopen_failed:
+    return status;
+}
+
+/* This routine sets the RGB values for the pixel_t structure using
+ * the colormap data heat_colormap. If the value is outside the
+ * acceptable png values 0,255 blue or red color is used instead. */
+void cmap(const int value, const int max_value, pixel_t *pix)
+{
+    int ival;
+
+    ival = ((double) value / (double) max_value) * 255;
+    /* Check for wrong values */
+    if (ival > 255) {
+        ival = 255;
+    }
+    if (ival < 0) {
+        ival = 0;
+    }
+
+    /* Pick color from heat_colormap */
+    pix->red = heat_colormap[ival][0];
+    pix->green = heat_colormap[ival][1];
+    pix->blue = heat_colormap[ival][2];
+}

tasks/fortran/pngwriter.h

+#ifndef PNGWRITER_H_
+#define PNGWRITER_H_
+
+int save_png(int *data, const int nx, const int ny, const char *fname);
+
+#endif

tasks/fortran/pngwriter_mod.F90

+! PNG writer
+module pngwriter
+
+contains
+
+  function save_png(data, nx, ny, fname) result(stat)
+
+    use, intrinsic :: ISO_C_BINDING
+    implicit none
+
+    integer, dimension(:,:), intent(in) :: data
+    integer, intent(in) :: nx, ny
+    character(len=*), intent(in) :: fname
+    integer :: stat
+
+    ! Interface for save_png C-function
+    interface
+       ! The C-function definition is
+       !   int save_png(double *data, const int nx, const int ny,
+       !                const char *fname)
+       function save_png_c(data, nx, ny, fname, order) &
+            & bind(C,name="save_png") result(stat)
+         use, intrinsic :: ISO_C_BINDING
+         implicit none
+         integer(kind=C_INT) :: data(*)
+         integer(kind=C_INT), value, intent(IN) :: nx, ny
+         character(kind=C_CHAR), intent(IN) :: fname(*)
+         character(kind=C_CHAR), value, intent(IN) :: order
+         integer(kind=C_INT) :: stat
+       end function save_png_c
+    end interface
+
+    stat = save_png_c(data, nx, ny, trim(fname) // C_NULL_CHAR, 'f')
+    if (stat /= 0) then
+       write(*,*) 'save_png returned error!'
+    end if
+
+  end function save_png
+
+end module pngwriter

tasks/fortran/solution/Makefile

+ifeq ($(COMP),)
+COMP=cray
+endif
+
+COMMONDIR=../../../heat/common
+
+ifeq ($(COMP),cray)
+FC=ftn
+CC=cc
+FCFLAGS=-O3
+CCFLAGS=-O3 -I/appl/opt/libpng/include -I$(COMMONDIR)
+LDFLAGS=-L/appl/opt/libpng/lib
+LIBS=-lpng -lz
+endif
+
+ifeq ($(COMP),gnu)
+FC=gfortran
+CC=gcc
+FCFLAGS=-fopenmp -O3 -Wall
+CCFLAGS=-fopenmp -O3 -Wall -I$(COMMONDIR)
+LDFLAGS=
+LIBS=-lpng
+endif
+
+EXE=mandelbrot
+OBJS=mandelbrot.o pngwriter_mod.o
+OBJS_PNG=pngwriter.o
+CORRECT_OBJS_PNG=$(COMMONDIR)/pngwriter.o
+
+all: $(EXE)
+
+$(COMMONDIR)/pngwriter.o: $(COMMONDIR)/pngwriter.c $(COMMONDIR)/pngwriter.h
+pngwriter_mod.o: pngwriter_mod.F90
+mandelbrot.o: mandelbrot.F90 pngwriter_mod.o
+
+pngwriter.o: pngwriter.c pngwriter.h
+	$(CC) $(CCFLAGS) -c -o $@ $<
+
+$(EXE): $(OBJS) $(OBJS_PNG)
+	$(FC) $(FCFLAGS) $(OBJS) $(OBJS_PNG) -o $@ $(LDFLAGS) $(LIBS)
+
+%.o: %.F90
+	$(FC) $(FCFLAGS) -c $< -o $@
+
+%.o: %.c
+	$(CC) $(CCFLAGS) -c $< -o $@
+
+.PHONY: clean
+clean:
+	-/bin/rm -f $(EXE) a.out *.o *.mod *.png *~ mandelbrot.png

tasks/fortran/solution/mandelbrot.F90

+program mandelbrot
+
+  use iso_fortran_env, only : REAL64
+  use pngwriter
+  use omp_lib
+
+  implicit none
+
+  integer, parameter :: dp = REAL64
+
+  integer, parameter :: max_iter_count = 512
+  integer, parameter :: max_depth = 6
+  integer, parameter :: min_size = 32
+  integer, parameter :: subdiv = 4
+
+  integer, parameter :: w = 2048
+  integer, parameter :: h = w
+
+  integer, pointer, dimension(:,:) :: iter_counts
+
+  integer :: stat
+
+  complex(dp) :: cmin, cmax
+  real(dp) :: t0, t1
+
+  allocate(iter_counts(w, h))
+
+
+  t0 = omp_get_wtime()
+
+  cmin = (-1.5, -1.0)
+  cmax = (0.5, 1.0)
+
+  !$omp parallel
+  !$omp single
+  call mandelbrot_block(iter_counts, w, h, cmin, cmax, 0, 0, w, 1)
+  !$omp end single
+  !$omp end parallel
+
+  t1 = omp_get_wtime()
+
+  stat = save_png(iter_counts, h, w, 'mandelbrot.png')
+
+  deallocate(iter_counts)
+
+  write(*,*) 'Mandelbrot set computed in', t1 - t0, 's'
+
+contains
+
+  ! Computes the Mandelbrot image recursively
+  ! At each call, the image is divided into smaller blocks (by a factor of
+  ! subdiv), and the function is called recursively with arguments corresponding
+  ! to subblock. When maximum recursion depth is reached or size of block
+  ! is smaller than predefined minimum, one starts to calculate actual pixel
+  ! values
+  !
+  ! - - - - - - - -           -----       -----
+  ! |             |           |   |       |   |
+  ! |             |           -----       -----
+  ! |             |  -->                         -->   ...
+  ! |             |           -----       -----
+  ! |             |           |   |       |   |
+  ! |             |           -----       -----
+  ! ---------------
+  !
+
+  recursive subroutine mandelbrot_block(iter_counts, w, h, cmin, cmax, &
+       x0, y0, d, depth)
+    implicit none
+
+    integer, pointer, dimension(:,:), intent(inout) :: iter_counts
+    integer, intent(in) :: w, h, x0, y0, d, depth
+    complex(dp), intent(in) :: cmin, cmax
+
+    integer :: block_size, i, j
+
+    block_size = d / subdiv
+    if ((depth + 1 < max_depth) .and. (block_size > min_size)) then
+       ! Subdivide recursively
+       do i=0, subdiv - 1
+          do j=0, subdiv - 1
+             !$omp task
+             call mandelbrot_block(iter_counts, w, h, cmin, cmax, &
+                  x0 + i*block_size, y0 + j*block_size, &
+                  block_size, depth + 1)
+             !$omp end task
+          end do
+       end do
+    else
+       ! Last recursion level reached, calculate the values
+       do j = y0 + 1, y0 + d
+          do i = x0 + 1, x0 + d
+             iter_counts(i, j) = kernel(h, w, cmin, cmax, i-1, j-1)
+          end do
+       end do
+    end if
+
+  end subroutine mandelbrot_block
+
+  ! Calculate iteration count for a pixel.
+  ! This function does not need to be edited
+  integer function kernel(h, w, cmin, cmax, x, y)
+    implicit none
+
+    integer, intent(in) :: h, w, x, y
+    complex(dp), intent(in) :: cmin, cmax
+
+    integer :: iteration
+    complex(dp) :: z, dc, c
+    real(dp) :: fx, fy
+
+    dc = cmax - cmin
+    fx = real(x, dp) / w
+    fy = real(y, dp) / h
+    c = cmin + fx * real(dc) +  fy * imag(dc) * (0.0, 1.0)
+
+    z = c
+    iteration = 0
+
+    do while (iteration < max_iter_count .and. abs(z)**2 < 4)
+       z = z**2 + c
+       iteration = iteration + 1
+    end do
+
+    kernel = iteration
+
+  end function kernel
+
+end program mandelbrot

tasks/fortran/solution/pngwriter.c

+#include <png.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <limits.h>
+#include "pngwriter.h"
+
+/* Datatype for RGB pixel */
+typedef struct {
+    uint8_t red;
+    uint8_t green;
+    uint8_t blue;
+} pixel_t;
+
+void cmap(const int value, const int max_val, pixel_t *pix);
+
+/* Heat colormap from black to white */
+// *INDENT-OFF*
+static int heat_colormap[256][3] = {
+    {  0,   0,   0}, { 35,   0,   0}, { 52,   0,   0}, { 60,   0,   0},
+    { 63,   1,   0}, { 64,   2,   0}, { 68,   5,   0}, { 69,   6,   0},
+    { 72,   8,   0}, { 74,  10,   0}, { 77,  12,   0}, { 78,  14,   0},
+    { 81,  16,   0}, { 83,  17,   0}, { 85,  19,   0}, { 86,  20,   0},
+    { 89,  22,   0}, { 91,  24,   0}, { 92,  25,   0}, { 94,  26,   0},
+    { 95,  28,   0}, { 98,  30,   0}, {100,  31,   0}, {102,  33,   0},
+    {103,  34,   0}, {105,  35,   0}, {106,  36,   0}, {108,  38,   0},
+    {109,  39,   0}, {111,  40,   0}, {112,  42,   0}, {114,  43,   0},
+    {115,  44,   0}, {117,  45,   0}, {119,  47,   0}, {119,  47,   0},
+    {120,  48,   0}, {122,  49,   0}, {123,  51,   0}, {125,  52,   0},
+    {125,  52,   0}, {126,  53,   0}, {128,  54,   0}, {129,  56,   0},
+    {129,  56,   0}, {131,  57,   0}, {132,  58,   0}, {134,  59,   0},
+    {134,  59,   0}, {136,  61,   0}, {137,  62,   0}, {137,  62,   0},
+    {139,  63,   0}, {139,  63,   0}, {140,  65,   0}, {142,  66,   0},  
+    {148,  71,   0}, {149,  72,   0}, {149,  72,   0}, {151,  73,   0},
+    {151,  73,   0}, {153,  75,   0}, {153,  75,   0}, {154,  76,   0},
+    {154,  76,   0}, {154,  76,   0}, {156,  77,   0}, {156,  77,   0},
+    {157,  79,   0}, {157,  79,   0}, {159,  80,   0}, {159,  80,   0},
+    {159,  80,   0}, {160,  81,   0}, {160,  81,   0}, {162,  82,   0},
+    {162,  82,   0}, {163,  84,   0}, {163,  84,   0}, {165,  85,   0},
+    {165,  85,   0}, {166,  86,   0}, {166,  86,   0}, {166,  86,   0},
+    {168,  87,   0}, {168,  87,   0}, {170,  89,   0}, {170,  89,   0},
+    {171,  90,   0}, {171,  90,   0}, {173,  91,   0}, {173,  91,   0},
+    {174,  93,   0}, {174,  93,   0}, {176,  94,   0}, {176,  94,   0},
+    {177,  95,   0}, {177,  95,   0}, {179,  96,   0}, {179,  96,   0},
+    {180,  98,   0}, {182,  99,   0}, {182,  99,   0}, {183, 100,   0},
+    {183, 100,   0}, {185, 102,   0}, {185, 102,   0}, {187, 103,   0},
+    {187, 103,   0}, {188, 104,   0}, {188, 104,   0}, {190, 105,   0},
+    {191, 107,   0}, {191, 107,   0}, {193, 108,   0}, {193, 108,   0},
+    {194, 109,   0}, {196, 110,   0}, {196, 110,   0}, {197, 112,   0},
+    {197, 112,   0}, {199, 113,   0}, {200, 114,   0}, {200, 114,   0},
+    {202, 116,   0}, {202, 116,   0}, {204, 117,   0}, {205, 118,   0},
+    {205, 118,   0}, {207, 119,   0}, {208, 121,   0}, {208, 121,   0},
+    {210, 122,   0}, {211, 123,   0}, {211, 123,   0}, {213, 124,   0},
+    {214, 126,   0}, {214, 126,   0}, {216, 127,   0}, {217, 128,   0},
+    {217, 128,   0}, {219, 130,   0}, {221, 131,   0}, {221, 131,   0},
+    {222, 132,   0}, {224, 133,   0}, {224, 133,   0}, {225, 135,   0},
+    {227, 136,   0}, {227, 136,   0}, {228, 137,   0}, {230, 138,   0},
+    {230, 138,   0}, {231, 140,   0}, {233, 141,   0}, {233, 141,   0},
+    {234, 142,   0}, {236, 144,   0}, {236, 144,   0}, {238, 145,   0},
+    {239, 146,   0}, {241, 147,   0}, {241, 147,   0}, {242, 149,   0},
+    {244, 150,   0}, {244, 150,   0}, {245, 151,   0}, {247, 153,   0},
+    {247, 153,   0}, {248, 154,   0}, {250, 155,   0}, {251, 156,   0},
+    {251, 156,   0}, {253, 158,   0}, {255, 159,   0}, {255, 159,   0},
+    {255, 160,   0}, {255, 161,   0}, {255, 163,   0}, {255, 163,   0},
+    {255, 164,   0}, {255, 165,   0}, {255, 167,   0}, {255, 167,   0},
+    {255, 168,   0}, {255, 169,   0}, {255, 169,   0}, {255, 170,   0},
+    {255, 172,   0}, {255, 173,   0}, {255, 173,   0}, {255, 174,   0},
+    {255, 175,   0}, {255, 177,   0}, {255, 178,   0}, {255, 179,   0},
+    {255, 181,   0}, {255, 181,   0}, {255, 182,   0}, {255, 183,   0},
+    {255, 184,   0}, {255, 187,   7}, {255, 188,  10}, {255, 189,  14},
+    {255, 191,  18}, {255, 192,  21}, {255, 193,  25}, {255, 195,  29},
+    {255, 197,  36}, {255, 198,  40}, {255, 200,  43}, {255, 202,  51},
+    {255, 204,  54}, {255, 206,  61}, {255, 207,  65}, {255, 210,  72},
+    {255, 211,  76}, {255, 214,  83}, {255, 216,  91}, {255, 219,  98},
+    {255, 221, 105}, {255, 223, 109}, {255, 225, 116}, {255, 228, 123},
+    {255, 232, 134}, {255, 234, 142}, {255, 237, 149}, {255, 239, 156},
+    {255, 240, 160}, {255, 243, 167}, {255, 246, 174}, {255, 248, 182},
+    {255, 249, 185}, {255, 252, 193}, {255, 253, 196}, {255, 255, 204},
+    {255, 255, 207}, {255, 255, 211}, {255, 255, 218}, {255, 255, 222},
+    {255, 255, 225}, {255, 255, 229}, {255, 255, 233}, {255, 255, 236},
+    {255, 255, 240}, {255, 255, 244}, {255, 255, 247}, {255, 255, 255}
+};
+// *INDENT-ON*
+/* Save the two dimensional array as a png image
+ * Arguments:
+ *   double *data is a pointer to an array of nx*ny values
+ *   int nx is the number of COLUMNS to be written
+ *   int ny is the number of ROWS to be written
+ *   char *fname is the name of the picture
+ *   char lang is either 'c' or 'f' denoting the memory
+ *             layout. That is, if 'f' is given, then rows
+ *             and columns are swapped.
+ */
+int save_png(int *data, const int height, const int width, const char *fname)
+{
+    FILE *fp;
+    png_structp pngstruct_ptr = NULL;
+    png_infop pnginfo_ptr = NULL;
+    png_byte **row_pointers = NULL;
+    int i, j, max_value;
+
+    /* Determine the maximum value */
+    max_value = INT_MIN;
+    for (i = 0; i < height * width; i++) {
+        max_value = data[i] > max_value ? data[i] : max_value;
+    }
+
+    /* Default return status is failure */
+    int status = -1;
+
+    int pixel_size = 3;
+    int depth = 8;
+
+    fp = fopen(fname, "wb");
+    if (fp == NULL) {
+        goto fopen_failed;
+    }
+
+    pngstruct_ptr =
+        png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+
+    if (pngstruct_ptr == NULL) {
+        goto pngstruct_create_failed;
+    }
+
+    pnginfo_ptr = png_create_info_struct(pngstruct_ptr);
+
+    if (pnginfo_ptr == NULL) {
+        goto pnginfo_create_failed;
+    }
+
+    if (setjmp(png_jmpbuf(pngstruct_ptr))) {
+        goto setjmp_failed;
+    }
+
+    png_set_IHDR(pngstruct_ptr, pnginfo_ptr, (size_t) width,
+                 (size_t) height, depth, PNG_COLOR_TYPE_RGB,
+                 PNG_INTERLACE_NONE,
+                 PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+
+    row_pointers = png_malloc(pngstruct_ptr, height * sizeof(png_byte *));
+
+    for (i = 0; i < height; i++) {
+        png_byte *row = png_malloc(pngstruct_ptr,
+                                   sizeof(uint8_t) * width * pixel_size);
+        row_pointers[i] = row;
+
+        // Branch according to the memory layout
+        for (j = 0; j < width; j++) {
+            pixel_t pixel;
+            // Scale the values so that values between
+            // 0 and 100 degrees are mapped to values
+            // between 0 and 255
+            cmap(data[j + i * width], max_value, &pixel);
+            *row++ = pixel.red;
+            *row++ = pixel.green;
+            *row++ = pixel.blue;
+        }
+    }
+
+    png_init_io(pngstruct_ptr, fp);
+    png_set_rows(pngstruct_ptr, pnginfo_ptr, row_pointers);
+    png_write_png(pngstruct_ptr, pnginfo_ptr,
+                  PNG_TRANSFORM_IDENTITY, NULL);
+
+    status = 0;
+
+    for (i = 0; i < height; i++) {
+        png_free(pngstruct_ptr, row_pointers[i]);
+    }
+    png_free(pngstruct_ptr, row_pointers);
+
+setjmp_failed:
+pnginfo_create_failed:
+    png_destroy_write_struct(&pngstruct_ptr, &pnginfo_ptr);
+pngstruct_create_failed:
+    fclose(fp);
+fopen_failed:
+    return status;
+}
+
+/* This routine sets the RGB values for the pixel_t structure using
+ * the colormap data heat_colormap. If the value is outside the
+ * acceptable png values 0,255 blue or red color is used instead. */
+void cmap(const int value, const int max_value, pixel_t *pix)
+{
+    int ival;
+
+    ival = ((double) value / (double) max_value) * 255;
+    /* Check for wrong values */
+    if (ival > 255) {
+        ival = 255;
+    }
+    if (ival < 0) {
+        ival = 0;
+    }
+
+    /* Pick color from heat_colormap */
+    pix->red = heat_colormap[ival][0];
+    pix->green = heat_colormap[ival][1];
+    pix->blue = heat_colormap[ival][2];
+}

tasks/fortran/solution/pngwriter.h

+#ifndef PNGWRITER_H_
+#define PNGWRITER_H_
+
+int save_png(int *data, const int nx, const int ny, const char *fname);
+
+#endif

tasks/fortran/solution/pngwriter_mod.F90

+! PNG writer
+module pngwriter
+
+contains
+
+  function save_png(data, nx, ny, fname) result(stat)
+
+    use, intrinsic :: ISO_C_BINDING
+    implicit none
+
+    integer, dimension(:,:), intent(in) :: data
+    integer, intent(in) :: nx, ny
+    character(len=*), intent(in) :: fname
+    integer :: stat
+
+    ! Interface for save_png C-function
+    interface
+       ! The C-function definition is
+       !   int save_png(double *data, const int nx, const int ny,
+       !                const char *fname)
+       function save_png_c(data, nx, ny, fname, order) &
+            & bind(C,name="save_png") result(stat)
+         use, intrinsic :: ISO_C_BINDING
+         implicit none
+         integer(kind=C_INT) :: data(*)
+         integer(kind=C_INT), value, intent(IN) :: nx, ny
+         character(kind=C_CHAR), intent(IN) :: fname(*)
+         character(kind=C_CHAR), value, intent(IN) :: order
+         integer(kind=C_INT) :: stat
+       end function save_png_c
+    end interface
+
+    stat = save_png_c(data, nx, ny, trim(fname) // C_NULL_CHAR, 'f')
+    if (stat /= 0) then
+       write(*,*) 'save_png returned error!'
+    end if
+
+  end function save_png
+
+end module pngwriter

work-sharing/LICENSE.txt

+
+Copyright (C) 2018  CSC - IT Center for Science Ltd.
+
+Licensed under the terms of the GNU General Public License as
+published by the Free Software Foundation, either version 3 of the
+License, or (at your option) any later version.
+
+Code is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+Copy of the GNU General Public License can be obtained from
+<http://www.gnu.org/licenses/>.
+

work-sharing/README.md

+## Work sharing for a simple loop ##
+
+In [sum.c](sum.c) (or [sum.F90](sum.F90) for Fortran) is a skeleton
+implementation for a simple summation of two vectors, C=A+B. Add the
+computation loop and add the parallel region with work sharing
+directives so that the vector addition is executed in parallel.

work-sharing/solution/sum.F90

+program vectorsum
+  implicit none
+  integer, parameter :: rk = kind(1d0)
+  integer, parameter :: ik = selected_int_kind(9)
+  integer, parameter :: nx = 102400
+
+  real(kind=rk), dimension(nx) :: vecA, vecB, vecC
+  real(kind=rk)    :: sum
+  integer(kind=ik) :: i
+
+  ! Initialization of vectors
+  do i = 1, nx
+     vecA(i) = 1.0_rk/(real(nx - i + 1, kind=rk))
+     vecB(i) = vecA(i)**2
+  end do
+
+  !$omp parallel do default(shared) private(i)
+  do i = 1, nx
+     vecC(i) = vecA(i) * vecB(i)
+  end do
+  !$omp end parallel do
+
+  ! Compute the check value
+  write(*,*) 'Reduction sum: ', sum(vecC)
+
+end program vectorsum

work-sharing/solution/sum.c

+#include <stdio.h>
+
+#define NX 102400
+
+int main(void)
+{
+    double vecA[NX], vecB[NX], vecC[NX];
+    double sum;
+    int i;
+
+    /* Initialization of the vectors */
+    for (i = 0; i < NX; i++) {
+        vecA[i] = 1.0 / ((double)(NX - i));
+        vecB[i] = vecA[i] * vecA[i];
+    }
+
+    #pragma omp parallel for default(shared) private(i)
+    for (i = 0; i < NX; i++) {
+        vecC[i] = vecA[i] * vecB[i];
+    }
+
+    sum = 0.0;
+    /* Compute the check value */
+    for (i = 0; i < NX; i++) {
+        sum += vecC[i];
+    }
+    printf("Reduction sum: %18.16f\n", sum);
+
+    return 0;
+}