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

До сих пор я работал со СЛАУ сравнительно небольших размерностей. Матрицы таких СЛАУ можно было размещать в оперативной памяти.
Матрицы плотно заполнены.
Вообще говоря я пока еще не знакомился с теорией, относящейся к данному вопросу, поскольку, как я уже писал, не было такой необходимости, а в университете я этого не проходил.
Пока у меня есть только два решения.
Первое - записывать получаемые массивы в файл, то есть размещать их на жестком диске, и запаршивать каждый элемент когда это необходимо. Тогда это может очень увеличить время счета.
Второе - вообще не считать коэффициенты матрицы до тех пор, пока соответствующий элемент не затребуется на определнном шаге алгоритма, решающего СЛАУ, и считать его только в этот момент. Тогда придется один и тот же элемент определять несколько раз, да и время счета, по-моему, не сильно будет отличаться от первого варианта (или я ошибаюсь?).
В общем, вопрос задан в названии темы. Прошу поделиться опытом.

Telc, а какая это размерность, что в оперативку не влазит?

Так, если даже 10000 уравнений, то 10^10 размерность... Вроде так-то немного, а вот "проклятие размерности".
В файл сохранять - это не выход, по-моему. Благо виртуальная память на то и придумано, чтоб файл подкачки раздувать. Ставьте 64-битную операционку, жёсткий диск побольше и должно влазить в нёё-то.

Порядка 50000 уравнений (матрица квадратная, плотно заполненная).

Круто. С наскока могу лишь предолжить опять 64-битную систему для решения сей проблемы, хотя наверняка есть более оптимальные способы.

В windows создаётся страница памяти нужного размера и в ней храните свой мегамассив. ОС самостоятельно будет принимать решение о выгрузке/загрузке страницы

Цитата AndreyK @ 05.11.07, 19:56

В windows создаётся страница памяти нужного размера и в ней храните свой мегамассив. ОС самостоятельно будет принимать решение о выгрузке/загрузке страницы

А как это сделать? Этот вариант является наиболее оптимальным?

Сообщение отредактировано: Telc - 05.11.07, 20:39

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


Telc
Массив на диск. Доступ по средство чатения. В память должно считываться только необходимые данные. Алгоритм обработки должен снизить число обращений к диску.

Решил полной ответ написать.
Для доступа к файлу есть несколько механизмов. Просто читать файл по частям в память или целеком.

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

Теперь о том почему так не сделать.
Во-первых в 32-битных системах есть ограничения на размер отоброжаемого файла.
Ограничение это то что пользовательскому процессу дается 2ГБ виртуальной памяти, в некоторых случаях 3ГБ. Но она разбита на куски. Так что максимальный размер ограничивается размера наибольшего куска.
Так что тут выход 64-битная ОС.

Во-вторых Если физическая память исчерпона, то программа будет просто ужастно тормазить. Загружая/выгружая страницы памяти.

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

Сообщение отредактировано: Pavia - 05.11.07, 21:15

Цитата Telc @ 05.11.07, 14:37

Порядка 50000 уравнений (матрица квадратная, плотно заполненная).

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

Нет, обойти такую проблему принципиально не получится.
Вопрос по сообщению nvm: как использовать?
Вопрос по сообщению Pavia: что будет работать быстрее: записать на диск массив и затем, по мере необходимости считывать тот или иной элемент, или вообще не запоминать массив (ни на диске ни в оперативке), а считать конкретный элемент непосредственно перед его использованием (например, надо умножить a_{i,j} на b_{j} я вызываю процедуру, которая мне посчитает этот элемент a_{i,j} и тогда я перемножу его на b_{j})?

Telc
Если ты будешь считывать из файла один элемент, то это медленно. Лучше так прочитал строчку умножил ее на элимент сохронил. Или считал две строчки произвел между ними операции сохронил результаты.

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

Telc
50000x50000x8=18,62 Гб. Тебе нужен хороший многопроцессорный сервер с как минимум 32 Гб оперативки и 64-х разрядная операционка. Если памяти не хватит, то время решения (и так немаленькое) увеличится еще на порядок даже при условии оптимального чтения с диска (то есть все строки уравнения хранятся на диске последовательно и обрабатываются тоже последовательно). Я вообще не понимаю, как ты собираешься решать такую систему, понадобятся очень эффективные алгоритмы (метод Гаусса слишком медленный).

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

Добавлено 06.11.07, 08:11
Предположим, что система решается методом Гаусса => это 2x50000 проходов по матрице, пусть скорость чтения с диска - 120 Мб/сек (рэйд там какой-нибудь). 18.62*1024*2*50000/120/86400=184 дня непрерывной работы. Так что нужно или долго ждать или использовать более эффективный метод...

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

Цитата Pacific @ 06.11.07, 07:53

Telc
50000x50000x8=18,62 Гб. Тебе нужен хороший многопроцессорный сервер с как минимум 32 Гб оперативки и 64-х разрядная операционка. Если памяти не хватит, то время решения (и так немаленькое) увеличится еще на порядок даже при условии оптимального чтения с диска (то есть все строки уравнения хранятся на диске последовательно и обрабатываются тоже последовательно). Я вообще не понимаю, как ты собираешься решать такую систему, понадобятся очень эффективные алгоритмы (метод Гаусса слишком медленный).

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

Добавлено 06.11.07, 08:11
Предположим, что система решается методом Гаусса => это 2x50000 проходов по матрице, пусть скорость чтения с диска - 120 Мб/сек (рэйд там какой-нибудь). 18.62*1024*2*50000/120/86400=184 дня непрерывной работы. Так что нужно или долго ждать или использовать более эффективный метод...

