PARDISO

Итак, как было обещано, сегодня я прочитаю лекцию по одному из подразделов глав MKL, по одному из основных глав, а именно, PARDISO Parallel Direct Solver. По сути это решатель, который решает вам систему ax = b. Кто-то с ним уже работал? Один, два, немного человек. Отлично! Сегодня мы будем подробно разбирать, что он из себя представляет. Итак, это некая библиотека, как я уже сказал, решает Ax = b. Основное условие, чтобы матрица А была разреженная. Что значит разреженная, из терминов понятно: небольшое количество ненулевых элементов в каждой строке. Сама матрица может быть огромная, это не принципиально. Главное, чтоб ненулевых элементов в строке было 5, 7, 10, 20 немного. Откуда возникают подобные системы? Подобные системы возникают, когда вы моделируете некие уравнения физики, дифференциальные уравнения, уравнения теплопроводности и т. д. и т. п. методом конечных элементов, конечных объемов и т. п.

Итак, на основе чего работает PARDISO. По сути система простая и тривиальная. Сначала делается разложения матрицы А по Холесскому, то есть представляется матрица А в виде LDU-разложения, и после этого получившаяся система по очереди решается, то есть матрица U обращается легко (верхнетреугольная матрица), диагональная обращается еще легче. Матрица L обращается тоже несложно, если вы уже научились обращать верхнетереугольную матрицу. Есть одно но: даже если матрица изначально была разреженная и количество ненулевых элементов в ней мало, после того, как вы проведете разложение Холесского, результирующие матрицы L и U у вас получатся практически плотное. Соответственно, хранить их невозможно, работать с ними невозможно. Понятно, если изначально матрица была миллиард на миллиард, то для плотной матрицы миллиард на миллиард никаких не хватит ни оперативных памятей, ни жестких дисков, не говоря уже о том, чтобы с нею работать. Поэтому предлагается изначальную систему преобразовать некими стандартными способами. Эти способы широко известны в Интернете, они достаточно перспективно развиваются. Матрица А умножается справа и слева на ортогональные матрицы, для того чтобы привести к более удобному виду. Соответственно, все это реализовано в PARDISO. Как я уже сказал, что разложение Холесского в общем виде - это представление матрицы А в виде LDU, где U – нижнее треугольное, D – нижнее треугольное, D – диагональное или блочно диагональное. Понятно, что после такого разложения систему достаточно удобно решать. Как делается подобное разложение? Достаточно просто. Пусть у вас есть матрица А, представленная в общем виде А 1 1 А 1 2, и справа вы записываете матрицу L и матрицу L, транспонированную в общем виде. Сходу выписывается уравнение, что А 1 = L 1 в квадрате, отсюда мы находим L 1. второе уравнение А 1 2 = L 1 1, которая, мы уже знаем, на L 1 2, находим L 1 2 . И так далее и т. д. пока нам не надоест или пока мы не найдем все элементы матрицы L. Вот вкратце, как это действует, как этому учат в университете.

Как я уже сказал, если у нас матрица большая и разреженная, то если мы попробуем применить разложение Холесского вот до такой матрицы, то у нас получится верхнетреугольная и нижнетреугольная матрица плотная. Однако, если мы сделаем некое ортогонально преобразование – умножим вот на такую матрицу и на такую матрицу слева и справа, очевидно, что матрицы ортогональные, результирующая матрица получается значительно проще и видно, что L и L-транспонированные матрицы получились разреженные. Количество ненулевых элементов мало - успех!

далее, после того, как мы разложили матрицу А на LDU, или как в нашем случае на L LT- транспонированное, для того, чтобы найти решение системы уравнения Аx , а у нас это основная задача, мы разбиваем систему уравнений на две. Сначала решаем Ly=z, потом L-транспонорованное z = f соответственно в обратном порядке. Ну, как решается такая система уравнений - это совсем очевидно. Из первого уравнения мы находим y1, далее знаем y1, подставляем во второе уравнение, находим y2 и так далее находим y3. Все просто, все линейно, по сути, спустились по матрице сверху вниз.

