Переосмысление программы расчета

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

1. Какие технологии и библиотеки использовать для программы?

Главное это решать огромную разреженную линейную систему.
Для этого попробуй использовать CUSP, под GPU пока не много
вариантов по работе с разреженными матрицами.

Всё остальное там не так важно.

Триангуляцию тоже сам писал?

Для решения СЛАУ я использую lapack для ленточных матриц, хотя можно еще увеличить скорость, если использовать разреженные матрицы.

Триангуляцию тоже сам писал?

У меня пока реализованы для регулярных сеток, нужно будет написать триангуляцию либо приспособить существующие решения.

ленточных матриц
для регулярных сеток

А у тебя там именно МКЭ или всё же МКР?
Твёрдые тела или fluid?

У меня МКЭ для твердых тел.

для регулярных сеток

Я имел в виду, что в качестве КЭ брал прямоугольный элемент и использовал для прямоугольных подобластей, поэтому нумерация узлов шла как для регулярной сетки.

Какой язык программирования

Не С++

Использовать ли boost и итераторы?

Нет.

Имхо.

Я бы писал на Pure C или на Фортране.

Плюсую за Си. Фортран не помню уже.

1. различные варианты blas'ов, cuda sdk ...

2. никакую, абстрагируйся, чтобы любую реализацию blas можно было подключать

3. c++

4. нет

Почему матрица ленточная? В общем случае там черти что получаться будет. Надо работать с разреженными матрицами и решать чему-нибудь типа gmres.

Фортран годен, но, вроде как сейчас все основные фортранщики под виндой кодят во всяких Fortran PowerStation и Visual Fortran.

Под Unix лучше и стабильнее Си для целей ТС ничего не придумаешь :)

2. никакую, абстрагируйся, чтобы любую реализацию blas можно было подключать

По мне, так ничего лучше MKL для CPU пока не сделали, или это не так?

Смотря для чего. Если рассматривать только blas часть, то это не так. У меня atlas работает не хуже.

Если рассматривать только blas часть, то это не так.

Только blas не совсем интересно. В MKL есть полный LAPACK, векторные функции, случайные числа, FFT и еще много чего в IPP.

У меня atlas работает не хуже.

О как. На AMD или Intel'e?

В MKL есть полный LAPACK, векторные функции, случайные числа, FFT и еще много чего в IPP.

а оно всегда надо?

О как. На AMD или Intel'e?

на intel

никакую, абстрагируйся, чтобы любую реализацию blas можно было подключать

Тогда что задавать/передавать в качестве матрицы, double* ?

Вот еще в фортране доступ к элементам матрицы идет по столбцам, а в с++ построчно, поэтому приходится использовать транспонирование. Как решается эта проблема?

Итераторы тоже не нужно?

Почему матрица ленточная? В общем случае там черти что получаться будет. Надо работать с разреженными матрицами и решать чему-нибудь типа gmres.

Т.к. реализовать тогда проще было, но нужно будет переделать с разреженными матрицами.

Разве MKL не платный?

Тогда что задавать/передавать в качестве матрицы, double* ?

У себя можешь работать с чем хочешь, а при передаче в blas используй double*. Я, например, использую std::vector<>.

Как решается эта проблема?

В blas всё предусмотрено. Для этого там флажки передаются 'T' или 'N'

Итераторы тоже не нужно?

На кой?

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

а оно всегда надо?

А для чего может понадобиться голый blas?

Разве MKL не платный?

Платный, кроме линукса для некоммерческого пользования. А что делать?

Например, для МКЭ. LAPACK тут не нужен, так как матрицы разреженные. Алгоритмы по работе с разреженными матрицами обычно пишутся свои, так как библиотек таких не много, на низком уровне лучше использовать blas для производительности.

Оно же под intel. Я не хочу делать прогу железоориентированным.

В blas всё предусмотрено. Для этого там флажки передаются 'T' или 'N'

Я этого не заметил, спасибо. И последний вопрос у меня часто повторяется вложенный цикл в разных подпрограммах (3-4 вложения), а само тело вложения меняется (всунуть в один цикл не получается). Как нибудь это более красиво записать можно?

Пример:

    for (int m=1; m<=iInt; m++)
      for (int n=1; n<=iInt; n++)
	for (int l=1; l<=iInt; l++)
	  for (int k=1; k<=iDef; k++) {
	    locSig(k,1) = Sig(iDef*(iInt*iInt*iInt*(i-1)+iInt*iInt*(m-1)+iInt*(n-1)+l-1)+k,1);
      ...
	  }

Оно же под intel.

Есть момент.

Я не хочу делать прогу железоориентированным.

Лично у меня нет не сил ни ресурсов поддерживать кучу платформ для счетных программ. Мне надо либо забить на производительность в пользу железо-свободы либо сесть на определенное железо. Я выбрал второе в пользу intel+nvida. Скоро Xeon Phi говорят будет, посмотрим, может и nvidia ненужна будет.

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

на низком уровне лучше использовать blas для производительности

И ATLAS быстрее MKL? Может это оверхед OpenMP?

