Я буквально скопировал и вставил из предоставленного исходного кода для Numeric Recipes for C для разложения LU Matrix на месте, проблема в том, что он не работает.
Я уверен, что делаю что-то глупое, но был бы признателен любому, кто сможет указать мне правильное направление на это;Я работал над этим весь день и не вижу, что я делаю неправильно.
ОБНОВЛЕНИЕ ПОСЛЕ ОТВЕТА: Проект завершен и работает .Спасибо всем за руководство.
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MAT1 3
#define TINY 1e-20
int h_NR_LU_decomp(float *a, int *indx){
//Taken from Numerical Recipies for C
int i,imax,j,k;
float big,dum,sum,temp;
int n=MAT1;
float vv[MAT1];
int d=1.0;
//Loop over rows to get implicit scaling info
for (i=0;i<n;i++) {
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i*MAT1+j])) > big)
big=temp;
if (big == 0.0) return -1; //Singular Matrix
vv[i]=1.0/big;
}
//Outer kij loop
for (j=0;j<n;j++) {
for (i=0;i<j;i++) {
sum=a[i*MAT1+j];
for (k=0;k<i;k++)
sum -= a[i*MAT1+k]*a[k*MAT1+j];
a[i*MAT1+j]=sum;
}
big=0.0;
//search for largest pivot
for (i=j;i<n;i++) {
sum=a[i*MAT1+j];
for (k=0;k<j;k++) sum -= a[i*MAT1+k]*a[k*MAT1+j];
a[i*MAT1+j]=sum;
if ((dum=vv[i]*fabs(sum)) >= big) {
big=dum;
imax=i;
}
}
//Do we need to swap any rows?
if (j != imax) {
for (k=0;k<n;k++) {
dum=a[imax*MAT1+k];
a[imax*MAT1+k]=a[j*MAT1+k];
a[j*MAT1+k]=dum;
}
d = -d;
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j*MAT1+j] == 0.0) a[j*MAT1+j]=TINY;
for (k=j+1;k<n;k++) {
dum=1.0/(a[j*MAT1+j]);
for (i=j+1;i<n;i++) a[i*MAT1+j] *= dum;
}
}
return 0;
}
void main(){
//3x3 Matrix
float exampleA[]={1,3,-2,3,5,6,2,4,3};
//pivot array (not used currently)
int* h_pivot = (int *)malloc(sizeof(int)*MAT1);
int retval = h_NR_LU_decomp(&exampleA[0],h_pivot);
for (unsigned int i=0; i<3; i++){
printf("\n%d:",h_pivot[i]);
for (unsigned int j=0;j<3; j++){
printf("%.1lf,",exampleA[i*3+j]);
}
}
}
WolframAlpha говорит, что ответ должен быть
1,3,-2
2,-2,7
3,2,-2
Я получаю:
2,4,3
0.2,2,-2.8
0.8,1,6.5
И до сих пор я нашел по крайней мере 3 разные версии «одного и того же» алгоритма, поэтому я совершенно запутался.
PS да, я знаю, что для этого есть как минимум дюжина разных библиотек, но яменя больше интересует понимание того, что я делаю неправильно, чем правильный ответ.
PPS, поскольку в LU-декомпозиции нижняя результирующая матрица равна единице, а при использовании алгоритма Краутса, как (я думаю) реализованного, доступ к индексу массивавсе еще безопасно, и L и U могут быть наложены друг на друга на месте;следовательно, единственная результирующая матрица для этого.