Оптимизированный способ найти M самых больших элементов в NxN array с использованием C ++

Мне нужен невероятно быстрый способ найти 2D позиции и значения M самых больших элементов в массиве NxN.

прямо сейчас я делаю следующее:

struct SourcePoint {
    Point point;
    float value;
}

SourcePoint* maxValues = new SourcePoint[ M ];

maxCoefficients = new SourcePoint*[
for (int j = 0; j < rows; j++) {
    for (int i = 0; i < cols; i++) {
        float sample = arr[i][j];
        if (sample > maxValues[0].value) {
            int q = 1;
            while ( sample > maxValues[q].value && q < M ) {
                maxValues[q-1] = maxValues[q];      // shuffle the values back
                 q++;
            }
            maxValues[q-1].value = sample;
            maxValues[q-1].point = Point(i,j);
        }
    }
}

Структура Point - это всего два ints - x и y.

Этот код в основном выполняет сортировку вставкой входящих значений. MaxValues ​​[0] всегда содержит SourcePoint с наименьшим значением, которое по-прежнему удерживает его в пределах максимальных значений M, использованных до сих пор. Это дает нам быструю и легкую помощь, если sample <= maxValues, мы ничего не делаем. Проблема, с которой я сталкиваюсь, - это перетасовка каждый раз, когда обнаруживается новое лучшее значение. Он работает полностью вниз по maxValues, пока не находит нужное место, перетасовывая все элементы в maxValues, чтобы освободить место для себя.

Я подхожу к тому моменту, когда я готов изучить решения SIMD или оптимизацию кеша, поскольку похоже, что происходит изрядная загрузка кеша. Снижение стоимости этой операции резко повлияет на производительность моего алгоритма в целом, поскольку он вызывается много раз и составляет 60-80% моей общей стоимости.

Я пробовал использовать std :: vector и make_heap , но я думаю, что накладные расходы на создание кучи перевесили экономию операций с кучей. Вероятно, это потому, что M и N обычно не большие. M обычно составляет 10-20 и N 10-30 (NxN 100-900). Проблема в том, что эта операция вызывается повторно, и ее нельзя предварительно вычислить.

Я просто подумал предварительнозагрузить первые M элементов maxValues, что может дать небольшую экономию. В текущем алгоритме первые M элементов гарантированно перетасовываются вниз только для первоначального заполнения maxValues.

Любая помощь от гуру оптимизации будет принята с благодарностью :)