А что тебе важнее, быстро обходить произвольный элемент или обойти все как можно быстрее? Я то говнокодом на питоне обходил, но у меня упора в производительность не было никакого. Так размерность понижай, если надо хитро ходить, но там ИМХО разная оптимизация будет для разных задач.

И то и то. Не, про питон речь вообще не идет.

Фиг знает, сходу ничего красивого в голову не приходит. Быстро загуглил, нашел это. Но там больше про GPU. Есть ещё википедия, как стартовая точка. https://en.wikipedia.org/wiki/Iterative_Stencil_Loops к сожалению конкретно эту задачу на скорость я не делал.

ЗЫ

Ты параллельно разные пиксели или воксели или что у тебя там обрабатывай

I - это некий гибрид Java-style итераторов и чего-то ещё:) Чаще всего такую сущность называют просто context, т.к. это он и есть, что-то что хранит сложное состояние отдельно от вычислений.

S - это просто адаптор в виде массива, в простонародии - сову на глобус.

Имхо, если задача стоит как сделать что-то с учётом соседей, эту сущность и надо инкапсулировать. Т.е. именно в элементе контейнера должна быть инкапсулированна эта взаимосвязь. Это можно сделать как храня ссылки на соседей так и лениво вычисляя в момент обращения.

Кстати, так и шаблон связи будет проще отслеживать, т.к. он будет локализован в одной точке.

Интерфейс видится как-то так:

#pragma omp parallel for firstprivate(I)
   for(size_t i=0; i<N; i++){ 
       auto& e = c[i];
       auto& n1 = e.neighbour1();
       ...
    }

Алыверды: с openmp никогда не работал за деньхи, как и с hpc кластерами, но на SO пишуть, что поддержка итераторов так-то есть в openmp 3.0, другое дело, что тебе не это нужно, а инкапсуляция связей элементов и итератор явно не то, что ищешь ты.

я, конечно, сварщик вообще не настоящий

но по-моему первая фраза «Есть многомерный массив, причем он может быть устроен достаточно сложным образом» - немного противоречива. Он массив, или он устроен сложным образом?

Те две дробилки, что я видел - были именно массивами. В одной из них для MPI-целей вся сетка делилась на блоки, блоки нумеровались Moreton curve. Связность и ближайшие соседи хранятся отдельными массивами (во второй вообще AMR не предполагалось, поэтому всё проще). И конечно никакие не плюсы.

По делу сказать нечего :)

Для обхода сложных контейнеров обычно используют итераторы, но это весьма криво параллелится средствами OpenMP

/thread

Иди, пили прямую параллельность для хер пойми чего средствами OpenMP.

Кроме того I может давать кучу дополнительной информации (координаты ячейки например), и тогда это уже больше чем итератор получается?

Итератор — это просто норматив интерфейса. Их, кстати, несколько разных: Input and output, Forward, Bidirectional, Random-access. Если можно ходить в обе стороны через ++, --, то это уже Bidirectional, а если ещё выполнять += offset, то Random-Access. По-моему Random-Access это как раз то что нужно для stenctil.

Дополнительно к интерфейсу итератора вы, конечно, можете предоставлять любой другой.

#pragma omp parallel for firstprivate(I)
   for(size_t i=0; i<N; i++){ 
       I.update(i); 
       ...
    }

Будучи настроенным, I ведет себя как умный указатель, то есть у него перегружены операции -> и *. М.б. его и инкрементировать можно - тогда это обычный итератор, но с возможностью update(i).

А что будет, если создать T=omp_get_max_threads() таких I, после чего создать массив Is = {I(0), I(N/T), ..., I(N - N/T)} и итерироваться как

#pragma omp parallel for 
for(int t = 0; t < T; ++t) {
    I_t = Is[t];
    for(size_t i=0; i<N/T; i++){ 
       I_t += i; 
       ...
    }
}

