I use C++ and CUDA/C and want to write code for a specific problem and I ran into a quite tricky reduction problem.

My experience in parallel programming isn't negligible but quite limited and I cannot totally forsee the specificity of this problem.
I doubt there is a convenient or even "easy" way to handle the problems I am facing but perhaps I am wrong.
If there are any resources (i.e. articles, books, web-links, ...) or key-words covering this or similar problems, please let me know.

I tried to generalize the whole case as good as possible and keep it abstract instead of posting too much code.

The Layout ...
--------------

I have a system of N inital elements and N result elements. (I'll use N=8 for example but N can be any integral value greater than three.)

```
static size_t const N = 8;
double init_values[N], result[N];

```

I need to calculate almost every (not all i'm afraid) unique permutation of the init-values without self-interference.

This means calculation `f(init_values[0],init_values[1])`, `f(init_values[0],init_values[2])`, ..., `f(init_values[0],init_values[N-1])`, `f(init_values[1],init_values[2])`, ..., `f(init_values[1],init_values[N-1])`, ... and so on.

This is in fact a virtual triangular matrix which has the shape seen in the following illustration.

```
 P     0    1    2    3    4    5    6    7
   |---------------------------------------
  0|   x 
   |
  1|   0    x
   |
  2|   1    2    x 
   |
  3|   3    4    5    x
   |
  4|   6    7    8    9    x
   |
  5|  10   11   12   13   14    x
   |
  6|  15   16   17   18   19   20    x
   |
  7|  21   22   23   24   25   26   27    x

```

Each element is a function of the respective column and row elements in `init_values`.

```
P[i] (= P[row(i)][col(i]) = f(init_values[col(i)], init_values[row(i)])

```

i.e.

```
P[11] (= P[5][1]) = f(init_values[1], init_values[5])

```

There are `(N*N-N)/2 = 28` possible, unique combinations (Note: `P[1][5]==P[5][1]`, so we only have a lower (or upper) triangular matrix) using the example `N = 8`.

The basic problem
-----------------

The result array is computed from P as a sum of the row elements minus the sum of the respective column elements.
For example the result at position 3 will be calculated as a sum of row 3 minus the sum of column three.

```
result[3] = (P[3]+P[4]+P[5]) - (P[9]+P[13]+P[18]+P[24])
result[3] = sum_elements_row(3) - sum_elements_column(3)

```

I tried to illustrate it in a picture with N = 4.

![Required triangluar reduction scheme](https://i.stack.imgur.com/HkZC4.png)

As a consequence the following is true:

* `N-1` operations (potential concurrent writes) will be performed on each `result[i]`
* `result[i]` will have `N-(i+1)` writes from subtractions and `i` additions
* Outgoing from each `P[i][j]` there will be a subtraction to `r[j]` and a addition to `r[i]`

This is where the main problems come into place: 

* Using one thread to compute each P and updating the result directly will result in multiple kernels trying to write to the same result location (N-1 threads each).
* Storing the whole matrix P for a subsequent reduction step on the other hand is very expensive in terms of memory consumption and therefore *impossible* for very large systems.

The idea of having a unqiue, shared result vector for each thread-block is impossible, too.
(N of 50k makes 2.5 billion P elements and therefore [assuming a maximum number of 1024 threads per block] a minimal number of 2.4 million blocks consuming over 900GiB of memory if each block has its own result array with 50k double elements.)

I think I could handle reduction for a more static behaviour but this problem is rather dynamic in terms of potential concurrent memory write-access.
(Or is it possible to handle it by some "basic" type of reduction?)

Adding some complications ...
-----------------------------

Unfortunatelly, depending on (arbitrary user) input, which is independant of the initial values, some elements of P need to be skipped.
Let's assume we need to skip permutations P[6], P[14] and P[18]. Therefore we have 24 combinations left, which need to be calculated.

How to tell the kernel which values need to be skipped?
I came up with three approaches, each having notable downsides if N is very large (like several ten thousands of elements).

### 1. Store all combinations ...

... with their respective row and column index `struct combo { size_t row,col; };`, that need to be calculated in a `vector<combo>` and operate on this vector. (used by the current implementation)

```
std::vector<combo> elements;
// somehow fill
size_t const M = elements.size();
for (size_t i=0; i<M; ++i)
{
    // do the necessary computations using elements[i].row and elements[i].col  
}

```

This solution consumes is consuming lots of memory since only "several" (may even be ten thousands of elements but that's not much in contrast to several billion in total) but it avoids 

* indexation computations
* finding of removed elements

for each element of P which is the downside of the second approach.

### 2. Operate on all elements of P and find removed elements

If I want to operate on each element of P and avoid nested loops (which i couldn't reproduce very well in cuda) I need to do something like this:

```
size_t M = (N*N-N)/2;
for (size_t i=0; i<M; ++i)
{
   // calculate row indices from `i`
   double tmp = sqrt(8.0*double(i+1))/2.0 + 0.5;
   double row_d = floor(tmp);
   size_t current_row = size_t(row_d);
   size_t current_col = size_t(floor(row_d*(ict-row_d)-0.5));
   // check whether the current combo of row and col is not to be removed
   if (!removes[current_row].exists(current_col))
   {
     // do the necessary computations using current_row and current_col
   }
}

```

The vector `removes` is very small in contrast to the `elements` vector in the first example but the additional computations to obtain `current_row`, `current_col` and the if-branch are very inefficient.
(Remember we're still talking about billions of evaluations.)

### 3. Operate on all elements of P and remove elements afterwards

Another idea I had was to calculate all valid and invalid combinations independently.
But unfortunately, due to summation errors the following statement is true:

```
calc_non_skipped() != calc_all() - calc_skipped()

```

**Is there a convenient, known, high performance way to get the desired results from the initial values?**

I know that this question is rather complicated and perhaps limited in relevance. Nevertheless, I hope some illuminative answers will help me to solve my problems.

---

The current implementation
--------------------------

Currently this is implemented as CPU Code with OpenMP. 
I first set up a vector of the above mentioned `combo`s storing every P that needs to be computed and pass it to a parallel for loop.
Each thread is provided with a private result vector and a critical section at the end of the parallel region is used for a proper summation.