6
$\begingroup$

I have a simple but long function that takes a vector x[10], and outputs a vector y[100]. It is an automatically generated eval function for a multivariate polynomial, ie, there is only (complex) multiplications and additions involved (5000 of them), with many variables to cache intermediate results that are reused. The structure goes like this:

static inline void
evaluate(const complex x[10], complex y[100])
{
  const complex &X0 = x[0];   // shortcuts to input
  const complex &X1 = x[1];
  ...
  const complex &X9 = x[9];

  static constexpr complex C0 = 1;   // constants
  static constexpr complex C1 = -1;

  const complex G0 = C1 * X8;        // multiply. G standds for 'gate'
  const complex G1 = C0 + G0;        // add.
  const complex G2 = G1 * X9;
  ...
  const complex G5000 = G3427 + G3430 + G3433 + G3436 + G3438;

  y[0] = G104; // output
  y[1] = G135;
  y[2] = G195;
  ....
}

How to go about optimizing this function to the max on an Intel processor? For instance, reducing the number of multiplications or approximating them would help, since complex multiplications are expensive. What I've tried:

  • I use -O3 -ffast-math -march=native and other flags to GCC 5-7
  • I profiled this and it has mostly L1 data write miss, and some L1 data read miss.
  • Idea: the gates (G's) form a DAG, so if I recursively pack each set of independent variables into vectors, then one might be able to exploit SIMD instructions. Compilers do this during SLP auto-vectorization, see recent work (1)(2), but I am unsure how much of it they handle.

This is a visualization of the DAG for a small example:

  G2 = G1 * X9;
  G3 = X8 * X6;
  G5 = G1 * X0;
  G6 = X8 * X7;
  G8 = X4 * X4;
  G9 = X5 * X5;
  G11 = X6 * X6;
  G12 = G11 + X6;

Example DAG with 8 gates

Here is an example visualization of 100 of these gates: Example DAG with 100 gates Each rank/line (and more) can be vectorized at a time. Do modern systems optimize this?

Note: the variables are const meaning this straightline program (no loops) is in Static Single Assignment (SSA) form, a term used in compiler optimization.

After ordering the nodes for vectorization from the DAG (selecting a suitable topological ordering), I get 21 layers/ranks of vectorizable operations (1 rank not in table corresponds to the seeds: constants or vector X): layer Each entry in this table can be vectorized with SIMD, running from left to right. I'll be executing this function 1 billion times in my application.

Improve this question
$\endgroup$
11
  • 1
    $\begingroup$ I suggest that you give the full formula for the smallest example which contains all features of the general case. $\endgroup$
    – Carl Christian
    Commented Mar 7, 2019 at 18:17
  • 3
    $\begingroup$ I agree with @CarlChristian but will add the following: Use loops. They make the code a lot more readable, less error-prone and the compiler is 'conditioned' to make them efficient on modern superscalar architectures. Also: Vectorization requires independence of data, whereas recursion has data dependencies by definition. $\endgroup$
    – Nox
    Commented Mar 8, 2019 at 9:30
  • $\begingroup$ You can calculate by traversing the graph. It already eliminates duplication of results. Perhaps with 2 adjacency matrices. One with the edge operator (e.g. Plus or Times), one with the symbolic formula, Then step through from the outermost edges inward. $\endgroup$
    – Edmund
    Commented Mar 9, 2019 at 10:56
  • 1
    $\begingroup$ Some suggestions: try putting the expressions in terms of FMA instructions (fused multiply-add) if the compiler isn't doing that already. This looks like a problem that would benefit from FPGA (slower) or ASIC (more expensive), but, of course, it complicates the logistics of getting your program to be run by other users. If you want to follow the parallelization path (i.e, for each of the 21 independent layers) consider going to the GPU (communication overhead might be a downside for a small problem like this). $\endgroup$
    – e.tadeu
    Commented Mar 17, 2019 at 10:56
  • 1
    $\begingroup$ continuing on the CPU: if you are using double, consider changing that to float (if the precision loss is affordable), this would reduce memory usage, which would reduce cache misses. Another thing that would reduce memory usage would be to reorder the instructions, assign to y[i] as soon as an output is ready and reuse G variables space, "recycling" them. $\endgroup$
    – e.tadeu
    Commented Mar 17, 2019 at 11:08

3 Answers 3

Reset to default
5
+50
$\begingroup$

OK, you have a very nice problem, I tried to run some benchmarks. First, I don't have your parameters so I used your small example. Second, since you do not specify the language, I used C + GSL (since I'm not that familiar with C++)