Итак, с математикой практически закончили. Способы вызова тривиальные: всего одна функция, либо вы зовете ее из фортрана, либо вы зовете ее из си. Единственная проблема этой функции - несмотря на то, что она одна, у нее очень много различных параметров, каждый из которых делает что-то свое. Соответственно, сегодня мы попытаемся разобраться, что делают эти параметры, чтобы потом не было нужды спрашивать, что это такое. Опять-таки, все это описано в мануале, просто я попытаюсь рассказать это более по-человечески, более понятным языком.В мануале это тоже все описано, только вкратце, более научно, а я попытаюсь расшифровать.

Итак, самый первый параметр PARDISО: Pd, основной…структура. В фортране - это массив integer*8 длины 64, в C - это массив 64 void. По сути, это основная структура PARDISO, которая внутри раскрывается во все внутренние данные, внутренние массивы. Как там что происходит, знать не надо. Главное знать, что для того, чтобы вызвать PARDISO, нужно вызвать, вам нужно аллоцировать коротенький массив.

В PARDISO поддерживается достаточно удобная вещь, то есть, если вам необходимо решить несколько систем уравнений с различными матрицами, но что у матриц абсолютно одинаковая структура, такое встречается достаточно часто, такое встречается в методах декомпозиции, такое встречается в итерационных метода, когда вам нужно несколько матриц определенной структуры. Вы все можете подсунуть одному и тому же handle и, соответственно, параметр maxfct показывает, какое максимальное количество матриц вы подадите в PARDISO, какое максимальное количество матриц вы будете решать, и параметр mnum: с какой матрицей вы хотите работать в данный момент, c какой матрицей или уже факторизованной матрицей, в зависимости от того, что вам интересно. Или вам интересно факторизовать, или вам интересно уже решать систему уравнений.

Mtype: еще более важный параметр. В PARDISO поддерживается решение девяти типов матриц, то есть матрицы, которые вы подаете, могут быть симметричные, несимметричные, эрмитовы, комплексные, ортогональные. В зависимости от каждого варианта типа матриц у вас по-разному пойдет алгоритм, по-разному будет производительность алгоритма. Если вы подаете матрицу и говорите, что она симметричная, PARDISO будет пытаться делать разложение LLT, если вы говорите, что она несимметричная, то уже придется делать LDLT. Понятно, что так как алгоритм разный, в вы получите разное время на выходе и разную погрешность в зависимости от алгоритма. И что еще более важное, в зависимости от разного типа матрицы каждая матрица подается чуть-чуть в другом формате. В частности, если вы подаете матрицу и утверждаете, что она симметричная, вам не надо подавать всю матрицу. Вы подаете только ее верхнетреугольную часть, потому что нижнетреугольная ее симметрична. Зачем? Чтоб вам не тратить лишний объем данных.

Параметр iface. Параметр - это двузначное или однозначное число, которое может быть в случае двухзначного это, например, 1.1 , 1.2, 1.3, 2.2, 2.3 и 3.3. Каждая цифра значит, с какой стадии PARDISO вам надо стартовать, последняя цифра – какой стадией PARDISO вам надо закончить.

Стадия 11 - это значит вам интересно только провести так называемый реодеринг матрицы А, то есть вам надо умножить матрицу А, построить ортогональную матрицу к вашей изначальной матрице для того, чтобы для нее легко применялось разложение Холесского, чтобы результат был достаточно разреженный. Это стадия 11- реодеринг.

Стадия 22, стадия факторизации. Это та стадия, на которой уже результирующая матрица после реодеринга представляется в виде LDLT или LLT, или LDU, в зависимости от типа матрицы.