А что будет, если создать T=omp_get_max_threads() таких I, после чего создать массив Is = {I(0), I(N/T), …, I(N - N/T)} и итерироваться

это будет очень похоже на то что делает openMP по умолчанию (один из вариантов), но без балансировки. Ну и больше букв.

Быстро загуглил, нашел это. Но там больше про GPU.

Спасибо. КОллеги такие статьи пишут (лучше на самом деле), но я им все равно пошлю.

Ты параллельно разные пиксели или воксели или что у тебя там обрабатывай

Ды тред тому и посвящен с-но - как это делать с минимальным оверхедом и минимальным кодом о стороны юзера (меня);-)

Чаще всего такую сущность называют просто context, т.к. это он и есть, что-то что хранит сложное состояние отдельно от вычислений.

Ок. Правда у меня оно видимо так и останется iterator-ом, что бы можно было синтаксическими сахаром пользоваться.

Т.е. именно в элементе контейнера должна быть инкапсулированна эта взаимосвязь. Это можно сделать как храня ссылки на соседей так и лениво вычисляя в момент обращения.

В неструктурированных сетках (из тетраэдров) так и делают (хранят ссылки на соседей), от безысходности. Но в описываемом случае это сожрет очень много памяти, которой и так нет.

Ну и шаблоны иногда по 100 с гаком элементов содержат, от численной схемы зависит.

поддержка итераторов так-то есть в openmp 3.0

может быть, но ХЗ что они там спрятали под капот. Но посмотрю, спасибо.

Еще вопрос. Вот есть итератор/контекст I. Он кроме всего прочего содержит указатель на контейнер, указатель на ячейку и еще всяко разно - он быстрый, но толстый. И есть счетчик, он медленный (итератор на его основе с нуля долго порождать) но маленький (8 байт). Мне нужна промежуточная сущность, эдакий продвинутый счетчик - она не содержит пойнтеров, но на ее основе можно быстро получить доступ к элементу произвольного контейнера если у него правильный размер/структура. И при этом она не сильно больше счетчика. Ну и работать с ней (на уровне интерфейса) хочется как со счетчиком.

это будет очень похоже на то что делает openMP по умолчанию (один из вариантов), но без балансировки

Вы правы. Но с двойным циклом вроде можно[1] итерироваться по I вместо i. Массив Is кстати лучше выкинуть, тогда будет не особо больше букв. Что-то типа:

#pragma omp parallel for 
for(int t = 0; t < T; ++t) {
    for(iterator_t I_t = I(t*N/T); I_t < I((t+1)*N/T); ++I_t){ 
       //Работаем с итератором I_t
    }
}

Если уж вы пользуетесь итераторами, надо стремиться избавлятся от обхода массива по int i. Иначе вообще непонятно, нафига городить огород.

[1] Судя по SO, OpenMP 3.0 поддерживает parallel for для итераторов, но каждый тред будет проходить весь массив, что не эффективно. Двойной цикл позволит избавиться от этого оверхеда.

P.S. Окончательно избавиться от int i в интерфейсе и обходиться только I тоже возможно.

Потеря балансировки загрузки - это критический недостаток. Обеспечить такую балансировку на самом деле можно и без i - но под капотом все равно i будет сидеть в том или ином виде. Чудес не бывает.

Ну и отказываться от счетчиков из «религиозных» соображений… Я выше Параллельный итератор для stencil-вычислений? (комментарий) писал, что кроме всего прочего нужен некий «суперсчетчик». Правда OpenMP его тоже нормально не обработает.

Да, если смущает синтаксис - ничего не мешается завернуть это в в свой метод foreach принимающий юзеровскую функцию (напр лямбду), которая в свою очередь принимает итератор.

Обеспечить такую балансировку на самом деле можно и без i - но под капотом все равно i будет сидеть в том или ином виде.

Именно так и нужно сделать: в публичном интерфейсе I, а в имплементации (под капотом) уже i. Весь смысл итераторов в С++ именно в стандартизации интерфейса обхода контейнеров.