#if __STDC_VERSION__ >= 199901L
#define _XOPEN_SOURCE 700
#else
#define _XOPEN_SOURCE 600
#endif /* __STDC_VERSION__ */
#include <time.h>
#include <gsl/gsl_vector.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>

gsl_vector_complex *A;
gsl_vector_complex *B;
gsl_vector_complex *C;

void evaluate(gsl_vector_complex *in, gsl_vector_complex *out)
{
    gsl_complex n_unity;
    gsl_complex unity;
    gsl_complex zero;

    GSL_SET_COMPLEX (&n_unity, -1, 0);
    GSL_SET_COMPLEX (&unity, 1, 0);
    GSL_SET_COMPLEX (&zero, 0, 0);

    const gsl_complex neg_x8 = gsl_complex_mul (n_unity, gsl_vector_complex_get (in, 8));

    gsl_vector_complex_set (A, 0, neg_x8);
    gsl_vector_complex_set (A, 1, neg_x8);
    gsl_vector_complex_set (A, 2, neg_x8);
    gsl_vector_complex_set (A, 3, gsl_vector_complex_get (in, 8));
    gsl_vector_complex_set (A, 4, gsl_vector_complex_get (out, 4));
    gsl_vector_complex_set (A, 5, neg_x8);
    gsl_vector_complex_set (A, 6, gsl_vector_complex_get (in, 8));
    gsl_vector_complex_set (A, 7, gsl_vector_complex_get (out, 7));
    gsl_vector_complex_set (A, 8, gsl_vector_complex_get (in, 4));
    gsl_vector_complex_set (A, 9, gsl_vector_complex_get (in, 5));
    gsl_vector_complex_set (A, 10, gsl_vector_complex_get (out, 10));
    gsl_vector_complex_set (A, 11, gsl_vector_complex_get (in, 6));
    gsl_vector_complex_set (A, 12, gsl_vector_complex_get (in, 6));

    gsl_vector_complex_set (B, 0, unity);
    gsl_vector_complex_set (B, 1, unity);
    gsl_vector_complex_set (B, 2, gsl_vector_complex_get (in, 9));
    gsl_vector_complex_set (B, 3, gsl_vector_complex_get (in, 6));
    gsl_vector_complex_set (B, 4, unity);
    gsl_vector_complex_set (B, 5, gsl_vector_complex_get (in, 0));
    gsl_vector_complex_set (B, 6, gsl_vector_complex_get (in, 7));
    gsl_vector_complex_set (B, 7, unity);
    gsl_vector_complex_set (B, 8, gsl_vector_complex_get (in, 4));
    gsl_vector_complex_set (B, 9, gsl_vector_complex_get (in, 5));
    gsl_vector_complex_set (B, 10, unity);
    gsl_vector_complex_set (B, 11, gsl_vector_complex_get (in, 6));
    gsl_vector_complex_set (B, 12, gsl_vector_complex_get (in, 6));

    gsl_vector_complex_set_basis (C, 1);
    gsl_vector_complex_set (C, 2, gsl_vector_complex_get (in, 9));
    gsl_vector_complex_set (C, 5, gsl_vector_complex_get (in, 0));
    gsl_vector_complex_set (C, 12, gsl_vector_complex_get (in, 6));

    gsl_vector_complex_mul (A, B);
    gsl_blas_zaxpy (unity, C, A);
    gsl_vector_complex_memcpy (out, A);

}