И, наконец, стадия 33, это стадия, которая показывает, что на данном этапе вы хотите именно решить систему уравнений уже с факторизованной матрицей, вы хотите получить решение системы уравнений. Соответственно, если вы вызываете фазу iface = 13, это значит, что вы хотите для данной матрицы сначала реодеринг, потом факторизацию, потом шаг решения без какого-либо прерывания. И, соответственно, последний вариант значения iface, iface=-1- освобождение памяти, деаллоцирование внутренних массивов и фактически окончание работы с PARDISO.

Параметры n, a, ia, ja. N – это размерность изначальной матрицы, тип у нее MKL INT. Массивы a, ia, ja – это тройка массивов, с помощью которой матрица А подается в CSR формате. PARDISO на данный момент поддерживает только CSR формат. Как он выглядит? Итак, как разреженная матрица А определяется в CSR формате с помощью массивов a, ia, ja. Массив ia: вы нумеруете все элементы подряд в матрице А по строчкам. То есть сначала вы нумеруете 1,2,3,4,5,6,7 – всего у вас 7 ненулевых элементов. Соответственно, массив ia показывает, с какого элемента начинается вполне конкретная строчка. Первая строчка всегда начинается с первого элемента, вторая строчка, в данном случае, начинается с третьего элемента, третья строчка – с пятого. Последний элемент - это всегда сумма ненулевых элементов плюс 1. Это 8. Это ia. Соответственно, ja массив значительно более простой, он показывает, в каком столбце находится данный элемент. То есть если вас интересует, в каком столбце находится четвертый элемент в матрице А, вы смотрите на ja четвертое и узнаете, что он находится во втором столбце. Ну и, соответственно, массив А - это уже массив дабл, сингл или комплекс, в котором записаны подряд все вот эти элементы.

Входной параметр perm - это вектор перестановок. То есть, как вариант, если вас интересует, на какую ортогональную матрицу, на какую матрицу перестановок PARDISO умножил изначальную матрицу, вы можете получить данные с помощью массива perm, в данном случае, в мануале более подробно описано, как и что чему соответствует. Либо вы можете PARDISO подать свой вариант перестановок, то есть задаете Перм в некоем стандартном типе, подаете ему, он говорит отлично. И, соответственно, факторизацию, разложение Холесского делает именно для вашей матрицы перестановок. Так же массив perm служит еще для кое-каких целей, на них я остановлюсь более подробно позже

Nrhs: PARDISO поддерживает, в частности, решение Ax=b  не только для одной правой части, но и для многих правых частей. Если вам надо решить одну и ту же систему уравнений, где правых частей у вас сто, соответственно, подаете nrhs=100, и PARDISO делает один раз реодеринг для данной матрицы, один раз разложение (факторизацию) и после этого решает более эффективно 100 уже факторизованных систем. Это будет эффективнее, чем сто раз вызывать PARDISO для каждой конкретной правой части.

Массивы b и x. Ну, здесь вообще все просто. Массив b – массив правых частей, если у вас их много, то массив размерности nrhs*N, где n – размерность матрицы. x – это вектор решения точно такой же размерности, в котором будет выводиться решение. Или не будет, в зависимости, опять-таки, от параметров.

Наконец, самый понятные параметры error и msglvl. msglvl это параметр, который говорит, хотите ли вы получать от PARDISO некую статическую информацию, в частности, число ненулевых элементов в матрице А, число ненулевых элементов получившееся в факторизованной матрице А, сколько времени было потрачено на факторизацию и тд. Error: вы правы, что если 0, то все отработало отлично, не совсем правы, что если не 0, то все отработало не отлично. Это значит, PARDISO вам говорит какую-то вполне определенную информацию, что произошло. Например, вы подали матрицу с очень большим числом обусловленности, но PARDISO вам возвращает информацию, что число обусловленности большое, и, может быть, поэтому некорректно получилось. Вариантов значений error достаточно много, они все описаны в MKL-мануале. Если вы какой-то получили, то можно там посмотреть, какой чему соответствует.

