How to optimise the re-ordering of an array?

I want to optimize the re-ordering of some data arrays that hold about 4 million unsigned shorts. The aim is to process a data stream by bringing values that should be similar to each other to be close to each other. The pseudo-code is like this:

  for( i=0; i<n; i++)
    dest[i] = src[ idx[i] ] ;

To optimize the code for a specific list of idx[i] I tried to compile a 4 million line c function with the idx values filled in:

void reorder( unsigned short * restrict i, unsigned short * restrict o) {
  o[0]=i[2075723];
  o[1]=i[2075724];
  o[2]=i[2075722];
  ...
  o[4194301]=i[4192257];
  o[4194302]=i[4192256];
  o[4194303]=i[4190208];
 }

I had hoped to get GCC to create a clever stream of pshufw/pblend/unpack instructions ... instead it hangs after using up a lot of memory (7 GB). I was trying to make copy based version to avoid the complications of doing swaps in place.

Would anyone be able to suggest good ways to produce optimized the code for this problem? So far I tried:

• ordered reading, random writes : 60 ms (openmp didn't help)
• ordered writing, random reads : 20 ms (openmp -> 4 ms)

I was hoping to end up getting closer to the memory bandwidth (order 0.4 ms). A scheme that takes into account cache size and does blocking should help but I don't know where to start for designing one to do it. I also wonder if there is a simple way to exploit SIMD instructions?

Making a toy example with transpose I couldn't even get gcc to output an SIMD version, see:

https://godbolt.org/z/bzGWad

Is this a difficult problem for compilers or am I missing something simple?

Edit 21/11/2018 Added a complete but minimal example of the problem

Here is a complete example of the problem I am trying to optimise. In reality the ordering is a more complicated function, but the point is just to order data pixels according to their distance from the image center, like unwinding a spiral.

#include <omp.h>
#include <vector>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <algorithm>

#define N 2048

// Sorting on output, one core
void reorder_simple( const std::vector<size_t> &indices,
             const unsigned short input[],
             unsigned short output[]){
  for( int i=0; i<N*N; i++)
    output[i] = input[ indices[i] ];
}
// Sorting on output write, many cores
void reorder_omp( const std::vector<size_t> &indices,
          const unsigned short input[],
          unsigned short output[]){
#pragma omp parallel for
  for( int i=0; i<N*N; i++)
    output[i] = input[ indices[i] ];
}
// Benchmark for memory throughput, one core
void copy_simple(  const std::vector<size_t> &indices,
           const unsigned short input[],
           unsigned short output[]){
  for( int i=0; i<N*N; i++)
    output[i] = input[i];
}
// Benchmark for memory throughput, many cores
void copy_omp (  const std::vector<size_t> &indices,
         const unsigned short input[],
         unsigned short output[]){
#pragma omp parallel for
  for( int i=0; i<N*N; i++)
    output[i] = input[i];
}

// Macro to avoid retyping
#define bench(func)                                          \
  func( indices, input, output);                             \
  start = omp_get_wtime();                                   \
  for( size_t i=0; i<100; i++)                               \
      func( indices, input, output );                        \
  end =  omp_get_wtime();                                    \
  std:: cout << std::setw(15) << #func <<                    \
     ", Time taken: "  << (end-start)/100 << " /s\n";

int main()
{
  std::vector<float> sort_order(N*N);
  std::vector<size_t> indices(N*N);
  float radius, azimuth, ci, cj;
  double start, end;
  unsigned short *input, *output;

  ci = N*0.496;  // changes according to calibration
  cj = N*0.4985;  // reality is more complicated (tilts etc)
  for( size_t i=0; i<N; i++){
    for( size_t j=0; j<N; j++){
      radius  = sqrt( (i-ci)*(i-ci) + (j-cj)*(j-cj) );
      azimuth = atan2( i-ci, j-cj ); // from -pi to pi
      sort_order[i*N+j] = round( radius ) + azimuth/2/M_PI;
      indices[i*N+j] = i*N+j;
    }
  }
  // Find the order to sort data onto a spiral 
  std::sort( indices.begin(), indices.end(),
         [&sort_order](int i, int j){
           return sort_order[i] < sort_order[j]; });
  // Invent some test data
  input = new unsigned short [N*N];
  output = new unsigned short [N*N];
  for( size_t i=0 ; i<N*N; i++){
    input[i] = i;
    output[i]= 0;
  }
  // some timing:
  bench(reorder_simple);
  bench(reorder_omp)   ;
  bench(copy_simple)   ;
  bench(copy_omp)      ;
}


   % g++ reorder.cpp -o reorder -std=c++11 -O3 -march=native -fopenmp -Wall
   % ./reorder
     reorder_simple, Time taken: 0.0179023 /s
        reorder_omp, Time taken: 0.00349932 /s
        copy_simple, Time taken: 0.00140805 /s
           copy_omp, Time taken: 0.000250205 /s

I would like to get the reorder_omp function to be closer to the speed of the copy_omp function. Detectors can run at 500 frames per second so 3.5 ms is bad in comparison to 0.25 ms.

Edit again: 21/11/2018 the code to write the function that does not compile

  //top of file
  #include <fstream>  
  ...
  //just before the end: 
  std::ofstream out;
  out.open("cfunc.c");
  out << "void cfunc( unsigned short * restrict input,\n" <<
         "            unsigned short * restrict output){ \n"; 
  for(int i=0;i<N;i++)
    for(int j=0;j<N;j++)
      out << "output[" << i*N+j << "] = input[" << indices[i*N+j] << "];\n";
  out << "}\n";
  out.close();

Testing this on a different machine I am getting compiler errors from both gcc (7.3.0) and clang (6.0.0). It compiles and runs with tcc (0.9.27) but finishes slower than the looping over the indices.