Multithreading gradient descent, required FFT to be of preset sizes, clean up timing, fix Windows bug

fast_tsne.m

@@ -245,6 +245,12 @@
             load_affinities = 0;
         end
     end
+
+    if (~isfield(opts, 'nthreads'))
+        nthreads = 0;
+    else
+        nthreads = opts.nthreads;
+    end
 
     X = double(X);
 
@@ -264,7 +270,9 @@
 
     disp('Data written');
     tic
-    [flag, cmdout] = system(fullfile(tsne_path,'/fast_tsne'), '-echo');
+    %[flag, cmdout] = system(fullfile(tsne_path,'/fast_tsne'), '-echo');
+    cmd = sprintf('%s temp/data.dat temp/result.dat %d',fullfile(tsne_path,'/fast_tsne'), nthreads);
+    [flag, cmdout] = system(cmd, '-echo');
     if(flag~=0)
         error(cmdout);
     end

src/nbodyfft.cpp

@@ -1,10 +1,10 @@
 #include "winlibs/stdafx.h"
 
+#include "parallel_for.h"
+#include "time_code.h"
 #include "nbodyfft.h"
 
 
-using namespace std;
-
 
 void precompute_2d(double x_max, double x_min, double y_max, double y_min, int n_boxes, int n_interpolation_points,
                    kernel_type_2d kernel, double *box_lower_bounds, double *box_upper_bounds, double *y_tilde_spacings,
@@ -72,7 +72,7 @@ void precompute_2d(double x_max, double x_min, double y_max, double y_min, int n
 
 void n_body_fft_2d(int N, int n_terms, double *xs, double *ys, double *chargesQij, int n_boxes,
                    int n_interpolation_points, double *box_lower_bounds, double *box_upper_bounds,
-                   double *y_tilde_spacings, complex<double> *fft_kernel_tilde, double *potentialQij) {
+                   double *y_tilde_spacings, complex<double> *fft_kernel_tilde, double *potentialQij, unsigned int nthreads) {
     int n_total_boxes = n_boxes * n_boxes;
     int total_interpolation_points = n_total_boxes * n_interpolation_points * n_interpolation_points;
 
@@ -113,6 +113,9 @@ void n_body_fft_2d(int N, int n_terms, double *xs, double *ys, double *chargesQi
         y_in_box[i] = (ys[i] - y_min) / box_width;
     }
 
+    INITIALIZE_TIME
+    START_TIME
+
     /*
      * Step 1: Interpolate kernel using Lagrange polynomials and compute the w coefficients
      */
@@ -143,6 +146,8 @@ void n_body_fft_2d(int N, int n_terms, double *xs, double *ys, double *chargesQi
         }
     }
 
+        END_TIME("Step 1");
+        START_TIME;
     /*
      * Step 2: Compute the values v_{m, n} at the equispaced nodes, multiply the kernel matrix with the coefficients w
      */
@@ -209,12 +214,13 @@ void n_body_fft_2d(int N, int n_terms, double *xs, double *ys, double *chargesQi
     delete[] fft_input;
     delete[] fft_output;
     delete[] mpol_sort;
-
+    END_TIME("FFT");
+    START_TIME
     /*
      * Step 3: Compute the potentials \tilde{\phi}
      */
-    for (int i = 0; i < N; i++) {
-        int box_idx = point_box_idx[i];
+    PARALLEL_FOR(nthreads,N, {
+        int box_idx = point_box_idx[loop_i];
         int box_i = box_idx % n_boxes;
         int box_j = box_idx / n_boxes;
         for (int interp_i = 0; interp_i < n_interpolation_points; interp_i++) {
@@ -223,14 +229,15 @@ void n_body_fft_2d(int N, int n_terms, double *xs, double *ys, double *chargesQi
                     // Compute the index of the point in the interpolation grid of points
                     int idx = (box_i * n_interpolation_points + interp_i) * (n_boxes * n_interpolation_points) +
                               (box_j * n_interpolation_points) + interp_j;
-                    potentialQij[i * n_terms + d] +=
-                            x_interpolated_values[interp_i * N + i] *
-                            y_interpolated_values[interp_j * N + i] *
+                    potentialQij[loop_i * n_terms + d] +=
+                            x_interpolated_values[interp_i * N + loop_i] *
+                            y_interpolated_values[interp_j * N + loop_i] *
                             y_tilde_values[idx * n_terms + d];
                 }
             }
         }
-    }
+    });
+    END_TIME("Step 3");
     delete[] point_box_idx;
     delete[] x_interpolated_values;
     delete[] y_interpolated_values;

src/nbodyfft.h

@@ -20,7 +20,7 @@ void precompute_2d(double x_max, double x_min, double y_max, double y_min, int n
 
 void n_body_fft_2d(int N, int n_terms, double *xs, double *ys, double *chargesQij, int n_boxes,
                    int n_interpolation_points, double *box_lower_bounds, double *box_upper_bounds,
-                   double *y_tilde_spacings, complex<double> *fft_kernel_tilde, double *potentialQij);
+                   double *y_tilde_spacings, complex<double> *fft_kernel_tilde, double *potentialQij, unsigned int nthreads);
 
 void precompute(double y_min, double y_max, int n_boxes, int n_interpolation_points, kernel_type kernel,
                 double *box_lower_bounds, double *box_upper_bounds, double *y_tilde_spacing, double *y_tilde,

src/parallel_for.h

@@ -0,0 +1,43 @@
+#ifndef PARALLEL_FOR_H
+#define PARALLEL_FOR_H
+#include<algorithm>
+#include <thread>
+#include <vector>
+#if defined(_OPENMP)
+#pragma message "Using OpenMP threading."
+#define PARALLEL_FOR(nthreads,LOOP_END,O) {          			\
+	if (nthreads >1 ) {						\
+		_Pragma("omp parallel num_threads(nthreads)")		\
+		{							\
+		_Pragma("omp for")					\
+		for (int loop_i=0; loop_i<LOOP_END; loop_i++) {		\
+			O;						\
+		}							\
+		}							\
+	}else{								\
+		for (int loop_i=0; loop_i<LOOP_END; loop_i++) {		\
+			O;						\
+		}							\
+	}								\
+}
+
+#else
+#define PARALLEL_FOR(nthreads,LOOP_END,O) {          			\
+	if (nthreads >1 ) {						\
+        std::vector<std::thread> threads(nthreads);		        \
+		for (int t = 0; t < nthreads; t++) {			\
+		    threads[t] = std::thread(std::bind(			\
+		    [&](const int bi, const int ei, const int t) { 		\
+			    for(int loop_i = bi;loop_i<ei;loop_i++) {   O;  }	\
+		    },t*LOOP_END/nthreads,(t+1)==nthreads?LOOP_END:(t+1)*LOOP_END/nthreads,t)); \
+		}							\
+		std::for_each(threads.begin(),threads.end(),[](std::thread& x){x.join();});\
+	}else{								\
+		for (int loop_i=0; loop_i<LOOP_END; loop_i++) {		\
+			O;						\
+		}							\
+	}								\
+}
+
+#endif
+#endif

src/time_code.h

@@ -0,0 +1,20 @@
+#ifndef TIME_CODE_H
+#define TIME_CODE_H
+#include <chrono>
+#if defined(TIME_CODE)
+#pragma message "Timing code"                                   
+#define INITIALIZE_TIME std::chrono::steady_clock::time_point STARTVAR;                       
+#define START_TIME           			\
+                STARTVAR = std::chrono::steady_clock::now();
+
+#define END_TIME(LABEL) {          			\
+                std::chrono::steady_clock::time_point ENDVAR = std::chrono::steady_clock::now();      \
+                printf("%s: %lld ms\n",LABEL, std::chrono::duration_cast<std::chrono::milliseconds>(ENDVAR-STARTVAR).count());       \
+}                                                                       
+#else 
+#define INITIALIZE_TIME                      
+#define START_TIME           			
+#define END_TIME(LABEL) {}                                                                       
+
+#endif
+#endif