И наконец, самый массив iparm. Он длины 64. Каждый элемент это массива = параметр, либо входной, либо выходной. Фактически, это некие рукоятки, изменяя которые, вы можете изменять работу PARDISO, и через которые вам возвращаются некие результаты работы. Например возвращается, в какой строке матрицы произошла какая-то ошибка, или что матрица, которую вы подали, не в csr формате. Разберем это еще более подробно. Итак, iparm[0]. Если вы не хотите возиться с этими параметрами, вникать в них, выставляете iparm[0] = 0 и ни про что остальное не беспокоитесь. PARDISO само вам выставит все остальные параметры, и вы как-то сможете работать. Это самый простой вариант.

Iparm [1]: В PARDISO поддерживается 3 варианта реодеринга. 3 варианта, каким способом PARDISO находит матрицу вращений перестановок для изначальной матрицы. Это либо minimum degree algorithm (алгоритм минимальных степеней), либо алгоритм вложенных сечений, либо параллельный алгоритм вложенных сечений, то же самое, только реализованное для процессоров с общей памятью. По умолчанию используется не параллельный алгоритм вложенных сечений. Почему? Потому что не смотря на то, что алгоритм параллельных сработает быстрее, результирующая матрица может получиться немного отличной, другой. И не факт, что она будет лучше или хуже, она будет другой. А у вас производительность факторизации шага 22 и шага 33 от этого резко зависит. Мы привыкли использовать iparm[1]=3 для своих локальных нужд, но вы можете попробовать, попробуйте все варианты, какой больше понравится. В зависимости от исходной матрицы могут быть эффективными разные алгоритмы.

Итак, iparm[3]. Я рассказываю про основные параметры. Итерационное уточнение: кроме того, что PARDISO это прямой решатель, достаточно часто оказывается, что в результате большого числа обусловленности или в результате погрешностей вычислений, если у вас собственные числа у изначальной матрицы были близки к нулю еще что-то, у вас на выходе невязка оказывается достаточно большая. В данном случае вы можете включить итерационное уточнение, и тогда полученное разложение Холесского PARDISO использует внутри как предобуславливатель для некоего итерационного процесса. Он досчитывает достаточно быстро, и невязка получается лучше. То есть вы попробовали запустить PARDISO, посмотрели, невязка вас не устроила, включили итерационное уточнение, получили ответ лучше. Как он определяется? Он определяется в виде 10 L +K, где К – это какой вид итерационного уточнения, а L – это до какой степени невязки надо гнать.

Iparm[4]. Управление матрицей перестановок. Если он 0, то вас не интересует эта матрица перестановок, что достаточно часто бывает, т.к. вас интересует чаще результирующий ответ. Либо, если он равен 1, то вы сами подаете матрицу перестановок, либо, если он равен 2, то вы просите, чтобы PARDISO вам матрицу перестановок вернул. Не смотря на то, что PARDISO поддерживает факторизацию и шаг решения матрицей перестановок от пользователя, я крайне не рекомендую этот способ. Потому что обычно без этого эффективность на порядки лучше. То есть, если вы совсем уверены, что ваша перестановка сделана успешно, тогда попробуйте, но скорее всего нет.

Iparm[5]. Вариант, куда сохранять решения. Либо вы сохраняете решения в вектор х, либо вы решения сохраняете поверх вектора b, что тоже часто нужно. Например, в итерационных методах, когда вы умножаете на предобуславливатель, если к предобуславливателю вы обращаетесь с помощью PARDISO, то, по сути, вам в этот вектор надо вернуть поправку, но необходимо отметить, что все равно вектор х нужен PARDISO для внутренних вычислений. То есть, в любом случае, его аллоцировать придется. По сути, iparm[5] это переключатель, который говорит, куда вам надо сохранять решение.