void evaluate_naive(gsl_vector_complex *in, gsl_vector_complex *out)
{
    gsl_complex n_unity;
    gsl_complex unity;
    gsl_complex zero;

    GSL_SET_COMPLEX (&n_unity, -1, 0);
    GSL_SET_COMPLEX (&unity, 1, 0);
    GSL_SET_COMPLEX (&zero, 0, 0);

    const gsl_complex neg_x8 = gsl_complex_mul (n_unity, gsl_vector_complex_get (in, 8));

    gsl_vector_complex_set (out, 0, neg_x8);
    gsl_vector_complex_set (out, 1, gsl_complex_add (unity, neg_x8));
    gsl_vector_complex_set (out, 2, gsl_complex_sub (gsl_vector_complex_get (in, 9), gsl_complex_mul (neg_x8, gsl_vector_complex_get (in, 9))));
    gsl_vector_complex_set (out, 3, gsl_complex_mul (gsl_vector_complex_get (in, 8), gsl_vector_complex_get (in, 6)));
    gsl_vector_complex_set (out, 5, gsl_complex_sub (gsl_vector_complex_get (in, 0), gsl_complex_mul (neg_x8, gsl_vector_complex_get (in, 0))));
    gsl_vector_complex_set (out, 6, gsl_complex_mul (gsl_vector_complex_get (in, 8), gsl_vector_complex_get (in, 7)));
    gsl_vector_complex_set (out, 8, gsl_complex_mul (gsl_vector_complex_get (in, 4), gsl_vector_complex_get (in, 4)));
    gsl_vector_complex_set (out, 9, gsl_complex_mul (gsl_vector_complex_get (in, 5), gsl_vector_complex_get (in, 5)));
    gsl_vector_complex_set (out, 11, gsl_complex_mul (gsl_vector_complex_get (in, 6), gsl_vector_complex_get (in, 6)));
    gsl_vector_complex_set (out, 12, gsl_complex_add (gsl_vector_complex_get (out, 11), gsl_vector_complex_get (in, 6)));
}

int main()
{
    struct timespec start, end;
    size_t time_spent;

    gsl_vector_complex *x;
    gsl_vector_complex *y;

    double x_,y_;
    gsl_complex val;

    const gsl_rng_type * T;
    gsl_rng * r;

    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    x = gsl_vector_complex_calloc (13);
    y = gsl_vector_complex_calloc (13);
    A = gsl_vector_complex_calloc (13);
    B = gsl_vector_complex_calloc (13);
    C = gsl_vector_complex_calloc (13);


    for (size_t iteration = 0; iteration < 100000000; iteration++)
    {
        for (size_t index = 0; index < 13; index++)
        {
            x_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            y_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            GSL_SET_COMPLEX(&val, x_, y_);

            gsl_vector_complex_set (x, index, val);
        }

        for (size_t index = 0; index < 13; index++)
        {
            x_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            y_ = gsl_ran_gaussian_ziggurat (r, 1.0);
            GSL_SET_COMPLEX(&val, x_, y_);

            gsl_vector_complex_set (y, index, val);
        }


        clock_gettime (CLOCK_MONOTONIC, &start);
#ifdef NAIVE
        evaluate_naive (x,y);
#else
        evaluate (x, y);
#endif
        clock_gettime (CLOCK_MONOTONIC, &end);
        time_spent = (end.tv_sec - start.tv_sec) * 1000000000 + (end.tv_nsec - start.tv_nsec);
        printf("%e\n", 1.0 * time_spent);

    }

    gsl_rng_free (r);
    gsl_vector_complex_free(x);
    gsl_vector_complex_free(y);
    gsl_vector_complex_free(A);
    gsl_vector_complex_free(B);
    gsl_vector_complex_free(C);
    return 0;
}

I tried both vectorized and naive implementations and here are my results:

gcc, no optimizations, naive version, mean = 270.6; std = 190.6 ns
gcc, no optimizations, vector version, mean = 318.1; std = 381.1 ns
gcc, optimized, naive version, mean = 106.9; std = 160.3 ns
gcc, optimized, vector version, mean = 154.5; std = 109.6 ns
gcc, optimized, machine-specific, naive version, mean = 105.2; std = 96.2 ns
gcc, optimized, machine-specific, vector version, mean = 154.8; std = 144.8 ns

icc, optimized, naive version, mean = 75.6; std = 99.4 ns
icc, optimized, vector version, mean = 164.6; std = 209.9 ns
icc, optimized, machine-specific, naive version, mean = 75.8; std = 189.5 ns
icc, optimized, machine-specific, vector version, mean = 166.9; std = 1850.2 ns // this std is somehow bad, however, I used my laptop to run the tests
icc, no optimizations, naive version, mean = 130.9; std = 153.6 ns
icc, no optimizations, vector version, mean = 266.6; std = 238.8 ns

Q.E.D. If you're not sure what to do, let the compiler (icc) do its dark magic. This case can be platform-specific, moreover, timings can change when you increase the number/size of arrays. P.S. There is one more possible solution to try: Intel provides its own functions for element-wise vector operations, that can improve the results, however, you have to rewrite the code. P.P.S. Actually, I completely forgot that icc by default doesn't follow IEEE standards (equivalent of -ffast-math) and can give you extremely incorrect results extremely quickly. Be aware.

