Алгоритм Краута (нижне-верхняя (LU) декомпозиция матрицы)
ОГЛАВЛЕНИЕ
Если мы можем разложить матрицу линейной системы A в произведение A=L*U, где L - нижняя, а U - верхняя треугольные матрицы, то решение системы уравнений с произвольной правой частью производится весьма просто, применением двух обратных подстановок. Более того, в отличие от известного метода Гаусса-Жордана, разложенная матрица позволяет быстро решать серии линейных уравнений с различными правыми частями при одной и той же матрице.
Метод Краута позволяет провести LU-декомпозицию матрицы примерно за то же число операций, что и "прямая" часть метода Гаусса-Жордана. Итоговые коэффициенты двух треугольных матриц упаковываются в матрицу того же размера, что и A, и на том же месте в памяти; при этом верхняя матрица U размещается в наддиагональной части и на диагонали; нижняя L в поддиагональной части, а диагональные элементы L считаются все равными 1 (без ограничения общности) и не выводятся.
Алгоритм Краута представляет из себя следующее:
- Положить Lii=1 i=0,...,N-1 (здесь ничего не производится, а только имеется в виду для дальнейших шагов;
- Для каждого j=0,1,...,N-1 проделать:
- Для всех i=0,...,j вычислить Uij:
Uij=Aij - SUM(0<=k<i)(Lik*Uik)
(при i=0 сумму полагать нулем); - Для всех i=j+1,...,N-1 вычислить:
Lij = Aij - SUM(0<=k<j)(Lik*Uik))/Uii
Обе процедуры выполняются до перехода к новому j.
- Для всех i=0,...,j вычислить Uij:
Устойчивая работа алгоритма Краута возможно только при условии выбора ведущего элемента для каждого обращения к j из п.2 алгоритма. Выбор производится перед выполнением п.2 для очередного j перестановки необработанных строк матрицы так, чтобы строка, подставляемая на место j, содержала наибольший по модулю элемент в колонке j. Порядок перестановки необходимо запомнить, для чего достаточно использовать дополнительный вектор целых величин. Еще один вектор используется внутри алгоритма для масштабирования.
Алгоритм Краута имеет несколько приложений, одно из которых - решение системы линейных уравнений (обратная подстановка с учетом порядка перестановки строк), другие - вычисление обратной матрицы и вычисление детерминанта. Подробнее вызовы функций, связанных с алгоритмом Краута и реализация деталей алгоритма находятся в комментариях к программе и в ней самой.
Программа LU-декомпозиции
Программа самой декомпозиции и некоторых ее приложений приводится ниже.
/*
LU-декомпозиция в соответствии с алгоритмом Краута.
Описание:
int LU_decompos(double **a,int n,int *indx,int *d,double *vv);
параметры:
a - исходная матрица (n x n) на входе;
n - размер матрицы;
indx - массив целых чисел (размера n) для запоминания перемещений;
d - на выходе, содержит +1 или -1 для четного или нечетного числа перемещений.
vv - временный массив (size n).
Результат:
0 - некорректная входная матрица (непригодна для декомпозиции),
1 - положительный результат.
Обратная замена, использующая LU-декомпозицию матрицы.
Описание:
void LU_backsub(double **a,int n,int *indx,double *b);
Параметры:
a - матрица, разложенная по Крауту;
n - размер матрицы;
indx - порядок перемещения, полученный алгоритмом декомпозиции;
b - вектор (размера n).
Примечание: a и indx не изменяются в этой функции и могут быть использованы повторно.
Инвертирование матрицы с использованием LU-разложенной матрицы.
Описание:
void LU_invert(double **a,int n,int *indx,double **inv,double *col);
Параметры:
a - матрица, разложенная по Крауту;
n - размер матрицы;
indx - порядок перестановки, полученный алгоритмом декомпозиции;
inv - матрица назначения;
col - временный массив (размера n).
Получение матрицы, полученной с использованием LU-разложенной матрицы
Описание:
double LU_determ(double **a,int n,int *indx,int *d);
Параметры:
a - матрица, разложенная по Крауту;
n - размер матрицы;
indx - порядок перемещения, полученный алгоритмом декомпозиции;
d - знак равенства (+1 или -1) полученный при декомпозиции.
Результат:
определяющее значение.
*/
/* Для fabs(); включение может быть перемещено, если fabs() реализована как inline */
#include <math.h>
/* "small number" для избежания переполнения в некоторых случаях */
#define TINY 1.e-30
/* декомпозиция */
int LU_decompos(double **a, int n, int *indx, int *d, double *vv) {
register int i,imax,j,k;
double big,sum,temp;
for(i=0;i<n;i++) {
big=0.;
for(j=0;j<n;j++) if((temp=fabs(a[i][j]))>big) big=temp;
if(big==0.) return(0);
vv[i]=big;
}
/* главный цикл алгоритма Краута */
for(j=0;j<n;j++) {
/* это часть a) алгоритма, исключая i==j */
for(i=0;i<j;i++) {
sum=a[i][j];
for(k=0;k<i;k++) sum-=a[i][k]*a[k][j];
a[i][j]=sum;
}
/* инициализация для поиска наибольшего элемента */
big=0.; imax=j;
for(i=j;i<n;i++) {
sum=a[i][j];
for(k=0;k<j;k++) sum-=a[i][k]*a[k][j];
a[i][j]=sum;
if((temp=vv[i]*fabs(sum))>=big) {big=temp;imax=i;}
}
/* обмен строк, если нужен, смена равенства и размера множителя */
if(imax!=j) {
for(k=0;k<n;k++) {
temp=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=temp;
}
*d=-(*d); vv[imax]=vv[j];
}
/* сохраняем индекс */
indx[j]=imax;
if(a[j][j]==0.) a[j][j]=TINY;
if(j<n-1) {
temp=1./a[j][j];
for(i=j+1;i<n;i++) a[i][j]*=temp;
}
}
return(1);
}
/* обратная замена */
void LU_backsub(double **a,int n,int *indx,double *b) {
register int i,j,ip,ii=-1;
double sum;
/* первый шаг обратной замены */
for(i=0;i<n;i++) {
ip=indx[i];
sum=b[ip];b[ip]=b[i];
if(ii>=0) for(j=ii;j<i;j++) sum-=a[i][j]*b[j];
else if(sum) ii=i; /* получен ненулевой элемент */
b[i]=sum;
}
/* второй шаг */
for(i=n-1;i>=0;i--) {
sum=b[i];
for(j=i+1;j<n;j++) sum-=a[i][j]*b[j];
b[i]=sum/a[i][i];
}
}
/* инвертирование матрицы */
void LU_invert(double **a,int n,int *indx,double **inv,double *col) {
register int i,j;
for(j=0;j<n;j++) {
for(i=0;i<n;i++) col[i]=0.;
col[j]=1.;
LU_backsub(a,n,indx,col);
for(i=0;i<n;i++) inv[i][j]=col[i];
}
}
/* определяющее вычисление*/
double LU_determ(double **a,int n,int *indx,int *d) {
register int j;
double res=(double)(*d);
for(j=0;j<n;j++) res*=a[j][j];
return(res);
}