# Propagation solver in C++ using discrete Hankel transform

I wrote a small program which propagates light pulses using a discrete hankel transformation (based on A quasi-discrete Hankel transform for nonlinear beam propagation, You Kai-Ming et al., 2009. It can be boiled down to

• a matrix creation at the start of the program
• Multiple vector-vector-multiplications and matrix-vector-multiplications.

I would like to optimize the second part, after it is the most expensive one. It is defined as (in C++, using armadillo)

//Doing init stuff
arma::cx_colvec f, g, g2_vec, propagation_vector;
//Initialize vectors

//Propagation
for(size_t i = 0; i < rounds; ++i)
{
dht(f, g);
propagate(g, propagation_vector, g2_vec);
idth(g2_vec, f);
}


with the function dht() defined as

void dht(const arma::cx_colvec &in, arma::cx_colvec &out)//F and G are pre-initialized as arma::cx_colvec(size, arma::fill::zeros)
{
if(out.size() != in.size())
out = arma::cx_colvec(in.size());
F = in % (r_max / bessel_zeros);//% denoting the element-wise multiplication, r_max a constant double value, and bessel_zeros a constant double vector
G = c * F;//With c a constant double matrix and * a Matrix-vector-product
out = G % (bessel_zeros / rho_max);//% as element-wise multiplication, rho_max a constant double and bessel_zeros a constant double vector
}


the function propagate() defined as

void propagate(const arma::cx_colvec &in, const arma::cx_colvec &propagation_vector, arma::cx_colvec &out)
{
out = in % propagation_vector;
}


and the function idht() defined as

void idht(const arma::cx_colvec &in, arma::cx_colvec &out)
{
if(out.size() != in.size())
out = arma::cx_colvec(in.size());
F = in % (rho_max / bessel_zeros);
G = c * F;
out = G % (bessel_zeros / r_max);
}


According to valgrind, the functions dht() and idht() take ~40 % of the total running time of the program each, which I would like to reduce. Which optimization possibilities do I have here? I can run the program on a cluster, and use CUDA. Would any of that help?

At the moment the program takes ~110 seconds for vectors of the length 1024 and rounds = 2048;, with rapidly increasing times for longer vectors. Pre-calculating the dividing vectors in dht and idth shaves off ~1 second in total.

The program itself is compiled using

g++ -O2 -g -march=native -std=gnu++17 -fopenmp main.cpp -o main -lfftw3_omp -lfftw3 -lm -larmadillo -lgomp -lpthread -lX11 -L/opt/boost/lib -lboost_system -L/opt/intel/mkl/lib/intel64 -lmkl_rt


with g++ as g++-7.2.1

The code for a MWE is

#include <iostream>
#include <fftw3.h>
#include <vector>
#include <complex>
#include <cmath>
#include <omp.h>
#include <gsl/gsl_dht.h>
#include <chrono>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/program_options.hpp>
#include </opt/intel/mkl/include/mkl.h>
#include <immintrin.h>
#include <boost/math/interpolators/cubic_b_spline.hpp>
#include <boost/math/interpolators/barycentric_rational.hpp>
#define DIMENSION D2

#define D1 1
#define D2 2
#define INTERP4 0.2
#define INTERP10 0.3
#define INTERPOLATE_RATIO 3

class physics_parameters{
public:
double wavelength;
const double speed_of_light = 299792458;
const double tau_c = 3.5e-15;
double focal_length;
const double epsilon_0 = 8.85418782e-12;
const double pulse_range = 1e-3;
const double x_range = 3 * pulse_range;
const int propagation_rounds = 2048;
const double num_points = 500;
const double cut_out_range = 15000;
const double refractive_index = 3.52;
const double nonlinear_refractive_index = 4.5e-18;
const double electron_charge = 1.60217662e-19;
const double hbar_val = 1.0545718e-34;
const double two_photon_absorption_coefficient = 2e-11;
const double pulse_intensity = 1e7;
const double pulse_t0 = 0;
const double reduced_electron_mass = 0.15 * 9.10938356e-31;
const double initial_carrier_density = 1e16;
//Constants based on fixed data
double wave_number = 2 * M_PI / wavelength;
double omega_in_free_space = wave_number * speed_of_light;
double focal_length_outside_material = 2 * focal_length;
double dz_outside = focal_length_outside_material / propagation_rounds;

physics_parameters(const double wavelength, const double focal_length) : wavelength(wavelength), focal_length(focal_length)
{}

~physics_parameters()
{}

};

double gaussian_pulse(const double x_pos, const double intensity, const double t_0, const double pulse_width)
{
return intensity * exp(-(((x_pos - t_0) * (x_pos - t_0))/(pulse_width * pulse_width)));
}

std::complex<double> pulse_lense(const double &x_pos, const double k_0, const double focal_length, const double intensity, const double t_0, const double pulse_width)
{
return gaussian_pulse(x_pos, intensity, t_0, pulse_width) * std::exp(std::complex<double>(0, -1) * k_0 * pow(x_pos, 2)/(2 * focal_length));
}

class HT
{
private:

public:
arma::colvec p_N, bessel_r_forward, bessel_rho_forward, bessel_r_backward, bessel_rho_backward;
arma::mat c;
double J_N, S, rho_max, alpha_k;
const double num_points, r_max;
int k;
arma::cx_colvec f, g, F, G;
arma::colvec alpha_N;
arma::colvec bessel_zeros;

HT(const double num_points, const double r_max) : num_points(num_points), r_max(r_max)
{
alpha_N = arma::colvec(num_points + 1);
bessel_zeros = arma::colvec(num_points + 1);
bessel_r_backward = arma::colvec(num_points + 1);
bessel_r_forward = arma::colvec(num_points + 1);
bessel_rho_backward = arma::colvec(num_points + 1);
bessel_rho_forward = arma::colvec(num_points + 1);
p_N = arma::colvec(num_points + 1);
f = arma::cx_colvec(num_points + 1);
g = arma::cx_colvec(num_points + 1);
F = arma::cx_colvec(num_points + 1);
G = arma::cx_colvec(num_points + 1);
J_N = boost::math::cyl_bessel_j_zero(0., num_points + 1);
alpha_N[0] = 0;
for(size_t i = 1; i < num_points + 1; ++i)
{
alpha_N[i] = boost::math::cyl_bessel_j_zero(1., double(i));
}
for(size_t i = 0; i < num_points + 1; ++i)
bessel_zeros[i] = std::abs(boost::math::cyl_bessel_j(0, alpha_N[i]));
k = int(num_points/4);
alpha_k = alpha_N[k];
double S_tmp = 0;
for(size_t i = 1; i < num_points + 1; ++i)
S_tmp += pow(boost::math::cyl_bessel_j(0, alpha_k * alpha_N[i] / J_N), 2)/pow(boost::math::cyl_bessel_j(0, alpha_N[i]), 2);
S_tmp += 1;
S = 2/std::abs(boost::math::cyl_bessel_j(0, alpha_k)) * sqrt(S_tmp);

c = arma::mat(num_points + 1, num_points + 1);
for(size_t i = 0; i < num_points + 1; ++i)
for(size_t j = 0; j < num_points + 1; ++j)
c(i, j) = (2 * boost::math::cyl_bessel_j(0, alpha_N[i] * alpha_N[j]/S)/(S*std::abs(boost::math::cyl_bessel_j(0, alpha_N[i])) * std::abs(boost::math::cyl_bessel_j(0, alpha_N[j]))));

rho_max = S/(2 * M_PI * r_max);
bessel_r_forward = r_max / bessel_zeros;
bessel_r_backward = bessel_zeros / r_max;
bessel_rho_forward = bessel_zeros / rho_max;
bessel_rho_backward = rho_max / bessel_zeros;
//Initial values for hankel done
for(size_t i = 0; i < num_points + 1; ++i)
p_N(i) = (alpha_N[i]/(2*M_PI*rho_max));
}

~HT(){}

void create_data_points(arma::colvec &p_N) const
{
p_N = this->p_N;
}

void dht2(const arma::cx_colvec &in, arma::cx_colvec &out);
void idht2(const arma::cx_colvec &in, arma::cx_colvec &out);
void propagate(const arma::cx_colvec &in, const arma::cx_colvec &propagation_vector, arma::cx_colvec &out);
};

void HT::dht2(const arma::cx_colvec &in, arma::cx_colvec &out)
{
if(out.size() != in.size())
out = arma::cx_colvec(in.size());
F = in % (r_max / bessel_zeros);
G = c * F;
out = G % (bessel_zeros / rho_max);
}

void HT::idht2(const arma::cx_colvec &in, arma::cx_colvec &out)
{
if(out.size() != in.size())
out = arma::cx_colvec(in.size());
F = in % (rho_max / bessel_zeros);
G = c * F;
out = G % (bessel_zeros / r_max);
}

void HT::propagate(const arma::cx_colvec &in, const arma::cx_colvec &propagation_vector, arma::cx_colvec &out)
{
out = in % propagation_vector;
}

inline void propagate_pulse(arma::cx_colvec &f_vec, arma::cx_colvec &g_vec, arma::cx_colvec &g2_vec, HT &ht_matrix, const arma::cx_colvec &linear_propagation_vector)
{
ht_matrix.dht2(f_vec, g_vec);
ht_matrix.propagate(g_vec, linear_propagation_vector, g2_vec);
ht_matrix.idht2(g2_vec, f_vec);
}

void print_vector(const arma::colvec &x_axis, const arma::cx_colvec &vector_data, std::string filename)
{
std::ofstream out(filename.c_str());
for(size_t i = 0; i < vector_data.size(); ++i)
out << x_axis[i] << '\t' << vector_data[i].real() << '\t' << vector_data[i].imag() << '\t' << abs(vector_data[i]) << '\n';
out.close();
}

int main(int argc, char *argv[])
{
{
using namespace arma;

//Used variables
mkl_set_dynamic(false);
int num_points = 1024;

double wavelength = 2100e-9;
double focal_length = 8e-3;
double boost_num_points = 0;

physics_parameters parameters(wavelength, focal_length);
cx_colvec f_vec(num_points + 1), f_orig_vec(num_points + 1), g_vec(num_points + 1), g2_vec(num_points + 1), linear_propagation_vector_outside(num_points + 1);
colvec p_N(num_points + 1);
HT ht_matrix = HT(num_points, parameters.x_range);
ht_matrix.create_data_points(p_N);
for(size_t i = 0; i < p_N.size(); ++i)
f_vec(i) = pulse_lense(p_N[i], parameters.wave_number, parameters.focal_length, parameters.pulse_intensity, parameters.pulse_t0, parameters.pulse_range);
//For testing, if the program worked correctly
for(size_t i = 0; i < p_N.size(); ++i)
f_orig_vec(i) = f_vec(i);
for(size_t i = 0; i < p_N.size(); ++i)
{
linear_propagation_vector_outside(i) = std::exp(-std::complex<double>(0, 1) / (parameters.wave_number) * 2. * M_PI * M_PI * pow(parameters.dz_outside, 1) * pow(ht_matrix.alpha_N[i]/(2*M_PI*parameters.x_range), 2));
}

for(size_t i = 0; i < parameters.propagation_rounds; ++i)
{
propagate_pulse(f_vec, g_vec, g2_vec, ht_matrix, linear_propagation_vector_outside);
}

//For check of results
print_vector(p_N, f_vec, "f_vec.txt");
print_vector(p_N, f_orig_vec, "f_orig_vec.txt");

}
return 0;
}


compiled with the makefile

CC=gcc
CXX=g++
CFLAGS=-I$(DIR) -c -fopenmp CXXFLAGS=-I$(DIR) -O3 -g -march=native -std=gnu++17 -c
SOURCES=fftw_example.cpp
OBJECTS=fftw_example.o
EXECUTABLE=fftw_example

.PHONY: default all clean

default: all

all: $(SOURCES)$(EXECUTABLE)

$(EXECUTABLE):$(OBJECTS)
$(CXX)$(LDFLAGS) $(OBJECTS) -o$@

.cpp.o:
$(CXX)$(CXXFLAGS) $< -o$@

clean:
rm -f *o fftw_example

• Yes, that is a typo.... And yes, I already looked at that, and linked with MKL Commented Nov 15, 2017 at 10:17
• Could you post a complete program? Even if not a full program, at least include the necessary headers! Commented Nov 15, 2017 at 12:44
• @TobySpeight: I added the headers, a full program containing all the init routines for the vectors would be too long. After I am not interested in optimizing this initial part, I assume it can be omitted. Is that correct? Commented Nov 21, 2017 at 19:57
• Including a full, working, compilable program allows the reviewers to test their improvements. Otherwise they have no idea if their implementation is really better than armadillo's implementation of the % operator, for example. It is fine if the data is just randomly generated example data. Commented Nov 23, 2017 at 8:22
• @Graipher: I added a MWE, which is a boiled-down version of my full code Commented Nov 23, 2017 at 21:40

# Don't resize out unnecessarily

There is no need to check for out.size() in dht2() and idht2(); you are assigning a vector to out at the end anyway, which will take care of resizing it.

# Avoid divisions

CPUs nowadays can do many floating point multiplications and additions in a single clock cycle, but divisions are still very slow. Therefore, try to avoid dividing where possible. Sometimes you can precalculate the reciprocal of a number or vector and then multiply by that instead of doing the division. In fact, you already have those reciprocals ready, you just are not using them. Consider writing this instead:

void HT::dht2(const arma::cx_colvec &in, arma::cx_colvec &out)
{
F = in % bessel_r_forward;
G = c * F;
out = G % bessel_rho_forward;
}

void HT::idht2(const arma::cx_colvec &in, arma::cx_colvec &out)
{
F = in % bessel_rho_backward;
G = c * F;
out = G % bessel_r_backward;
}


# Use in-place operations where possible

You have vectors F and G that are used to hold intermediate results. However, you don't need to use them if you use in-place operations. If you cannot use in-place operations, then F and G might be useful as they avoid having to allocate memory, but that's about it. So in dht2(), I would write:

void HT::dht2(const arma::cx_colvec &in, arma::cx_colvec &out)
{
F = in % bessel_r_forward; // Avoid memory allocation for intermediate result
out = c * F;               // Store result in out
out %= bessel_rho_forward; // Multiply in-place
}


This avoids some memory usage and potentially also a bit of memory bandwidth.

# Precalculate as much as possible

Instead of multiplying the input and output vectors by bessel_*, consider pre-multiplying the matrix c such that you avoid the multiplications on the vectors.

# Avoid unnecessary calculations

If you inline all function calls inside propagate_pulse(), and make it a member of HT, you end up with:

void HT::propagate_pulse(arma::cx_colvec &f_vec, arma::cx_colvec &g_vec, arma::cx_colvec &g2_vec, const arma::cx_colvec &linear_propagation_vector) {
// dht2()
F = f_vec % bessel_r_forward;
G = c * F;
g_vec = G % bessel_rho_forward;

// propagate()
g_vec2 = g_vec % linear_propagation_vector;

// idht2()
F = g_vec2 % bessel_rho_backward;
G = c * F;
f_vec = G % bessel_r_backward;
}


There are some quite obvious optimizations to be done here. Consider the three element-wise multiplications in the middle; these commute, and bessel_rho_forward is the reciprocal of bessel_rho_backward, so these two cancel each other out:

void HT::propagate_pulse(arma::cx_colvec &f_vec, const arma::cx_colvec &linear_propagation_vector) {
F = f_vec % bessel_r_forward;
G = c * F;
F = G % linear_propagation_vector;
G = c * F;
f_vec = G % bessel_r_backward;
}


This already gets rid of two intermediate vectors, but you can avoid one more by doing the element-wise vector multiplications in-place:

void HT::propagate_pulse(arma::cx_colvec &f_vec, const arma::cx_colvec &linear_propagation_vector) {
f_vec %= bessel_r_forward;
F = c * f_vec;
F %= linear_propagation_vector;
f_vec = c * F;
f_vec %= bessel_r_backward;
}


And then there is the fact that you run propagate_pulse() in a loop. You could then make use of the fact that the multiplications between bessel_r_forward and bessel_r_backward cancel each other out, and that you multiply a vector by c twice in a row. So I would precalculate c2 = c * c, and then write:

void HT::propagate_pulse(arma::cx_colvec &f_vec, const arma::cx_colvec &linear_propagation_vector, std::size_t rounds) {
if (!rounds)
return;

// First round
f_vec %= bessel_r_forward;
F = c * f_vec;
F %= linear_propagation_vector;

// Subsequent rounds
for (std::size_t i = 1; i < rounds; ++i) {
F = c2 * F;
F %= linear_propagation_vector;
}

// Cleanup
f_vec = c * F;
f_vec %= bessel_r_backward;
}


Now consider that you can multiply the columns of c2 with the elements of linear_propagation_vector, and then raise the resulting matrix to the power of rounds - 1. This way you can get rid of the loop.