Iparm[7]. Тоже понятно. Сколько максимально итераций можно сделать в итерационном уточнении. Что необходимо отметить, если этот параметр больше 0, то делается заданное количество итераций в той же точности, если меньше, то итерационное уточнение делается в двойной точности. То есть если вы PARDISO вызывали к системе в одинарной точности, то у вас итерационное уточнение будет в дабл, если вы PARDISO вызывали в дабл точности, то итерационное уточнение у вас пойдет в quad precision. Это будет дольше, но есть шанс получить ответ точнее. Опять таки сами выбираете, что вы хотите.

Массивы и параметры iparm[9], iparm[10] и iparm[12], массивы, которые работают с малыми числами, возникающими на диагонали матрицы при факторизации. Если вы подадите знако-неопределенную матрицу, у которой есть как положительные, так и отрицательные собственные числа, вы имеете хороший шанс встретить ноль на диагонали при факторизации. Поскольку вы встретили ноль, сразу возникает проблема. Соответственно, вот эти параметры iparm[9], iparm[10] и iparm[12], они отвечают за то, чтобы этой ситуации как-то избежать. Но, опять-таки, как избежать? Вы запустили PARDISO, которое вернуло Вам error не равное 0, равное "мы встретили ноль на диагонали.", после чего вы один из этих параметров меняете. То есть я могу описать, объяснить что каждый из этих параметров делает, но в любом случае все это надо делать экспериментально. Запустили с одними параметрами – не получилось получить решение, прогнали с другими параметрами - получилось.

Iparm[17], если вы подадите его отрицательным, то после этапа факторизации он выдаст вам количество ненулевых элементов в разложении Холесского. Это вам даст возможность посмотреть, сколько памяти у вас потребовало PARDISO, потому что основная память, которая уходит, это количество ненулевых элементов в разложении Холесского или в LU-разложении.

Ну, и Iparm[18] вернет вам приблизительные мегафлопы после этапа реодеринга. Зачем нужны параметры iparm[17]- iparm[18]. Вы знаете размер вашей оперативной памяти, провели этап реодеринга, посмотрели, оказалось, что число ненулевых элементов у вас точно не влезает в оперативную память. Значит надо что-то менять. Чтобы не гнать факторизацию и по пути не получать "у вас кончилась память".

Iparm[20]. Тут надо остановиться поподробнее. Если iparm[20] равен нулю, то все в разложении Холесского, как мы привыкли: один элемент - одно число, считаем и считаем. Если iparm[20] равен 1, то используется так называемый Баунч-Кафман-алгоритм, когда в разложении Холесского мы встречаем 0, то мы говорим: на самом деле, элемент нашего диагонального разложения - это не элемент, а блок 2 на 2. И тогда понятно: у нас проблема не в элементе нулевом, а если определитель этой матрицы 2 на 2 нулевой, что шансов меньше. И тогда наша матрица D в разложении LDU становится не совсем диагональной, а блочно диагональной, у которой на диагонали блоки 1x1 и 2x2. Баунч-Кафман-алгоритм используется по умолчанию, как более стабильный, но производительность чуть-чуть отличается от случая, когда Баунч-Кафман-алгоритм отключен, причем нельзя сказать в какую сторону. И, опять-таки, стабильность отличается. Попробуйте, что больше понравится.

Iparm[23] - одна из наших последних разработок. Это параметр, который говорит, по какому алгоритму надо факторизовать матрицу. Либо по стандартному алгоритму, достаточно хорошо описанному в большом количестве публикаций, либо по так называемому Two-Level Factorization, или двухслойный алгоритм. И он действительно показывает значительное улучшение по сравнению со стандартным алгоритмом на многих матрицах, но только для большого числа ядер на процессоре. Если у вас процессор 8-ядерный, 16-ядерный или 32-ядерный, включаете двухслойный алгоритм, и, скорее всего, производительность факторизации у вас будет лучше.