Да и вообще Cuda тоже не панацея. У меня есть несколько железок (например, ноут на amd+ati, а домашний декстоп intel+nvidia). Иногда приходится использовать ноут, чтобы прогить саму программу. Поэтому в данный момент я не хочу лишать возможности поработать вне дома.

Быстрее. Сравнивал на функции dgemm. Естественно, и там и там использовал многопоточный вариант.

Красиво не знаю. Можно писать коротко. Например, обернуть все циклы в шаблонную функцию matrix_generator и передавать в эту функцию код с помощью boost::bind или c++11 lambda.

Сравнивал на функции dgemm.

Каких размеров матрицы, если не секрет? Как ATLAS собран был?

Все спасибо. Пойду переписывать.

Да и вообще Cuda тоже не панацея.

CUDA изначально вендор-лок, в отличие от.

лишать возможности поработать вне дома

Ну оно работать будет и на AMD, просто не самым быстрым образом.

Каких размеров матрицы, если не секрет

Любых до 10000x10000.

Как ATLAS собран был?

Под проц q9550

Ну оно работать будет и на AMD, просто не самым быстрым образом.

Ясно. Мне кажется, что ради нескольких процент разницы в скорости не стоит затачивать программу под конкретное железо.

Под Unix лучше и стабильнее Си для целей ТС ничего не придумаешь :)

ifort? Не, не слышал.

Сколько фортранщиков им пользуется? А сколько фортранщиков пользуется gfortran'ом? У вас есть такие сведения? Мне действительно интересно, я в этом вопросе не слишком сведущ.

По поводу стиля:

1. Что используешь ++i или i++?
2. using namespace std; vector<> ... или std::vector<> ...

правильных алгоритмов и технологий

Каков целевой размер сетки (числа КЭ)?

какая библиотека blas лучше подойдет?

ЕМНИП, blas для плотных матриц, а МКЭ создает разреженные...

Число КЭ для трехмерного тела будет много, минимум 1000 элементов.

Чтобы получить разреженную матрицу, нужно перемножить еще много плотных матриц друг на друга, только при рассылке получается разреженная.

эээ, 1000 — это разве много? В МКР «много» начинается с 1000000000...
А сколько базисных функций в конечном элементе?

нужно перемножить еще много плотных матриц друг на друга

каков размер этих матриц?

Пример:

конкретно этот «пример» сворачивается до чего-то вроде:

const int Z=iInt*iInt*iInt*iDef;
for(int z=0, iz=i*Z; z<Z; z++) {
  locSig(z%iDef,0) = Sig(iz+z,0);
  ...
}

Почти все.

Почти все что? Пользуются ifort? Хм, это интересно, спасибо.

эээ, 1000 — это разве много? В МКР «много» начинается с 1000000000...

Не много, но меньше брать вообще нет смысла. А лучше всего 10^6. В МКР же это число узлов. А у одного элемента 8,20,... в зависимости от элемента.

Какую реализацию работы с матрицами выбрать?

никакую, абстрагируйся, чтобы любую реализацию blas можно было подключать

Подскажи как все это делается. Вроде умом понимаю, а в реализации затруднения.

Очень просто, я написал обертки над разными blas (включая cublas), которые на вход принимают double*. Выбор обертки происходит в compile-time с помощью страшных скриптов на cmake. Есть также небольшие шамаства с cuda, так как double*, который находится в видюхе нельзя использовать в обычном user-space.

ссылка на код

Спасибо. Буду разбираться.

Спасибо за хороший пример. Сразу видно, что проделана кропотливая работа. Возникло несколько вопросов:

1) Как ты применяешь свою библиотеку?
2) Ты реализуешь матрицу как обертку над *double через класс?
3) В каком месте у тебя используется std::vector?
4) Что такое Alloc?
5) Какая лицензия на код?
6) Что было вдохновителем для создания кода, boost?

1. как-то так. Это решение ур-я баратропного вихря на сфере (Навье-Стокс в терминах функции тока в сферических координатах). Где-то у меня в примерах для уравнения Лапласа был код, который может работать везде, в том числе и на GPU (nvidia).

2. без оберток

3. вместо выделения/освобождения памяти использую нечто похожее на std::vector, который сам выделяет и освобождает (exception-safe между прочим)

4. это хрень, которая может выделять память в оперативке или в видео-памяти, если используется cuda

5. sleepycat license, где не указано другого (в заголовках каждого исходника есть)

6. Я Александреску с Саттером обчитался и в те времена еще не отошел. Сейчас бы я такую жесть не написал.

Вообще, сейчас бы я наверно ничего писать не стал, потому что все изобретенные мной велосипеды вдруг появились в cuda-sdk, в том числе работа с разреженными матрицами.

2. без оберток

Тогда что это?

template < typename T, typename Alloc >
class Array
{
	typedef Array < T, Alloc > my_type;

	size_t size_;
	T * data_;
	Alloc alloc_;
...

Это нечто похожее на std::vector. Обертка над T * (обычно в качестве T выступает double). Служит для удобного автоматического выделения памяти на нужном устройстве, а также для автоматического освобождения при выходе из области видимости. Как с матрицей я с этим не работаю, у меня все алгоритмы на вход получают сырой указатель.