Еще вопрос. Вот есть итератор/контекст I. Он кроме всего прочего содержит указатель на контейнер, указатель на ячейку и еще всяко разно - он быстрый, но толстый. И есть счетчик, он медленный (итератор на его основе с нуля долго порождать) но маленький (8 байт). Мне нужна промежуточная сущность, эдакий продвинутый счетчик - она не содержит пойнтеров, но на ее основе можно быстро получить доступ к элементу произвольного контейнера если у него правильный размер/структура. И при этом она не сильно больше счетчика. Ну и работать с ней (на уровне интерфейса) хочется как со счетчиком.

Ваш промежуточный счётчик это функтор generic_I, возвращающий «нормальный» итератор, если его вызывать на контейнере I = generic_I(array). Интерфейс функтору вы можете приделать какой хотите. И да, ещё понадобится обратная функция generic_I = make_generic(I), «забывающая» конкретный массив.

Это как то сложно. В моем представлении суперсчетчик - это кусок итератора (некоторые поля итератора как структуры) + метод делания итератора и конструктор принимающий итератор. Да, конечно это можно и как функтор оформить - но можно и как аргумент в [] контейнера передавать.

В неструктурированных сетках (из тетраэдров) так и делают (хранят ссылки на соседей), от безысходности. Но в описываемом случае это сожрет очень много памяти, которой и так нет.

Дык тут как всегда, плати либо памятью либо ленивым вычислением по месту, вариантов особо то и нет, но на уровне интерфейса то вообще пофиг.

Еще вопрос. Вот есть итератор/контекст I. Он кроме всего прочего содержит указатель на контейнер, указатель на ячейку и еще всяко разно - он быстрый, но толстый. И есть счетчик, он медленный (итератор на его основе с нуля долго порождать) но маленький (8 байт). Мне нужна промежуточная сущность, эдакий продвинутый счетчик - она не содержит пойнтеров, но на ее основе можно быстро получить доступ к элементу произвольного контейнера если у него правильный размер/структура. И при этом она не сильно больше счетчика. Ну и работать с ней (на уровне интерфейса) хочется как со счетчиком.

Без конкретики - очень сложно. Есть ощущение, что ты просто пока не понял чего именно нужно, ну либо донести не можешь. Было бы очень круто сделать минимальный пример демонстрирующий проблемы которые решаются. В идеале в терминах розовых понни или чего-то ещё куда не придётся дополнительно вникать. Тогда будет гораздо проще манипулировать абстракциями.

Чисто на слух это звучит так, что ты хочешь выделить состояние из контекста. Т.е. контекст будет хранителем знания о том как из состояния почерпнуть связи, а его собственным состоянием будет только счётчик.

Плюсую

Вот такой минимальный пример:

template <typename T, int D> class Arr{
   std::vector<T> data;
   Ind<D> N, mul; // размеры и некий контекст, Ind<D> это int[D] обвешанный всякими плюшками
public:
   void init(Ind<D> N_){ 
      N = N_;  data.resize(N.prod());
      mul[0] = 1; for(int i=1; i<D; i++) mul[i] = mul[i-1]*N[i-1];
       
   }
   T& operator [](Ind<D> ijk){ return data[ijk*mul]; }

   struct iterator{
       size_t i; Ind<D> ijk;
       Arr *arr;
       T* operator -> (){ return &(arr.data[i]); }
       void operator ++(){ 
           i++; ijk[0]++; // это быстро
           if(ijk[0]==arr->N[0]) ijk = i%arr->N;  // юзаем перегруженный %, это долго
       }
       // доступ к соседям
       iterator operator + (Ind<D> delta){
           iterator ret = *this; ret.ijk += delta;
           ret.i = ret.ijk%arr->N;
           return ret;
       }
       bool is_bound_down(int axe){ return ijk[axe]==0; }
       bool is_out_down(int axe){ return ijk[axe]<0; }
       bool is_bound_up(int axe){ return ijk[axe]==arr->N[axe]-1; }
       bool is_out_up(int axe){ return ijk[axe]>=arr->N[axe]; }
   };
};

