|
Страница 2 из 2
Программа 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);
}
|