3
\$\begingroup\$

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 <armadillo>
#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_num_threads(mkl_get_max_threads());
        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
LDFLAGS=-lfftw3 -lm -larmadillo -lpthread -lX11 -L/opt/boost/lib -lboost_system -L/opt/boost/lib -lboost_program_options -L/opt/intel/mkl/lib/intel64 -lmkl_rt
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
\$\endgroup\$
  • \$\begingroup\$ Yes, that is a typo.... And yes, I already looked at that, and linked with MKL \$\endgroup\$ – arc_lupus Nov 15 '17 at 10:17
  • 2
    \$\begingroup\$ Could you post a complete program? Even if not a full program, at least include the necessary headers! \$\endgroup\$ – Toby Speight Nov 15 '17 at 12:44
  • \$\begingroup\$ @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? \$\endgroup\$ – arc_lupus Nov 21 '17 at 19:57
  • \$\begingroup\$ 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. \$\endgroup\$ – Graipher Nov 23 '17 at 8:22
  • 2
    \$\begingroup\$ @Graipher: I added a MWE, which is a boiled-down version of my full code \$\endgroup\$ – arc_lupus Nov 23 '17 at 21:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.