здесь состоянием является Ind ijk. Но можно к нему добавить счетчик i, где то удобнее через i, где то через ijk.

Не, это уже решение, надо поставить задачу и обозначить узкие места(узкие в плане дизайна). Другое дело, если достаточно хорошо тут сформулировать задачу, то и решение может оказаться очевидным.

А в приведённом чёрт ногу сломает, если не особо понимаешь, что в итоге должно получиться. Так-то всё ок, кроме того, что is_* функи надо унести в другую сущность, например какой-то «зерокост» мог бы быть в таком варианте:

template<typename D>
struct Iterator {
private:
Ind<D> ijk_;
public:
const Ind<D>& ijk() const {return ijk_;}
};

struct Bounds {
 Bounds(const Iterator& it) : it_(it) {}
 private:
 const Iterator& it_;
};

Или, можно просто bounds сделать другом iterator. Или даже вообще оставить как у тебя есть и не париться, пока не понятна проблема решаемая в общем.

P.S.: подсветка синтаксиса чуть чуть, но облегчает чтение.

Не, это уже решение, надо поставить задачу и обозначить узкие места(узкие в плане дизайна). Другое дело, если достаточно хорошо тут сформулировать задачу, то и решение может оказаться очевидным.

Дык если бы я бы мог это сделать то я бы и тред не создавал!

Давай я попробую набросать интерфейс, а ты покритикуешь?;-)

// это какой то D-мерный массив из ячеек типа T
// предполагается что таких массивов может быть много разных, 
// но они все имеют унифицированный интерфейс
template <typename T, int D> Arr{
...
public:
   size_t size() const;
   
   struct Index; // обеспечивает быстрый доступ к произвольному элементу
   T& operator [](Index); // вот его можно использовать
   Index get_index(Ind<D>) const;   // так их можно создавать
   Index get_index(size_t i) const; // или так 

   struct StencilPoint; // это точка шаблона (сдвиг относительно ячейки)
   StencilPoint get_stencil_point(Ind<D> delta);

   struct Cell: protected Index{ // дает доступ к ячейке и ее окружению
   protected:
      Arr *arr;
   public:
      T& operator *();
      T* operator ->();
      Ind<D> get_pos() const;

      // битовая маска указывающая на принадлежность к границе
      // по два бита на ось, первый слева второй справа
      uint32_t is_bound() const; 

      // медленный доступ к соседним ячейкам
      // это очень удобный классический интерфейс, но 
      // непонятно как сигнализировать о выходе за границу
      T& operator [](Ind<D> delta);

      // быстрый доступ к соседним ячейкам
      // но проблема с сигнализацией о выходе за границу та же
      T& operator [](StencilPoint delta);

      // проверка выхода за границу
      // с точки зрения производительности лучше
      // бы ее объединить с доступом к соседней ячейке
      bool check_out(Ind<D> delta);
      bool check_out(StencilPoint delta);

      // это неуклюжая попытка такого объединения
      // возвращается битовая маска в которой отмечены
      // промахи по осям (по два бита на ось)
      uint32_t get_nb(Ind<D> delta, T** ptr);
      uint32_t get_nb(StencilPoint delta, T** ptr);       
   };
   Cell get_cell(Ind<D>); // доступ к произвольной ячейке
   Cell get_cell(Index);  

   struct Iterator: public Cell{
      void operator ++ (); // для синтаксического сахара
      void operator = (size_t i); // для обхода обычными циклами
   };
   Iterator begin();
   Iterator end();
};

Во, кстати индекс, это сама по себе отличная абстракция. Посмотри на те же индексы к кутишным моделям, это как раз итератор и мета инфа. Притом и по многомерным массивам и по деревьям.

