Как программисты работают с большими массивами? , Решение СЛАУ, не помещающихся в оперативной памяти

[Pacific]

experimenter
Ну еще можно сделать кластер из 16 компов, на каждом из которых будет по 2 Гб оперативки, объединить их в локальную сеть, и этот кластер тоже довольно быстро решит систему. Объем передаваемой между компами информации будет довольно небольшим. Так что автору прямая дорога на www.parallel.ru

[experimenter]

Цитата Pacific @ 06.11.07, 08:30

experimenter
Ну еще можно сделать кластер из 16 компов, на каждом из которых будет по 2 Гб оперативки, объединить их в локальную сеть, и этот кластер тоже довольно быстро решит систему. Объем передаваемой между компами информации будет довольно небольшим. Так что автору прямая дорога на www.parallel.ru

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

[Pacific]

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

Сообщение отредактировано: Pacific - 06.11.07, 08:47

[experimenter]

Мне самому стало интересно и я вбил в гугл следущие слова

huge system of linear equations

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

Добавлено 06.11.07, 08:49

Цитата Pacific @ 06.11.07, 08:46

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

Да, пожалуй, это самый пока перспективный вариант. другое дело, что в универах компы обычно не очень-то мощные... по 2 ГБ оперативки там точно не стоит.

[Pacific]

Здесь главное суммарный объем памяти. Но все равно, компов понадобится много, не во всяком универе столько будет

[experimenter]

Цитата Pacific @ 06.11.07, 08:52

Здесь главное суммарный объем памяти. Но все равно, компов понадобится много, не во всяком универе столько будет

Да, согласен. В общем, Telc, надо руководителю объяснять, что нужно для решения такой офигенной системы уравнений, а он уж пусть разруливает ситуацию дальше.
Pacific, а вы не могли бы прямую ссылку дать, где объясняется, как создать такую сетку?

[Pacific]

experimenter
А что здесь такого, обычная локальная сеть. Алгоритмы распределенных вычислений можно найти на том же parallel.ru или в гугле.

[experimenter]

Цитата Pacific @ 06.11.07, 09:11

experimenter
А что здесь такого, обычная локальная сеть. Алгоритмы распределенных вычислений можно найти на том же parallel.ru или в гугле.

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

[Pacific]

experimenter
Прямой ссылки я не знаю, но представляю себе распределенный метод Гаусса так:

1) Основная операция в методе Гаусса - сложение двух векторов C=A*t+B, где t - некоторый коэффициент.
2) Главный компьютер рассылает (UDP-броадкастом) вектор A всем остальным компьютерам, плюс он сообщает им, какой элемент сейчас ведущий. Остальные компьютеры определяют коэффициенты t для всех строк B своей части матрицы и выполняют операцию (1)
3) Когда шаг 2 закончен (остальные компы прислали сообщения главному), главный компьютер рассылает броадкаст с требованием предоставить ему (i+1)-ую строчку матрицы. Соответствующий компьютер отсылает данные.
4) Теперь можно повторить шаги 2-3 со следующей строчкой.

Объем передаваемой информации по сети - 2 строки матрицы на каждой итерации, что совсем немного, и обычная локалка это легко потянет.

Сообщение отредактировано: Pacific - 06.11.07, 09:26

[experimenter]

Цитата Pacific @ 06.11.07, 09:25

experimenter
Прямой ссылки я не знаю, но представляю себе распределенный метод Гаусса так:

1) Основная операция в методе Гаусса - сложение двух векторов C=A*t+B, где t - некоторый коэффициент.
2) Главный компьютер рассылает (UDP-броадкастом) вектор A всем остальным компьютерам, плюс он сообщает им, какой элемент сейчас ведущий. Остальные компьютеры определяют коэффициенты t для всех строк B своей части матрицы и выполняют операцию (1)
3) Когда шаг 2 закончен (остальные компы прислали сообщения главному), главный компьютер рассылает броадкаст с требованием предоставить ему (i+1)-ую строчку матрицы. Соответствующий компьютер отсылает данные.
4) Теперь можно повторить шаги 2-3 со следующей строчкой.

Объем передаваемой информации по сети - 2 строки матрицы на каждой итерации, что совсем немного, и обычная локалка это легко потянет.

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

[Telc]

Все написанное меня удивило.
Дело в том, что я на одной научно-практической конференции, к сожалению не в моем вузе и не в моем городе она проходила, слушал доклад, посвященный прикладной задаче из электромагнетизма. Людям приходится работать с системами нелинейных алгебраических уравнений. Решают они их методом Ньютона. На каждом шаге этого метода приходтся решать систему линейных уравнений. Вот они утверждают, что на обычном компе (с двухъядерным процессором) решают системы порядка 40000 уравнений и более за двое суток.
Используют модификацию метода Гаусса. Весь массив хронят в файле на винте и загружают столько строк, сколько помещается в оперативку, затем преобразовывают матрицу (приводят к трапецивидной форме) и записывают в файл. Так они в итоге получают треугольную матрицу.
(тезисы конференции к сожалению еще не вышли; да и ничего конкретного там расписано не будет, поскольку основные результаты их работы сводятся к физике).

