Bad C++ code in benchmark

[oleid]

With great interest I read the following blog post:

https://lpython.org/blog/2023/07/lpython-novel-fast-retargetable-python-compiler/

But to my surprise I discovered the C++ code in the benchmark of dijkstra_shortest_path is pessimized rather than optimized.

The usage of maps makes this Code much, much slower as it could.
On my system (some ancient i7 laptop), the following is more than 40 times faster (see below) than the original code.

#include <iostream>
#include <vector>

int32_t dijkstra_shortest_path(int32_t n, int32_t source) {
    int32_t i, j, v, u, mindist, alt, dummy, uidx;
    std::vector<int32_t>  dist, prev;
    std::vector<bool> visited;
    std::vector<int32_t> Q, graph(n*n);
    Q.reserve(n);
    dist.reserve(n);
    prev.reserve(n);
    
    for(i = 0; i < n; i++) {
        for(j = 0; j < n; j++) {
            graph[n * i + j] = std::abs(i - j);
        }
    }

    for(v = 0; v < n; v++) {
        dist.push_back(2147483647);
        prev.push_back(-1);
        Q.push_back(v);
        visited.push_back(false);
    }
    dist[source] = 0;

    while(Q.size() > 0) {
        u = -1;
        mindist = 2147483647;
        for(i = 0; i < Q.size(); i++) {
            if( mindist > dist[Q[i]] ) {
                mindist = dist[Q[i]];
                u = Q[i];
                uidx = i;
            }
        }
        Q.erase(Q.begin() + uidx);
        visited[u] = true;

        for(v = 0; v < n; v++) {
            if( v != u and not visited[v] ) {
                alt = dist[u] + graph[n * u + v];

                if( alt < dist[v] ) {
                    dist[v] = alt;
                    prev[v] = u;
                }
            }
        }
    }

    int32_t dist_sum = 0;
    for(i = 0; i < n; i++) {
        dist_sum += dist[i];
    }
    return dist_sum;
}


int main() {
    int32_t n = 4000;
    int32_t dist_sum = dijkstra_shortest_path(n, 0);
    std::cout<<dist_sum<<std::endl;
    return 0;
}

For convenience, this is the diff:

--- a.cpp	2023-07-30 10:59:17.482487750 +0200
+++ f.cpp	2023-07-30 14:57:58.543699027 +0200
@@ -1,13 +1,15 @@
 #include <iostream>