Я вот не понимаю, зачем человеку 32 ГБ оперативки. Таких компов днём с огнём не сыскать. Вот, например, взять 32-битную винду. Там можно захапать теоретически 4 ГБ оперативки, но никто же не говорит, что нужно, чтоб все 4 ГБ были вставлены в комп. Этим и занимается виртуальная память, раздувая оперативку до максимально возможного значения за счёт файла подкачки.

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

experimenter
Если матрица не влезет в память, то все упрется в скорость считывания с диска => полгода на решение одной системы. Если в память влезет, то все упрется в ПСП памяти, а это уже гораздо более разумное время счета (одна-две недели). Так что ему нужен такой комп (а они продаются, вот только стоят соответственно).

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

Если матрица не влезет в память, то все упрется в скорость считывания с диска => полгода на решение одной системы. Если в память влезет, то все упрется в ПСП памяти, а это уже гораздо более разумное время счета (одна-две недели). Так что ему нужен такой комп (а они продаются, вот только стоят соответственно).

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

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

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

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

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

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

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

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

huge system of linear equations

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Цитата 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 строки матрицы на каждой итерации, что совсем немного, и обычная локалка это легко потянет.

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

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

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

Попробуйте покопать в сторону блочно-матричных алгоритмов (на таких алгоритмах основан 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
По поводу двух стульев. Когда я говорил про операции с матрицами, помещающимися в кэш, то имел в виду, что обновление остатка также можно оптимизировать и реализовать в блочно-матричном стиле. Но, подумав, понял, что на самом деле здесь не тот случай, и определяющим фактором будет скорость работы с диском. Так что стул только один.

Цитата 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.

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

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

Не уверен, что это подойдет - ведь предыдущие строки наверняка надо хранить в памяти, на каждой новой строке проходить по ним.

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

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

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

Еще про память.
Пусть мы имеем m систем n-го порядка с одинаковой матрицей, то есть расширенная матрица имеет размер n*(n+m)
При решении системы получается последовательность матриц размеров:
1*(n+m-1), 2*(n+m-2), 3*(n+m-3), ..., (n-1)*(m+1), n*(m), которые могут быть размещены по одному и тому же адресу (при построчном порядке, как в C). Последняя матрица будет содержать решения.

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

Цитата experimenter @ 06.11.07, 08:25

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

Гм... Ну зачем же покупать? можно арендовать.
Мне стало интересно - порылся по теме в инете и нашол некоторые цифры:

Цитата с какого то форума

... на западе, там стоимость одного гигагерца маширнного времени за час: 1 доллар, у нас (в России) 0,15 - 0,20, то есть в 5-7 раз меньше...

Сообщение отредактировано: Alexus - 08.11.07, 05:10

Цитата Alexus @ 08.11.07, 05:08

Гм... Ну зачем же покупать? можно арендовать.
Мне стало интересно - порылся по теме в инете и нашол некоторые цифры:

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

В любом случае, прежде необходимо отаработать сам алгоритм на простом ПК.

Цитата Telc @ 08.11.07, 12:45

В любом случае, прежде необходимо отаработать сам алгоритм на простом ПК.

А до этого лучше на бумажке!

Не влезая в дебри теории, личн я бы подошёл к проблеме чисто технически

Преобразовал бы массив в список - т.е. паре (i, j) сопоставляется значение. И есть небольшой классик, который в себе умеет держать набор этих значений. Да, матрицы не разреженные, экономии это не даст
Зато есть возможность контролировать размер занимаемой памяти - (sizeof(i)*2 + sizeof(value))*current_loaded_items
Ну а далее уже можно тупо считывать с диска все записи, пока забьём отведенный размер. А после - считывать по надобности, высвобождая самые "старые" значения. Дешёво и сердито. И время на реализацию - 1-2 часа

Нет я конечно могу понять желание решить в лоб 50000 уравнений, но даже если удастья представить в памяти такой объём данных, то врядли стоит завидовать человеку, которому предстоит дождаться решения данной системы, к тому же хранить только те данные которые могут быть вычеслены это ещё более нагрузка на ЦП, да и где гарант, что система окажется совместна или не будет имет бесконечное множество решении. В жизни не видел ситуации где бы могло понадобится такое количество корней для чего либо. С другой стороны использовать жесткий диск конечно можно, но тогда необходимо минимизировать количество обращений к нему. Например если для решения системы используется всеми горячо любимый метод Гаусса, то можно хранить отдеьные строки (в количестве 2 штук) на прямом и обратном ходе алгоритма, после этого, в случае если система является невырожденной(по сути каждой переменной соотвестует едиснтвенное решение), то необходимо, считать только диагональные элементы, в противном случае лучше сразу бить в барабаны, потомучто даже если можно будет наити аналитическую форму записи для свободных переменных, то анализ такой системы в большинстве случаем сведётся к поиску оптимума, а вот эта задача скорее всего уже приведёт к ожиданию конца вычислении, ну как минимум до пенсии уважаемого Telc. Другого способа просто нет.

Блин. Ну почему постоянно у всех есть желание изобретать велосипед? Блочно-матричные алгоритмы - самое простое и естественное решение проблемы обращений к жесткому диску. Не надо никаких списков (для линейной алгебры вообще смерть - замедление в десятки раз), не надо ничего вообще. Матрица обрабатывается не по строкам, а по блокам. Соотношение между размером задачи и количеством доступной оперативной памяти естественным образом преобразуется в размер обрабатываемого блока. Число считываний элемента матрицы с жесткого диска не N (как в обычном алгоритме), а N/K, где K - размер блока.

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