multithread C++ solver

python/frep-C.py

 #!/usr/bin/env python3
 #
 # frep-C.py
-#    functional representation to C++ solver
+#    functional representation to C++ multi-threaded solver
 #
 # usage:
 #    pcb.py | frep-C.py [dpi [filename]]
 #
-# Neil Gershenfeld 8/23/21
-# (c) Massachusetts Institute of Technology 2021
+# Neil Gershenfeld 11/26/22
 #
 # This work may be reproduced, modified, distributed,
 # performed, and displayed for any purpose, but must
@@ -71,12 +70,11 @@ file.write(
 f"""
 #include <iostream>
 #include <cmath>
+#include <thread>
 #include <png.h>
+//
 using namespace std;
-int fn(float X,float Y,float Z) {{
-   return ({fn});
-   }}
-int main(int argc, char** argv) {{
+//
 float xmin = {xmin};
 float xmax = {xmax};
 float ymin = {ymin};
@@ -90,15 +88,25 @@ int main(int argc, char** argv) {{
 int ny = (ymax-ymin)/delta;
 int *m = (int*) calloc(nx*ny,sizeof(int));
 float layers[] = {{{layers}}};
+int nthreads = std::thread::hardware_concurrency();
+//
+int fn(float X,float Y,float Z) {{
+   return ({fn});
+   }}
+//
+void calc(int nx,int ny,int nthreads,int thread) {{
    int intensity;
-   for (int layer = 0; layer < {nlayers}; ++layer) {{
+   for (int layer = 0; layer < 2; ++layer) {{
       float z = layers[layer];
+      if (thread == 0)
          cout << "   z = " << z << endl;
       if (zmin == zmax)
-         intensity = 0xFFFFFF;
+         intensity = 0xffffff;
       else
          intensity = ((int) (255*(z-zmin)/(zmax-zmin))) | (255 << 8) | (255 << 16);
-      for (int iy = 0; iy < ny; ++iy) {{
+      int iystart = thread*ny/nthreads;
+      int iyend = (thread+1)*ny/nthreads;
+      for (int iy = iystart; iy < iyend; ++iy) {{
          float y = ymin+iy*delta;
          for (int ix = 0; ix < nx; ++ix) {{
             float x = xmin+ix*delta;
@@ -106,6 +114,16 @@ int main(int argc, char** argv) {{
             }}
          }}
       }}
+   }}
+//
+int main(int argc, char** argv) {{
+   cout << "   calculate " << nx << "x" << ny << " with " << nthreads << " threads" << endl;
+   std::thread threads[nthreads];
+   for (int i = 0; i < nthreads; ++i)
+      threads[i] = std::thread(calc,nx,ny,nthreads,i);
+   for (int i = 0; i < nthreads; ++i)
+      threads[i].join();
+   //
    FILE *file;
    file = fopen("{filename}","wb");
 	png_structp pngfile;
@@ -138,7 +156,7 @@ file.close()
 # compile
 #
 print("compile ...")
-os.system("time g++ frep-C.cpp -o frep-C -lm -lpng -O -ffast-math")
+os.system("time g++ frep-C.cpp -o frep-C -lm -lpng -O -ffast-math -pthread")
 #
 # execute
 #
@@ -147,4 +165,4 @@ os.system("time ./frep-C")
 #
 # clean up
 #
-os.system("rm -f frep-C.cpp frep-C")
+#os.system("rm -f frep-C.cpp frep-C")