LINUX.ORG.RU

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

 , , ,


2

5

У меня есть программа на С++ для расчетов методом конечных элементов.
В силу того, что написана была не очень удачно (имеется плохие структуры данных и некоторые недостатки в алгоритме всей программы),
нужно ее переписать с использованием правильных алгоритмов и технологий.
В программе нужно перейти на использование blas для увеличения скорости работы с матрицами.
Возможно в будущем придется добавить возможность использования MPI или Cuda/OpenCL.
Также ее нужно сделать более универсальной.

Поэтому возникли следующие вопросы:

1. Какие технологии и библиотеки использовать для программы?
2. Какую реализацию работы с матрицами выбрать?
3. Какой язык программирования (c++, fortran, python) выбрать и какая библиотека blas лучше подойдет?
4. Использовать ли boost и итераторы?

★★★★★

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

Zodd ★★★★★ ()

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

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

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

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

Sphinx ★★☆☆ ()
Ответ на: комментарий от Sphinx

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

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от Sphinx

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

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

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

Zodd ★★★★★ ()

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

Не С++

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

Нет.

Имхо.

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

buddhist ★★★★★ ()

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

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

3. c++

4. нет

Reset ★★★★★ ()
Ответ на: комментарий от Zodd

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

Reset ★★★★★ ()
Ответ на: комментарий от Deleted

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

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

buddhist ★★★★★ ()
Ответ на: комментарий от Reset

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

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

ebantrop ()
Ответ на: комментарий от Reset

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

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

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

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

ebantrop ()
Ответ на: комментарий от ebantrop

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

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

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

на intel

Reset ★★★★★ ()
Ответ на: комментарий от Reset

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

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

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

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

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от Zodd

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

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

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

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

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

На кой?

Reset ★★★★★ ()
Ответ на: комментарий от Zodd

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

Reset ★★★★★ ()
Ответ на: комментарий от ebantrop

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

Reset ★★★★★ ()
Ответ на: комментарий от Reset

В 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);
      ...
	  }
Zodd ★★★★★ ()
Ответ на: комментарий от Zodd

Оно же под intel.

Есть момент.

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

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

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

ebantrop ()
Ответ на: комментарий от ebantrop

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

Zodd ★★★★★ ()
Ответ на: комментарий от Zodd

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

Reset ★★★★★ ()
Ответ на: комментарий от Zodd

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

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

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

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

ebantrop ()
Ответ на: комментарий от ebantrop

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от aedeph_

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

buddhist ★★★★★ ()
Ответ на: комментарий от Reset

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

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

Zodd ★★★★★ ()

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

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

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

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

VLev ()
Ответ на: комментарий от VLev

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от Zodd

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

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

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

VLev ()
Ответ на: комментарий от Zodd

Пример:

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

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);
  ...
}
и если iDef — степень двойки, то и работать будет куда быстрее. Вообще, такое впечатление, что это из fortran-а. В C пределы обычно меняются от 0 до iInt-1 (в моем «примере» уже учтено).

VLev ()
Ответ на: комментарий от aedeph_

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

buddhist ★★★★★ ()
Ответ на: комментарий от VLev

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от Reset

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

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от Zodd

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

ссылка на код

Reset ★★★★★ ()
Последнее исправление: Reset (всего исправлений: 2)
Ответ на: комментарий от Reset

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от Zodd

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

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

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

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

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

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

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

Reset ★★★★★ ()
Последнее исправление: Reset (всего исправлений: 1)
Ответ на: комментарий от Reset

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

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

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

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

Zodd ★★★★★ ()
Ответ на: комментарий от Zodd

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

Reset ★★★★★ ()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.