Iparm[26]: Этот параметр проверит, действительно ли CSR формат, который вы подали, корректный. Он проследит, пронумеровали ли вы все элементы матрицы последовательно, в случае симметричной матрицы подали ли вы только верхний портрет. И еще, несколько тонкостей: что в массиве iа-ja нет отрицательных элементов, и т.д. Подобные тонкости PARDISO проверит самостоятельно в случае iparm[26]=1, и если у вас там ошибка, PARDISO выйдет с кодом ошибки и скажет, (если msglvl=1 ): вот эта строка в CSR формате не корретна. Но, опять-таки, каждый раз включать этот режим не стоит, потому что на него тоже тратится какое-то время, соответственно, производительность PARDISO уменьшается.

Iparm[27] и iparm[34]. ну, здесь все просто. Если iparm[27] = 0 или 1, то у вас либо дабл, либо сингл входные данные. И если у вас iparm[34] = 0, то у вас фортран стиль входных данных, то есть первый элемент массима ia равен 1, ну и все столбцы матрицы а нумеруются с единицы. Если iparm[34]=1, то у вас с стиль входных данных. Первый столбец - это нулевой по счету столбец, первая строка - это нулевая по счету строка, ia[0] = 0.

Iparm[30] – так же недавняя разработка. Допустим вас не интересует все решение система Ax=b, а интересует только последняя компонента. Или если вы точно знаете, что у вас в правой части все нули, кроме 1-го, 2-х, 3-х элементов, то вы с помощью массива perm и переключателя iparm[30] можете объяснить это PARDISO. И тогда у вас сильно увеличится эффективность этапа 33 (solving step). Но что надо понимать: у вас резко может упасть эффективность этапов 11 и 22., то есть это имеет смысл делать только в случае, когда вам действительно нужно найти несколько компонен вектора решения. Или вам надо решить очень большое количество систем с одной и той же матрицей много раз, когда только несколько компонент правой части не ноль. Например, вы решаете систему уравнений Ax=b для функции Грина.

Iparm[59] – ooc алгоритм. Если вам не хватает оперативной памяти, то PARDISO предлагает решение использовать жесткий диск, использовать его достаточно эффективно. То есть PARDISO самостоятельно работает с жестким диском: сбрасывает туда файлы, читает оттуда файлы. Все это делается быстро и эффективно. Называется этот алгоритм OOC – алгоритм. Есть некий набор переменных окружения, с помощью которых вы можете задать, какой объем оперативной памяти вы хотите использовать. Куда вы хотите, чтобы PARDISO писало свои рабочие временные файлы, хранить эти файлы потом или не хранить и т.д.. Но, в любом случае, с помощью этого алгоритма вы можете решать значительно большие задачи, чем вам позволяла бы оперативная память. Допустим, оперативной памяти у вас всего 4 Гб, а жесткого диска Гб 500 найдется, это значит, можно решать матрицы значительное большего размера. Но надо понимать, что этот алгоритм менее эффективный, чем если бы у вас вся память требуемая PARDISO поместилась в оперативную память. Но, тем не менее, иногда и достаточно часто приходится с этим жить.

Как я уже сказал, есть интерфейс все один интерфейс - PARDISO. Но если вы привыкли работать С++, то PARDISO можно вызывать с помощью DSS интерфейса. DSS расшифровывается как Direct Sparse Solver. Допустим, если вы хотите решить систему уравнений с уже факторизованной комплексной матрицей, то вы можете вызвать просто вызвать такой вариант DSS. Набор этих интерфейсов достаточно большой. Всем основным вызовам PARDISO соответствует DSS интерфейс, но надо помнить: если вам надо решить комплексную матрицу, вы можете вызвать DSS_SOLVE_COMPLEX, если вы там хотите сделать еще и итерационное уточнение, задать еще и в C-формате, то такое DSS не поддерживается, вам придется вызывать PARDISO. То есть если вы хотите покрутить что-то в PARDISO, поисследовать, полюбопытствовать, как-то улучшить производительность, то вам нужно работать через PARDISO интерфейс.