[Pacific]

Telc
Это как раз то, что я описывал, только вместо кластера один компьютер, обрабатывающий матрицу кусками. Значит, кластер справится еще быстрее. Кстати, действительно, все строки матрицы, которые есть в оперативке, можно привести к такому виду, что останутся элементы только на главной диагонали, при небольшом числе обращений к диску (за один проход по файлу можно привести к такому виду весь кусок, находящийся в памяти).

[shadeofgray]

Попробуйте покопать в сторону блочно-матричных алгоритмов (на таких алгоритмах основан LAPACK 3). Это алгоритмы линейной алгебры, которые оперируют не со строками/столбцами, а с подматрицами матрицы, помещающимися в кэш. За счет этого достигается 2-3 кратное ускорение работы на больших матрицах, которые помещаются в оперативке, но не влезают в кэш процессора. В вашем случае ускорение будет ещё выше.

Я не могу сейчас в деталях описать эти алгоритмы, попробуйте поискать в lapack working notes, там есть pdf-Файлы для скачивания. Если кратко, то смысл в том, что большая матрица разбивается на две части - небольшая угловая подматрица, и огромный остаток. Сначала мы проводим операции над угловой подматрицей (например, строим её LU-разложение), затем обновляем остаток. После чего берем угловую подматрицу у остатка и обновляем уже её, и так далее. Смысл в том, что угловая подматрица влезает в кэш (в вашем случае - в оперативку) и считается очень быстро, а обновление остатка - по сути умножение огромной матрицы на матрицу небольшого размера, что можно очень эффективно проделать за один проход по остатку. Скажем, если у вас матрица NxN, а размер блока 10, то вам потребуется N/10 проходов по ней, чтобы всё рассчитать - N/10 считываний в оперативку с жесткого диска вместо N. Операции те же, просто переупорядочены так, чтобы работать более эффективно.

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

Я не уверен, что реализация, приведенная в LAPACK, подойдет вам без изменений (во-первых, она предполагает, что вся матрица помещается в оперативку; во-вторых, возможно, вы захотите оптимизировать структуру хранения матрицы на диске с учетом того, в каком порядке будут обходиться её элементы - чтобы по возможности считывать элементы в последовательных секторах). Кроме того, надо очень хорошо подумать над вопросом обусловленности вашей задачи - блочно-матричные алгоритмы менее стабильны, чем обычные, т.к. там имеют место определенные сложности с выбором ведущего элемента.

В общем, мне кажется, надо копать в этом направлении. Попробуйте разобраться, если будут вопросы - спрашивайте. Я с этой темой неплохо знаком, просто нет времени расписать всё в деталях Но на конкретные вопросы отвечу.

Добавлено 06.11.07, 16:02

Цитата Telc @ 06.11.07, 10:32

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

Да, это почти оно. По духу то же самое, немного другая реализация.

Добавлено 06.11.07, 16:06
По поводу двух стульев. Когда я говорил про операции с матрицами, помещающимися в кэш, то имел в виду, что обновление остатка также можно оптимизировать и реализовать в блочно-матричном стиле. Но, подумав, понял, что на самом деле здесь не тот случай, и определяющим фактором будет скорость работы с диском. Так что стул только один.

[AndreyK]

Цитата Pavia @ 05.11.07, 20:43

AndreyK
И где это вы такую глупость вычитали?

Здесь:
The VirtualAlloc function reserves or commits a region of pages in the virtual address space of the calling process. Memory allocated by this function is automatically initialized to zero, unless MEM_RESET is specified.

To allocate memory in the address space of another process, use the VirtualAllocEx function.

LPVOID VirtualAlloc(
LPVOID lpAddress, // region to reserve or commit
SIZE_T dwSize, // size of region
DWORD flAllocationType, // type of allocation
DWORD flProtect // type of access protection
);
Parameters
lpAddress
[in] Starting address of the region to allocate. If the memory is being reserved, the specified address is rounded down to the next 64-kilobyte boundary. If the memory is already reserved and is being committed, the address is rounded down to the next page boundary. To determine the size of a page on the host computer, use the GetSystemInfo function. If this parameter is NULL, the system determines where to allocate the region.
dwSize
[in] Size of the region, in bytes. If the lpAddress parameter is NULL, this value is rounded up to the next page boundary. Otherwise, the allocated pages include all pages containing one or more bytes in the range from lpAddress to (lpAddress+dwSize). This means that a 2-byte range straddling a page boundary causes both pages to be included in the allocated region.
flAllocationType
[in] Type of memory allocation. This parameter must contain one of the following values.

[shadeofgray]

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