Русская Википедия:LUP-разложение
LUP-разложение (LUP-декомпозиция) — представление данной матрицы <math>A</math> в виде произведения <math>PA = LU</math> где матрица <math>L</math> является нижнетреугольной с единицами на главной диагонали, <math>U</math> — верхнетреугольная общего вида, а <math>P</math> — т. н. матрица перестановок, получаемая из единичной матрицы путём перестановки строк или столбцов. Такое разложение можно осуществить для любой невырожденной матрицы. LUP-разложение используется для вычисления обратной матрицы по компактной схеме, вычисления решения системы линейных уравнений. По сравнению с алгоритмом LU-разложения алгоритм LUP-разложения может обрабатывать любые невырожденные матрицы и при этом обладает более высокой вычислительной устойчивостью.
Алгоритм LUP-разложения
Пусть <math>A=(a_{ij})</math>, <math>L=(l_{ij})</math>, <math>U=(u_{ij})</math>. На практике как правило вместо матрицы перестановок P используют вектор перестановок получаемый из вектора <math>p_{i}=(1, 2, ..., n)</math> путём перестановки элементов соответствующих номерам строк переставляемых в матрице P. Например, если
<math>P = \begin{pmatrix}
0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\
\end{pmatrix} </math>
то <math>p = (2, 1, 3)</math> так как матрица P получена путём перестановки первой и второй строки. Вычисление LUP-разложения ведётся в несколько шагов. Пусть матрица C = A. На каждом i-м шаге сначала производится поиск опорного (ведущего) элемента — максимального по модулю элемента среди элементов i-го столбца, находящихся не выше i-й строки, после чего строка с опорным элементом меняется местами с i-й строкой. Одновременно производится такой же обмен в матрице P. При этом, если матрица невырождена, то опорный элемент гарантированно будет отличен от нуля. После этого все элементы текущего i-го столбца, находящиеся ниже i-й строки, делятся на опорный. Далее из всех элементов <math>c_{jk}</math> находящихся ниже i-й строки и i-го столбца (то есть таких что j>i и k>i) вычитается произведение <math>c_{ji}c_{ik}</math>. После этого счётчик i увеличивается на единицу и процесс повторяется пока i<n где n — размерность исходной матрицы. После того как все шаги будут выполнены матрица C будет представлять собой следующую сумму:
<math>C = L+U-E</math>
где E — единичная матрица.
В алгоритме используется три вложенных линейных цикла так что общую сложность алгоритма можно оценить как O(n³).
Реализация алгоритма на языке C++
Ниже представлен программный код приведённого выше алгоритма на языке C++. Здесь Matrix — некоторый контейнер, поддерживающий операцию индексирования. Обратите внимание, что отсчёт ведётся с нуля, а не с единицы.
void LUP(const Matrix &A, Matrix &C, Matrix &P) {
//n - размерность исходной матрицы
const int n = A.Rows();
C = A;
//загружаем в матрицу P единичную матрицу
P = IdentityMatrix();
for( int i = 0; i < n; i++ ) {
//поиск опорного элемента
double pivotValue = 0;
int pivot = -1;
for( int row = i; row < n; row++ ) {
if( fabs(C[ row ][ i ]) > pivotValue ) {
pivotValue = fabs(C[ row ][ i ]);
pivot = row;
}
}
if( pivotValue != 0 ) {
//меняем местами i-ю строку и строку с опорным элементом
P.SwapRows(pivot, i);
C.SwapRows(pivot, i);
for( int j = i+1; j < n; j++ ) {
C[ j ][ i ] /= C[ i ][ i ];
for( int k = i+1; k < n; k++ )
C[ j ][ k ] -= C[ j ][ i ] * C[ i ][ k ];
}
}
}
}
//теперь матрица C = L + U - E
Литература