lup разложение

Ответ от МИР[гуру]
Обратная матрица — такая матрица A-1, при умножении на которую исходная матрица A даёт в результате единичную матрицу E: AA − 1 = A − 1A = E
Использование LU/LUP-разложения
Если матрица A невырождена то для нее можно рассчитать LUP-разложение PA = LU. Пусть PA = B, B − 1 = D. Тогда из свойств обратной матрицы можно записать: D = U − 1L − 1. Если умножить это равенство на U и L то можно получить два равенства вида UD = L − 1 и DL = U − 1. Первое из этих равенств представляет собой систему из n² линейных уравнений для frac{n(n+1)}{2} из которых известны правые части (из свойств треугольных матриц) . Второе представляет также систему из n² линейных уравнений для frac{n(n-1)}{2} из которых известны правые части (также из свойств треугольных матриц) . Вместе они престаляют собой систему из n² равенств. С помощью этих равенств можно реккурентно определить все n² элементов матрицы D. Тогда из равенства (PA)-1 = A-1P-1 = B-1 = D. получаем равенство A-1 = DP.
Ниже представлен код на языке C++ алгоритма обращения матрицы с помощью LUP-разложения с комментариями. Matrix Inverse(const Matrix &A) {
//n - размерность квадратной матрицы A
const int n=A.Rows();
//при инициализации задается размерность nxn
Matrix X(n, n);
Matrix P(n, n);
Matrix C(n, n);
//предполагается что в результате следующего вызова матрица C = L + U - E
LUP(A, C, P);
for(int k = n-1; k >= 0; k--) {
X[ k ][ k ] = 1;
for( int j = n-1; j > k; j--) X[ k ][ k ] -= C[ k ][ j ]*X[ j ][ k ];
X[ k ][ k ] /= C[ k ][ k ];
for( int i = k-1; i >= 0; i-- ) {
for( int j = n-1; j > i; j-- ) {
X[ i ][ k ] -= C[ i ][ j ]*X[ j ][ k ];
X[ k ][ i ] -= C[ j ][ i ]*X[ k ][ j ];
}
X[ i ][ k ] /= C[ i ][ i ];
}
}
X = X*P;
return( X );
}
Советую специальную программу по нахождению обратных матриц LU-методом - программа ДЕЛФИ
И есть хороший сайт в подмогу: Источник: матрица