Improve this answer
$\endgroup$
3
$\begingroup$

Once you're at the level of a long list of expressions, there is little you can still do other than hope that the compiler finds opportunities for optimization at the assembly level. This may bring you 10 or 20% but not something that's going to make a meaningful difference.

What you ought to do is go back to the original formula you're trying to implement. There may be more efficient ways to look at the whole problem, for example through the use of loops, iterations, reformulations of the mathematical statement, etc. This is the kind of thing that might be able to accelerate your function by a factor of ten -- factors that matter. But none of this is accessible any more in the code you show, and no compiler will be able to undo the transformations that were made to come up with the code you show.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I was wondering if compilers really optimize this automatically. I can't go back to the function as this is automatically generated. The generation is optimal in the sense that it is the cheapest evaluation of the nonlinear function (eg, like Horner's method). However, are there are other topological sorting of this chain of instructions that could yield better performance? I've seen compilers unable to optimize simpler situations. $\endgroup$
    – rfabbri
    Commented Mar 9, 2019 at 3:17
  • 5
    $\begingroup$ I wholeheartedly disagree with the statement that only orders of magnitude in speed-up matter. A speed up of 20%, or let's make it 30% for simplicity's sake, is a difference between starting a numerical experiment on Monday and having the results the following Monday, vs having them the following Thursday. $\endgroup$
    – Nox
    Commented Mar 9, 2019 at 8:59
  • $\begingroup$ @rfabbri: Yes, compilers do optimize such functions. They seek for common subexpressions, they vectorize operations, they interleave operations to keep the various parts of processors that can work in parallel busy, etc. Whether that finds anything of great use in the current context is a different matter. $\endgroup$
    – Wolfgang Bangerth
    Commented Mar 11, 2019 at 13:03
  • 2
    $\begingroup$ @Nox: I don't disagree with your statement. What I was trying to point out is that once you're at a level as low as the one shown in the post, there is really rather little you can do other than hope that your compiler finds opportunities. It is a rare program where you can find 30% of performance optimization without changing the overall algorithm, let alone one where that kind of performance can be found in just a single function. $\endgroup$
    – Wolfgang Bangerth
    Commented Mar 11, 2019 at 13:05
  • $\begingroup$ @WolfgangBangerth I am at the point where I need to extract all possible speed from this function. I already tuned the algorithm, profiled, etc. This function and similar ones are evaluated beyond a billion times in a typical application. I need to squeeze time at this point just as a game programmer needs to increase FPS once the game is ready for this. $\endgroup$
    – rfabbri
    Commented Mar 16, 2019 at 18:31
1
$\begingroup$

Would it be possible to do the 1 billion runs method-parallel? Even if there is no way to parallelize the interior of the method body, depending on your problem, you might reach some form of concurrency by evoking the method in parallel. If you do not re-use your output y[100] as an input x[10] to the function then nobody stops you from executing it within different threads/cores/machines etc.

pseudocode:

open_one_zillion_threads(); x = get_threadlocal_x(); y = get_threadlocal_y(); evaluate(const x[10], y[100]); store_y_in_proper_place(y[100]); rejoin_threads();

also:

What are your precision constraints concerning the output? Is the result necessary in double precision? Do you really need all aspects of the output? In many problems there is a reasonable tradeoff between precision and speed and even switching from double to single might boost performance.

another idea:

It should be possible to express the operations in a form using matrices and vectors. If it is all simple operations you might interpret your input array as a vector, and apply a series of matrix ops to it which represent the layers. These matrices will not be square matrices:

$$ y[100] = A( B( C( D ~x[10]) ) ) $$

This might be counter intuitive as there will be loads of zeros in your matrix, but it opens your problem up to use powerful and well optimized LA-Packages (Lapack / Blas etc.) which possibly bring you out-of-the-box parallelization on those matrix vector products and they might be better at using SIMD / AVX etc. If you somehow can make full use of AVX registers in the matrix-vector operations, then having a couple of them contain zeros does not hurt you too much. The Matrix approach would put you in a place where there is loads of online material and libraries and you are more flexible if you want to change certain parts of that method.

If you are on an intel cpu, then using their homemade intel compiler might give you another 50%.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.