-#include <unordered_map>
 #include <vector>
 
 int32_t dijkstra_shortest_path(int32_t n, int32_t source) {
     int32_t i, j, v, u, mindist, alt, dummy, uidx;
-    std::unordered_map<int32_t, int32_t> graph, dist, prev;
-    std::unordered_map<int32_t, bool> visited;
-    std::vector<int32_t> Q;
-
+    std::vector<int32_t>  dist, prev;
+    std::vector<bool> visited;
+    std::vector<int32_t> Q, graph(n*n);
+    Q.reserve(n);
+    dist.reserve(n);
+    prev.reserve(n);
+    
     for(i = 0; i < n; i++) {
         for(j = 0; j < n; j++) {
             graph[n * i + j] = std::abs(i - j);
@@ -15,10 +17,10 @@
     }
 
     for(v = 0; v < n; v++) {
-        dist[v] = 2147483647;
-        prev[v] = -1;
+        dist.push_back(2147483647);
+        prev.push_back(-1);
         Q.push_back(v);
-        visited[v] = false;
+        visited.push_back(false);
     }
     dist[source] = 0;

Furthermore, I implemented a rust version for comparison.

fn dijkstra_shortest_path(n: i32, source: i32) -> i32 {
    let mut dist = Vec::with_capacity(n as usize);
    let mut prev = Vec::with_capacity(n as usize);
    let mut visited = Vec::with_capacity(n as usize);
    let mut q = Vec::with_capacity((n * n) as usize);
    let mut graph = vec![0; (n * n) as usize];

    for i in 0..n {
        for j in 0..n {
            graph[(n * i + j) as usize] = (i - j).abs();
        }
    }

    for v in 0..n {
        dist.push(2147483647);
        prev.push(-1);
        q.push(v);
        visited.push(false);
    }
    dist[source as usize] = 0;

    let mut u;
    let mut uidx = 0;
    let mut mindist;
    let mut alt;
    while q.len() > 0 {
        u = -1;
        mindist = 2147483647;
        for i in 0..q.len() {
            let d = dist[q[i] as usize];
            if mindist > d {
                mindist = d;
                u = q[i];
                uidx = i;
            }
        }
        q.remove(uidx);
        visited[u as usize] = true;

        for v in 0..n {
            if v != u && !visited[v as usize] {
                alt = dist[u as usize] + graph[(n * u + v) as usize];

                if alt < dist[v as usize] {
                    dist[v as usize] = alt;
                    prev[v as usize] = u;
                }
            }
        }
    }

    dist.iter().sum()
}

fn main() {
    let n = 4000;
    let dist_sum = dijkstra_shortest_path(n, 0);
    println!("{}\n", dist_sum);
}

Benchmark

$ clang++ -std=c++20 -ffast-math -O3 -march=native old.cpp   -o old
$ clang++ -std=c++20 -ffast-math -O3 -march=native better.cpp   -o better
$ rustc main.rs -C opt-level=3  -o better-rust
$ hyperfine --warmup=20 --min-runs=10  ./old ./better ./better-rust 
Benchmark 1: ./old
  Time (mean ± σ):      1.690 s ±  0.112 s    [User: 1.380 s, System: 0.299 s]
  Range (min … max):    1.407 s …  1.790 s    10 runs
 
Benchmark 2: ./better
  Time (mean ± σ):      40.3 ms ±   2.7 ms    [User: 31.5 ms, System: 8.2 ms]
  Range (min … max):    34.4 ms …  45.8 ms    71 runs
 
Benchmark 3: ./better-rust
  Time (mean ± σ):      40.4 ms ±   1.7 ms    [User: 32.8 ms, System: 7.5 ms]
  Range (min … max):    37.7 ms …  44.2 ms    68 runs
 
Summary
  ./better ran
    1.00 ± 0.08 times faster than ./better-rust
   41.93 ± 3.97 times faster than ./old

(compilers used: rustc 1.71.0-nightly (2f2c438dc 2023-05-08) and clang version 15.0.7 on ArchLinux)

[certik]

@oleid awesome, beautiful! Thank you for taking the time to optimize the benchmark, we'll use it.

We'll benchmark against a list based version in LPython as well. I am curious what the results will be.

[czgdp1807]

@oleid Thank you for the optimisation. I would like to mention some points,

1. I chose to use dict because I wanted to benchmark hash-maps of LPython and C++.
2. LPython doesn't support pre-allocating memory for list the way C++ does (by calling reserve). I am not sure Python has the ability to pre-allocate memory for list. If Python hasn't, then LPython can't support that. By pre-allocation here I mean actually allocating memory in Python for a list. If we want to keep reserve as a no-operation thing (as @certik suggested) in Python then its possible to support it in LPython.

So, now I would post 3 codes below,

Without calling reserve

#include <iostream>
#include <vector>

int32_t dijkstra_shortest_path(int32_t n, int32_t source) {
    int32_t i, j, v, u, mindist, alt, dummy, uidx;
    std::vector<int32_t>  dist, prev;
    std::vector<bool> visited;
    std::vector<int32_t> Q, graph;

    for(i = 0; i < n; i++) {
        for(j = 0; j < n; j++) {
            graph.push_back(std::abs(i - j));
        }
    }

    for(v = 0; v < n; v++) {
        dist.push_back(2147483647);
        prev.push_back(-1);
        Q.push_back(v);
        visited.push_back(false);
    }
    dist[source] = 0;

    while(Q.size() > 0) {
        u = -1;
        mindist = 2147483647;
        for(i = 0; i < Q.size(); i++) {
            if( mindist > dist[Q[i]] ) {
                mindist = dist[Q[i]];
                u = Q[i];
                uidx = i;
            }
        }
        Q.erase(Q.begin() + uidx);
        visited[u] = true;

        for(v = 0; v < n; v++) {
            if( v != u and not visited[v] ) {
                alt = dist[u] + graph[n * u + v];

                if( alt < dist[v] ) {
                    dist[v] = alt;
                    prev[v] = u;
                }
            }
        }
    }

    int32_t dist_sum = 0;
    for(i = 0; i < n; i++) {
        dist_sum += dist[i];
    }
    return dist_sum;
}


int main() {
    int32_t n = 20000;
    int32_t dist_sum = dijkstra_shortest_path(n, 0);
    std::cout<<dist_sum<<std::endl;
    return 0;
}

(lp) 20:22:17:~/lpython_project/lpython % clang++ -ffast-math -funroll-loops -O3 /Users/czgdp1807/lpython_project/dijkstra_cpp.cpp
(lp) 20:22:25:~/lpython_project/lpython % time ./a.out
199990000
./a.out  0.70s user 0.47s system 92% cpu 1.274 total

With calling reserve

#include <iostream>
#include <vector>

int32_t dijkstra_shortest_path(int32_t n, int32_t source) {
    int32_t i, j, v, u, mindist, alt, dummy, uidx;
    std::vector<int32_t>  dist, prev;
    std::vector<bool> visited;
    std::vector<int32_t> Q, graph(n*n);
    Q.reserve(n);
    dist.reserve(n);
    prev.reserve(n);

    for(i = 0; i < n; i++) {
        for(j = 0; j < n; j++) {
            graph[n * i + j] = std::abs(i - j);
        }
    }

    for(v = 0; v < n; v++) {
        dist.push_back(2147483647);
        prev.push_back(-1);
        Q.push_back(v);
        visited.push_back(false);
    }
    dist[source] = 0;

    while(Q.size() > 0) {
        u = -1;
        mindist = 2147483647;
        for(i = 0; i < Q.size(); i++) {
            if( mindist > dist[Q[i]] ) {
                mindist = dist[Q[i]];
                u = Q[i];
                uidx = i;
            }
        }
        Q.erase(Q.begin() + uidx);
        visited[u] = true;

        for(v = 0; v < n; v++) {
            if( v != u and not visited[v] ) {
                alt = dist[u] + graph[n * u + v];

                if( alt < dist[v] ) {
                    dist[v] = alt;
                    prev[v] = u;
                }
            }
        }
    }

    int32_t dist_sum = 0;
    for(i = 0; i < n; i++) {
        dist_sum += dist[i];
    }
    return dist_sum;
}


int main() {
    int32_t n = 20000;
    int32_t dist_sum = dijkstra_shortest_path(n, 0);
    std::cout<<dist_sum<<std::endl;
    return 0;
}

(lp) 20:24:22:~/lpython_project/lpython % clang++ -ffast-math -funroll-loops -O3 /Users/czgdp1807/lpython_project/dijkstra_cpp.cpp
(lp) 20:24:26:~/lpython_project/lpython % time ./a.out
199990000
./a.out  0.48s user 0.13s system 99% cpu 0.623 total

Equivalent LPython code (without calling something like reserve)

from lpython import i32

def dijkstra_shortest_path(n: i32, source: i32) -> i32:
    i: i32; j: i32; v: i32; u: i32; mindist: i32; alt: i32; dummy: i32; uidx: i32
    dist_sum: i32;
    graph: list[i32] = []
    dist: list[i32] = []
    prev: list[i32] = []
    visited: list[bool] = []
    Q: list[i32] = []

    for i in range(n):
        for j in range(n):
            graph.append(abs(i - j))

    for v in range(n):
        dist.append(2147483647)
        prev.append(-1)
        Q.append(v)
        visited.append(False)
    dist[source] = 0

    while len(Q) > 0:
        u = -1
        mindist = 2147483647
        for i in range(len(Q)):
            if mindist > dist[Q[i]]:
                mindist = dist[Q[i]]
                u = Q[i]
                uidx = i
        dummy = Q.pop(uidx)
        visited[u] = True

        for v in range(n):
            if v != u and not visited[v]:
                alt = dist[u] + graph[n * u + v]

                if alt < dist[v]:
                    dist[v] = alt
                    prev[v] = u

    dist_sum = 0
    for i in range(n):
        dist_sum += dist[i]
    return dist_sum


def test():
    n: i32 = 20000
    print(dijkstra_shortest_path(n, 0))

test()

(lp) 20:26:35:~/lpython_project/lpython % time ./a.out
199990000
./a.out  0.67s user 0.15s system 99% cpu 0.822 total

Final Results

1. C++ code without calling reserve - 1.274 (relative - 1.55)
2. C++ code with calling reserve - 0.623 (relative - 0.76)
3. LPython code (equivalent to without calling reserve) - 0.822 (relative - 1.0)

[certik]

It looks like reserve makes a big difference here, so we should support it in LPython as well. For example by from lpython import list_reserve and just call this function. The CPython implementation will not do anything. I opened up an issue for this at #2232.

[oleid]

Instead of reserving and calling push_back in this case one could as well simply initialize directly, as pointed out here: https://www.phoronix.com/forums/forum/phoronix/latest-phoronix-articles/1400748-lpython-is-the-latest-python-implementation-aiming-to-be-very-fast-multiple-backends?p=1400898#post1400898

The following is c++20 due to usage of views

std::vector<int32_t> dist(n, std::numeric_limits<int32_t>::max());
std::vector<int32_t> prev(n, -1);
std::vector<int8_t> visited(n, false); // could omit false, as that is the default initialized value of bool, but does not hurt
auto integer_sequence = std::views::iota(0, n);
std::vector<int32_t> Q(integer_sequence.begin(), integer_sequence.end());

I didn't benchmark it, but should be faster.
Also, please note that vector<bool> is slow in C++, it's a trap.

[certik]

@oleid great. We'll benchmark the best C++ implementation that you (we) can come up with.

[NobodyXu]

@certik I also have a more optimized Rust-version with optimized Vec initialization, avoid array bound checking and use of conditional move to avoid branch mis-prediction:

Details

fn dijkstra_shortest_path(n: i32, source: i32) -> i32 {
    let mut dist = vec![i32::MAX; n as usize];
    let mut prev = vec![-1; n as usize];
    let mut visited = vec![false; n as usize];
    let mut q = Vec::with_capacity(n as usize);
    let mut graph = Vec::with_capacity((n as usize) * (n as usize));

    for i in 0..n {
        for j in 0..n {
            graph.push((i - j).abs());
        }
    }

    q.extend(0..n);
    dist[source as usize] = 0;

    while !q.is_empty() {
        let mut uidx = 0;
        let mut u = -1;
        let mut mindist = i32::MAX;

        q.iter().copied().enumerate().for_each(|(i, q_value)| {
            let d = dist[q_value as usize];

            let is_smaller = mindist > d;

            // Prompt compiler to use conditional move instead of branching,
            // since CPU cannot make correct prediction about this.
            mindist = mindist.min(d);
            u = if is_smaller { q_value } else { u };
            uidx = if is_smaller { i } else { uidx };
        });

        let u = u as usize;
        let n = n as usize;

        q.remove(uidx);
        visited[u] = true;

        // The loop below never modifies `dist[u as usize]` since `v != u`
        // and the loop only modifies `dist[v as usize]`, so it can be
        // considered a constant in the loop and load outside it.
        let dist_u = dist[u];

        for v in 0..n {
            if v != u && !visited[v] {
                let alt = dist_u + graph[n * u + v];

                let is_smaller = alt < dist[v];

                // Prompt compiler to use conditional move instead of branching,
                // since cpu cannot make correct prediction about this.
                dist[v] = if is_smaller { alt } else { dist[v] };
                prev[v] = if is_smaller { u as i32 } else { prev[v] };
            }
        }
    }

    dist.iter().sum()
}

fn main() {
    let n = 4000;
    let dist_sum = dijkstra_shortest_path(n, 0); // 7998000
    println!("{}\n", dist_sum);
}

[oleid]

@certik I also have a more optimized Rust-version with optimized Vec initialization, avoid array bound checking and use of conditional move to avoid branch mis-prediction:

Interesting, on my system that's actually about 50% slower on my system 🤔

For some reason pushing elements to the graph (although being pre-allocated) is slower than setting the elements to zero and overwriting 🤷

The whole while loop didn't change a thing for me. The compiler seems to be smart enough to eliminate all the bound checks by itself.

Anyway, the following is the fastest rust version on my system (and also in-line with the c++ version posted below):

fn dijkstra_shortest_path(n: i32, source: i32) -> i32 {
    let mut dist = vec![i32::MAX; n as usize];
    let mut prev = vec![-1; n as usize];
    let mut visited = vec![false; n as usize];
    let mut q = Vec::with_capacity(n as usize);
    let mut graph = vec![0; (n * n) as usize];

    for i in 0..n {
        for j in 0..n {
            graph[(n * i + j) as usize] = (i - j).abs();
        }
    }

    q.extend(0..n);
    dist[source as usize] = 0;

    let mut u;
    let mut uidx = 0;
    let mut mindist;
    let mut alt;
    while q.len() > 0 {
        u = -1;
        mindist = i32::MAX;
        for i in 0..q.len() {
            let d = dist[q[i] as usize];
            if mindist > d {
                mindist = d;
                u = q[i];
                uidx = i;
            }
        }
        q.remove(uidx);
        visited[u as usize] = true;

        for v in 0..n {
            if v != u && !visited[v as usize] {
                alt = dist[u as usize] + graph[(n * u + v) as usize];

                if alt < dist[v as usize] {
                    dist[v as usize] = alt;
                    prev[v as usize] = u;
                }
            }
        }
    }

    dist.iter().sum()
}

fn main() {
    let n = 4000;
    let dist_sum = dijkstra_shortest_path(n, 0);
    println!("{}\n", dist_sum);
}

[oleid]

@oleid great. We'll benchmark the best C++ implementation that you (we) can come up with.

The version with the view didn't work for me. The follwing did:

--- f.cpp	2023-07-30 14:57:58.543699027 +0200
+++ g.cpp	2023-07-31 12:36:02.138319588 +0200
@@ -1,27 +1,23 @@
 #include <iostream>
 #include <vector>
+#include <numeric>
+#include <limits>
 
 int32_t dijkstra_shortest_path(int32_t n, int32_t source) {
     int32_t i, j, v, u, mindist, alt, dummy, uidx;
-    std::vector<int32_t>  dist, prev;
-    std::vector<bool> visited;
-    std::vector<int32_t> Q, graph(n*n);
-    Q.reserve(n);
-    dist.reserve(n);
-    prev.reserve(n);
-    
+    std::vector<int32_t> dist(n, std::numeric_limits<int32_t>::max());
+    std::vector<int32_t> prev(n, -1);
+    std::vector<int8_t> visited(n, false); // could omit false, as that is the default initialized value of bool, but does not hurt
+    std::vector<int32_t> Q(n);
+    std::vector<int32_t> graph(n*n);
+    std::iota(Q.begin(), Q.end(), 0);
+        
     for(i = 0; i < n; i++) {
         for(j = 0; j < n; j++) {
             graph[n * i + j] = std::abs(i - j);
         }
     }
 
-    for(v = 0; v < n; v++) {
-        dist.push_back(2147483647);
-        prev.push_back(-1);
-        Q.push_back(v);
-        visited.push_back(false);
-    }
     dist[source] = 0;
 
     while(Q.size() > 0) {

My original version runs 15% slower than this one.

[NobodyXu]

Interesting, on my system that's actually about 50% slower on my system 🤔

Could it be that conditional move is actually slower than branching in this case?

For some reason pushing elements to the graph (although being pre-allocated) is slower than setting the elements to zero and overwriting 🤷

Perhaps there's just too many if-else in Vec::push for the compiler to optimize them

[certik]

@oleid, @NobodyXu, thank you both! We are implementing the "reserve" feature into LPython and we'll redo the benchmarks with your updated codes. I am curious what we will see!

[oleid]

I am curious what we will see!

Me too! And thanks again for your work on lpython! Very interesting stuff!

[maxbachmann]

LPython doesn't support pre-allocating memory for list the way C++ does (by calling reserve). I am not sure Python has the ability to pre-allocate memory for list. If Python hasn't, then LPython can't support that. By pre-allocation here I mean actually allocating memory in Python for a list. If we want to keep reserve as a no-operation thing (as @certik suggested) in Python then its possible to support it in LPython.

Python allows you to pre allocate the list at least from the C-API using PyList_New (https://docs.python.org/3/c-api/list.html#c.PyList_New). However it will be in an invalid state until you actually set all the elements. So either you will need to ensure the proper usage of this partially initialized list, or you simply initialize them to something like None:

PyObject*   list_reserve(Py_ssize_t len)
{
    PyObject* list = PyList_New(len);
    if(!list)
        return list;
    
    for(Py_ssize_t i = 0; i < len; ++i)
    {
        Py_INCREF(Py_None); // this is not needed in Python 3.12 since Py_None is now an immortal object
        PyList_SET_ITEM(list, i, Py_None);
    }
    return list;
}

Maybe from Python something like [None] * len generates somewhat similar bytecode 🤔
Cython optimizes from [None] * len to the code above as long as len is known to be an integer.

[certik]

@maxbachmann interesting, thanks! We could later see if we should internally transform [None] * len into what we are currently implementing in #2238, which will work like this:

from lpython import reserve
reserve(l, 100)

We don't support None in list elements, so it might create more confusion to use the syntax [None] * len, so the function reserve is probably less confusing.

[maxbachmann]

We don't support None in list elements, so it might create more confusion to use the syntax [None] * len, so the function reserve is probably less confusing.

The optimization is not specific to the None case. Basically whenever it is able to derive the output length it will use it for this kind of optimization. So the same would work for things like [1] * len. Another case where this can be interesting are comprehensions of known size like e.g. [i for i in range(n)]. The case of list comprehensions is not implemented in Cython yet either (cython/cython#2844).