#include<stdio.h>
#include<stdlib.h>
#include<math.h>
void inverse(double **A,unsigned int n, unsigned int *error)
{
unsigned int i,j,l,k;
double maxA,d,e;
unsigned int *M;
M=(unsigned int*)malloc(n*sizeof(unsigned int));
(*error)=0;
for(i=1;i<=n;i++)
{
// Czesciowy wybor elementu glownego wg wzoru (1.49)
maxA=0.0;
for(j=i;j<=n;j++)
{
d=A[j - 1][i - 1];
if(fabs(maxA)<fabs(d))
{
maxA=d;
k=j;
}
}
if(maxA==0.0)
{
(*error)=1;
free(M);
return;
}
/* Zapisywanie wskaznikow wierszy wystepowania elementu
ekstremalnego w i-tej iteracji w postaci wektora M[i] */
M[i - 1]=k;
A[k - 1][i - 1]=1.0;
for(j=1;j<=n;j++)
{
// Przestawienie i-tego wiersza z k-tym
d=(double)(A[k - 1][j - 1]/maxA);
A[k - 1][j - 1]=A[i - 1][j - 1];
A[i - 1][j - 1]=d;
}
// Generacja ciagu macierzy (1.42) wg wzoru rekurencyjnego (1.41)
for(j=1;j<=n;j++)
if(j!=i)
{
d=A[j - 1][i - 1];
A[j - 1][i - 1]=0.0;
for(l=1;l<=n;l++)
{
e=d*A[i - 1][l - 1];
A[j - 1][l - 1]-=e;
}
}
}
// Przestawianie kolumn macierzy zgodnie z wektorem wskaznikow M[i] (1.46)
for(i=n;i>=1;i--)
{
k=M[i - 1];
if(k!=i)
for(j=1;j<=n;j++)
{
d=A[j - 1][i - 1];
A[j - 1][i - 1]=A[j - 1][k - 1];
A[j - 1][k - 1]=d;
}
}
free(M);
}
int main(int argc, char **argv)
{
FILE *in;
FILE *out;
unsigned int i,j,n;
unsigned int error;
double **A;
if ((in = fopen(argv[1], "rt"))
== NULL)
{
fprintf(stderr, "Cannot open input file.\n");
return 1;
}
if ((out = fopen(argv[2], "wt"))
== NULL)
{
fprintf(stderr, "Cannot open output file.\n");
return 1;
}
while(!feof(in))
{
fscanf(in,"%d",&n);
A=(double**)malloc(n*sizeof(double*));
for(i=0;i<n;i++)
A[i]=(double*)malloc(n*sizeof(double));
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
fscanf(in,"%lf",&A[i - 1][j - 1]);
inverse(A,n,&error);
if(error==0)
{
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
fprintf(out,"%.12lf ",A[i - 1][j - 1]);
fprintf(out,"\n");
}
}
else
fprintf(out,"Macierz osobliwa \n");
for(i=0;i<n;i++)
free(A[i]);
free(A);
}
fclose(in);
fclose(out);
return 0;
}
Powyższy kod znalazłem w książce Bernarda Barona Metody numeryczne w Turbo Pascalu
(chociaż książka jest już niedostępna to kody są dostępne na stronie Helionu)
i przetłumaczyłem go na C
Znalazłem także w książce Knutha Sztuka programowania następujący algorytm
https://pdfhost.io/v/s7g7JnkFf_Knuth_inverse_matrix
Czy funkcja odwracania macierzy z kodów dołączonych do książki Bernarda Barona jest
zapisem algorytmu z książki Knutha ?
Jaki algorytm został tutaj przedstawiony