Т.е. Cell даже не обязательно делать адаптором с наследованием реализации, это может быть тупо индекс, т.к. это суть природа сущности хранить информацию для эффективного доступа. Просто многие забывают, что простой целочисленный индекс это не всегда самое эффективное представление такого состояния(в этот момент обычно где-то плачет разработчик баз данных).

Про оси мне уже сложно судить, т.к. с такими вещами я никогда не работал. Но по идее, я бы отказывался от адаптора в виде отдельного класса итератора, и работал бы всегда в терминах индексов. С моей колокольни Ind<D> это тоже индекс, просто другого типа.

Если продолжать аналогию, кажется, что если тебе нужна эмуляция плэйн доступа, то тебе сначала нужно построить простой массив индексов, а вот уже по нему итерироваться.

если тебе нужна эмуляция плэйн доступа, то тебе сначала нужно построить простой массив индексов, а вот уже по нему итерироваться.

Ну собственно так и делается, но неявным (ленивым) образом.

Насчет иерархии Index<-Cell<-Iterator.

Index (в моем тутошнем представлении, термин мб не вполне удачный) - это максимально компактная фиговина которая обеспечивает доступ к элементу. При этом Index не самодостаточен - в нем нет указателя на контейнер/элемент, это следствие компактности. Но подумав - да, ты прав - индексов может быть несколько разных, различающихся соотношением компактность/скорость доступа/удобство работы.

Cell - это прямой доступ к ячейке и ее соседям + выдача всякой информации о ячейке - расположение в пространстве, является ли граничной и тд. Но Cell не имеет интерфейса для перехода на следующую ячейку, это важно. Это попытка разделить действия над элементом и обход контейнера.

Iterator - это Cell которым можно обойти весь массив.

Собственно у меня такая идеология уже была для сложных контейнеров вроде AMR, но тут появилось три новых момента:

1. Индексов может быть несколько разных, зависит от контейнера и от задачи. Видимо им надо будет давать описательные имена.

2. Итератор может эффективно отслеживать обычный целочисленный индекс и такой индекс должен быть всегда, это дает параллельный обход.

3. Эту конструкцию надо единообразно прикрутить ко всем моим контейнерам.

Есть такая книженция - верёвка достаточной длинны, чтобы выстрелить себе в ногу. Там концепция граней очень хорошо имхо ложится на идею множественных индексов. Она конечно скорее художественная, но идеи у товарища интересные, как раз за вдохновением можно туда сгонять.

это максимально компактная фиговина которая обеспечивает доступ к элементу

Имхо это как раз и есть индекс. С тем лишь дополнением, что она позволяет доступиться элемента по какому-то аспекту, отличному от прямого обхода (например сортировка). Более того, индексы могут порождать другие индексы, другого типа. Т.е. индекс на элемент может содержать или вычислять индексы на соседей.

в нем нет указателя на контейнер/элемент, это следствие компактности

Так если ты например будешь строить тот же индекс сортировки, ты тоже возьмёшь второй массив, и заполнишь его например указателями на элементы но в этом массиве они будут лежать в порядке сортировки. Эти индексы (их роль тут будет выполнять просто указатель) тоже ничего не будут знать про контейнер. С указателем аналогия не очень, я её выбрал, чтобы не упоминать слово индекс два раза в одном предложении в разных значениях. Но если заменить указатели на индексы элементов в изначальном массиве, то получиться как раз такая же ситуация, индекс это просто чиселко, и если ты не знаешь как его интерпретировать то толку от него ноль.

Итератор может эффективно отслеживать обычный целочисленный индекс и такой индекс должен быть всегда, это дает параллельный обход.

В зависимости от задачи разумеется, но я бы обязательно рассматривал вариант построения именно вектора сложных индексов, видел такой подход вживую, уровня неделю подготовиться и за час долететь. Хотя хорошо оно работает когда вычисления надо повторять с минимальными изменениями контейнера конечно.

Ок, спасибо большое. Ушел дальше думать и книжку читать:-)