Thread local storage: add example

From 50f91c868e803b0519be2bae1862d1742dec9afa Mon Sep 17 00:00:00 2001
From: Heinrich Schuchardt <[email protected]>
Date: Tue, 10 Jan 2017 22:37:54 +0100
Subject: [PATCH 1/1] Add multithreading example

Signed-off-by: Heinrich Schuchardt <[email protected]>
---
 examples/tls/Makefile       |   5 +
 examples/tls/clustering.mod | 109 +++++++++++++++++++++
 examples/tls/multiseed.c    | 226 ++++++++++++++++++++++++++++++++++++++++++++
 examples/tls/thread.h       |  53 +++++++++++
 4 files changed, 393 insertions(+)
 create mode 100644 examples/tls/Makefile
 create mode 100644 examples/tls/clustering.mod
 create mode 100644 examples/tls/multiseed.c
 create mode 100644 examples/tls/thread.h

diff --git a/examples/tls/Makefile b/examples/tls/Makefile
new file mode 100644
index 0000000..61063d3
--- /dev/null
+++ b/examples/tls/Makefile
@@ -0,0 +1,5 @@
+all:
+	gcc multiseed.c -I. -lglpk -pthread -o multiseed
+
+check:
+	./multiseed clustering.mod 20
diff --git a/examples/tls/clustering.mod b/examples/tls/clustering.mod
new file mode 100644
index 0000000..21c17b0
--- /dev/null
+++ b/examples/tls/clustering.mod
@@ -0,0 +1,109 @@
+/*
+ * Author: Heinrich Schuchardt <[email protected]>
+ *
+ * This model solves a clustering problem:
+ *
+ * Out of 50 towns select 3 to be cluster centers and assign the other
+ * towns to the clusters such that the sum of the population weighted
+ * euclidian distances between towns and centers is minimized.
+ *
+ * The solution is saved as a scalable vector graphic file with a
+ * pseudo-random file name.
+ */
+
+# Output file
+param fn, symbolic := "00000" & 100000 * Uniform01();
+param f, symbolic := "ct" & substr(fn, length(fn) - 4) & ".svg";
+
+# Centers
+param nc := 3;
+set C := {1 .. nc};
+
+# Towns
+param nt := 50;
+set T := {1 .. nt};
+param xt{T} := Uniform01();
+param yt{T} := Uniform01();
+param pt{T} := ceil(1000 * Uniform01());
+
+# Image size
+param scale := 1000;
+
+# Colors
+# saturation [0, 255]
+param sat := 192;
+param hue{c in C} := 6 * (c - 1) / nc;
+param red{c in C} := 
+  if hue[c] <= 1 or hue[c] >= 5 then 255
+  else (if hue[c] >=2 and hue[c] <= 4 then 255 - sat
+  else (if hue[c] <=2 then 255 - sat + sat * (2-hue[c])
+  else 255 - sat + sat * (hue[c]-4) ));
+param green{c in C} := 
+  if hue[c] >= 1 and hue[c] <= 3 then 255
+  else (if hue[c] >= 4 then 255 - sat
+  else (if hue[c] <=1 then 255 - sat + sat * hue[c]
+  else 255 - sat + sat * (4-hue[c]) ));
+param blue{c in C} := 
+  if hue[c] >= 3 and hue[c] <= 5 then 255
+  else (if hue[c] <=2 then 255 - sat
+  else (if hue[c] <=3 then 255 - sat + sat * (hue[c]-2)
+  else 255 - sat + sat * (6-hue[c]) ));
+
+var x{T};
+var y{T,T}, binary;
+
+minimize obj : sum{c in T, t in T} y[c,t] * pt[t]
+               * sqrt((xt[c] - xt[t])^2 + (yt[c] - yt[t])^2);
+
+s.t. sumx : sum{c in T} x[c] = nc;
+s.t. cxy{c in T, t in T} : y[c,t] <= x[c];
+s.t. sumy{t in T} : sum{c in T} y[c,t] = 1;
+
+solve;
+
+for {c in T : x[c] > .5} {
+  printf "Center %5.4f %5.4f\n", xt[c], yt[c];
+  for {t in T : y[c,t] > .5} {
+    printf "  Town %5.4f %5.4f (%5.0f)\n", xt[t], yt[t], pt[t];
+  }
+}
+
+# Output the solution as scalable vector graphic
+
+# header
+printf "<?xml version=""1.0"" standalone=""no""?>\n" > f;
+printf "<!DOCTYPE svg PUBLIC ""-//W3C//DTD SVG 1.1//EN"" \n" >> f;
+printf """http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"">\n" >> f;
+printf "<svg width=""%d"" height=""%d"" version=""1.0"" \n",
+  1.2 * scale, 1.2 * scale >> f;
+printf "xmlns=""http://www.w3.org/2000/svg"">\n" >> f;
+
+# background
+printf "<rect x=""0"" y=""0"" width=""%d"" height=""%d""" &
+       " stroke=""none"" fill=""white""/>\n",
+       1.2 * scale, 1.2 * scale>> f;
+
+# border
+printf "<rect x=""%d"" y=""%d"" width=""%d"" height=""%d""" &
+       " stroke=""black"" stroke-width="".5"" fill=""white""/>\n",
+       .1 * scale, .1 * scale, scale, scale >> f;
+
+# circles for towns
+for {t in T}
+  printf {s in T, c in C : y[s,t] > .5 
+         && c = floor( .5 + sum{u in T : u <= s} x[u])}
+         "<circle cx=""%f"" cy=""%f"" r=""%f"" stroke=""black"" " &
+         "stroke-width=""1"" fill=""rgb(%d,%d,%d)""/>\n",
+         (.1 + xt[t]) * scale, (.1 + yt[t]) * scale, .001 * sqrt(pt[t]) * scale,
+         red[c], green[c] , blue[c] >> f;
+
+# lines from towns to assigned centers
+for {t in T, c in T : y[c,t] > .5}
+  printf "<line x1=""%f"" y1=""%f"" x2=""%f"" y2=""%f""" &
+         " style=""stroke:black;stroke-width:.5""/>\n",
+         (.1 + xt[c]) * scale, (.1 + yt[c]) * scale,
+         (.1 + xt[t]) * scale, (.1 + yt[t]) * scale >> f;
+
+printf "</svg>\n" >> f;
+
+end;
diff --git a/examples/tls/multiseed.c b/examples/tls/multiseed.c
new file mode 100644
index 0000000..68a4344
--- /dev/null
+++ b/examples/tls/multiseed.c
@@ -0,0 +1,226 @@
+/* multiseed.d (multithreading demo) */
+
+/***********************************************************************
+*
+*  This code is part of GLPK (GNU Linear Programming Kit).
+*
+*  Author: Heinrich Schuchardt <[email protected]>
+*
+*  Copyright (C) 2000-2017 Andrew Makhorin, Department for Applied
+*  Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
+*  reserved. E-mail: <[email protected]>.
+*
+*  GLPK is free software: you can redistribute it and/or modify it
+*  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.
+*
+*  GLPK 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.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+
+/*
+ * This program demonstrates running the GLPK library with multiple threads.
+ *
+ * When called the program requires two arguments:
+ *
+ * filename - the name of the MathProg model to be solved
+ * threads  - the count of parallel threads to be run.
+ *
+ * Each thread is run with a different seed for the random number generator
+ * provided by the GLPK library.
+ */
+
+#include <glpk.h>
+#include <malloc.h>
+#include <setjmp.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "thread.h"
+
+#define BUFLEN 256
+
+struct task {
+   pthread_t tid;
+   char *filename;
+   int seed;
+   size_t pos;
+   char buf[BUFLEN + 1];
+   int line;
+   jmp_buf jmp;
+};
+
+pthread_mutex_t mutex;
+
+int term_hook(void *info, const char *text)
+{
+   struct task *task = (struct task *) info;
+   size_t len = strlen(text);
+   char *pos;
+
+   pthread_mutex_lock(&mutex);
+   if (task->pos + len > BUFLEN) {
+      printf("%02d-%05d %s%s", task->seed, ++task->line, task->buf, text);
+      task->pos = 0;
+      task->buf[0] = 0;
+   } else {
+      strcpy(task->buf + task->pos, text);
+      task->pos += len;
+   }
+   if (strchr(task->buf, '\n')) {
+      printf("%02d-%05d %s", task->seed, ++task->line, task->buf);
+      task->pos = 0;
+      task->buf[0] = 0;
+   }
+   pthread_mutex_unlock(&mutex);
+   return -1;
+}
+
+void error_hook(void *info)
+{
+   struct task *task = (struct task *) info;
+
+   term_hook(task, "Error caught\n");
+   glp_free_env();
+   longjmp(task->jmp, 1);
+}
+
+void worker(void *arg)
+{
+   struct task *task = (struct task *) arg;
+   int ret;
+   glp_prob *lp;
+   glp_tran *tran;
+   glp_iocp iocp;
+
+   glp_error_hook(error_hook, task);
+
+   if (setjmp(task->jmp)) {
+      return;
+   }
+
+   glp_term_hook(term_hook, arg);
+
+   glp_printf("Seed %02d\n", task->seed);
+
+   lp = glp_create_prob();
+   tran = glp_mpl_alloc_wksp();
+   glp_mpl_init_rand(tran, task->seed);
+
+   ret = glp_mpl_read_model (tran, task->filename, GLP_OFF);
+   if (ret != 0) {
+      glp_mpl_free_wksp (tran);
+      glp_delete_prob(lp);
+      glp_printf("Model %s not valid\n", task->filename);
+      glp_free_env();
+   }
+
+   ret = glp_mpl_generate(tran, NULL);
+   if (ret != 0) {
+      glp_mpl_free_wksp (tran);
+      glp_delete_prob(lp);
+      glp_printf("Cannot generate model %s\n", task->filename);
+      glp_free_env();
+   }
+
+   glp_mpl_build_prob(tran, lp);
+
+   glp_init_iocp(&iocp);
+   iocp.presolve = GLP_ON;
+
+   ret = glp_intopt(lp, &iocp);
+
+   if (ret == 0) {
+      glp_mpl_postsolve(tran, lp, GLP_MIP);
+   }
+
+   glp_mpl_free_wksp (tran);
+   glp_delete_prob(lp);
+
+   if (0 == task->seed % 3) {
+      glp_error("Voluntarily throwing an error in %s at line %d\n",
+                __FILE__, __LINE__);
+   }
+
+   glp_term_hook(NULL, NULL);
+
+   glp_error_hook(NULL, NULL);
+
+   glp_free_env();
+}
+
+#ifdef __WOE__
+DWORD run(void *arg)
+{
+#else
+void *run(void *arg)
+{
+#endif
+   worker(arg);
+   pthread_exit(NULL);
+#ifdef __WOE__
+   return 0;
+#else
+   return NULL;
+#endif
+}
+
+int main(int argc, char *argv[])
+{
+   pthread_t t;
+   int i, n, rc;
+   struct task *tasks;
+
+   if (argc != 3) {
+      printf("Usage %s filename threads\n"
+             "  filename - MathProg model file\n"
+             "  threads  - number of threads\n",
+             argv[0]);
+      exit(EXIT_FAILURE);
+   }
+
+   n = atoi(argv[2]);
+   if (n > 50) {
+      printf("Number of threads is to high (> 50).\n");
+      exit(EXIT_FAILURE);
+   }
+   if (n <= 1) {
+      printf("Need positive number of threads\n");
+      exit(EXIT_FAILURE);
+   }
+
+   tasks = calloc(n, sizeof(struct task));
+   if (!tasks) {
+      printf("Out of memory");
+      exit(EXIT_FAILURE);
+   }
+
+   pthread_mutex_init(&mutex, NULL);
+
+   for (i = 0; i < n; ++i) {
+      tasks[i].filename = argv[1];
+      tasks[i].seed = i + 1;
+      tasks[i].pos = 0;
+      tasks[i].buf[0] = 0;
+      tasks[i].line = 0;
+      rc = pthread_create(&tasks[i].tid, NULL, run, &tasks[i]);
+      if (rc) {
+         printf("ERROR; return code from pthread_create() is %d\n", rc);
+         exit(EXIT_FAILURE);
+      }
+   }
+   for (i = 0; i < n; ++i) {
+      pthread_join(tasks[i].tid, NULL);
+   }
+
+   pthread_mutex_destroy(&mutex);
+
+   return EXIT_SUCCESS;
+}
diff --git a/examples/tls/thread.h b/examples/tls/thread.h
new file mode 100644
index 0000000..7853c36
--- /dev/null
+++ b/examples/tls/thread.h
@@ -0,0 +1,53 @@
+/* thread.h (Threading) */
+
+/***********************************************************************
+*  This code is part of GLPK (GNU Linear Programming Kit).
+*
+*  Author: Heinrich Schuchardt <[email protected]>
+*
+*  Copyright (C) 2000-2017 Andrew Makhorin, Department for Applied
+*  Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
+*  reserved. E-mail: <[email protected]>.
+*
+*  GLPK is free software: you can redistribute it and/or modify it
+*  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.
+*
+*  GLPK 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.
+*
+*  You should have received a copy of the GNU General Public License
+*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
+***********************************************************************/
+
+#ifndef THREAD_H
+
+#define THREAD_H 1
+
+#ifdef HAVE_CONF_H
+#include "config.h"
+#endif // HAVE_CONF_H
+
+#ifdef __WOE__
+#include <windows.h>
+typedef CRITICAL_SECTION pthread_mutex_t;
+typedef DWORD pthread_t;
+//@todo The return type of routine C is "DWORD" for Windows and "void *" for Posix.
+#define pthread_create(A,B,C,D) \
+  (int)(CreateThread(NULL,0,(LPTHREAD_START_ROUTINE)&C,D,0,A)==NULL)
+#define pthread_exit(A)          ExitThread(0)
+#define pthread_mutex_destroy(A) DeleteCriticalSection(A)
+#define pthread_mutex_init(A,B)  (InitializeCriticalSection(A),0)
+#define pthread_mutex_lock(A)    (EnterCriticalSection(A),0)
+#define pthread_mutex_unlock(A)  (LeaveCriticalSection(A),0)
+#define pthread_self()           GetCurrentThreadId()
+#define pthread_join(A, B) \
+  (WaitForSingleObject(A, INFINITE),CloseHandle(A),0)
+#else
+#include <pthread.h>
+#endif
+
+#endif // THREAD_H
-- 
2.1.4