W
Fragt mich bloß nicht wo ich den Code her habe... Irgendwo mal Cut&Paste gemacht!!! Muss mit Sicherheit noch optimiert werden
#include <math.h>
#define MAT_AT(A, COLW, J, I) (A[(J)*(COLW)+(I)])
static double invert_matrix(double *A, int dim)
{
int i, j, k;
int *pivot_i, *pivot_j; // Ort des Pivot
double pivot; // das Pivot selbst
double tmp;
double det = 1.0;
pivot_i = (int*)malloc(dim*sizeof(int));
pivot_j = (int*)malloc(dim*sizeof(int));
for(k = 0; k < dim; ++k)
{
// k-tes Pivot finden
pivot = MAT_AT(A, dim, k, k); // suche init.
pivot_i[k] = k;
pivot_j[k] = k;
for(i = k; i < dim; ++i)
for(j = k; j < dim; ++j)
if(fabs(MAT_AT(A, dim, i, j)) > fabs(pivot))
{
pivot_i[k] = i;
pivot_j[k] = j;
pivot = MAT_AT(A, dim, i, j);
}
// aufmultiplizierte Pivots ergeben später die Determinante
det *= pivot;
if(det == 0.0)
{
// singulaere Matrix
free(pivot_i);
free(pivot_j);
return 0.0;
}
// Zeilen vertauschen, Vorzeichen beachten
i = pivot_i[k];
if(i != k) // wenn sich die Zeilen unterscheiden...
for(j = 0; j < dim; ++j)
{
tmp = -MAT_AT(A, dim, k, j);
MAT_AT(A, dim, k, j) = MAT_AT(A, dim, i, j);
MAT_AT(A, dim, i, j) = tmp;
}
// Spalten vertauschen
j = pivot_j[k];
if(j != k) // wenn sich die Spalten unterscheiden...
for(i = 0; i < dim; ++i)
{
tmp = -MAT_AT(A, dim, i, k);
MAT_AT(A, dim, i, k) = MAT_AT(A, dim, i, j);
MAT_AT(A, dim, i, j) = tmp;
}
// Spalte durch -pivot teilen
for(i = 0; i < dim; ++i)
if(i != k) // Pivot nicht einbeziehen
MAT_AT(A, dim, i, k) /= -pivot;
// Matrix reduzieren...
for(i = 0; i < dim; ++i)
{
tmp = MAT_AT(A, dim, i, k);
for (j = 0; j < dim; ++j)
if (i!=k && j!=k) // Pivot nicht einbeziehen
MAT_AT(A, dim, i, j) += tmp * MAT_AT(A, dim, k, j);
}
// Zeile durch Pivot teilen
for(j = 0; j < dim; ++j)
if(j != k) // Pivot nicht einbeziehen
MAT_AT(A, dim, k, j) /= pivot;
// Pivot durch seinen Kehrwert ersetzen
MAT_AT(A, dim, k, k) = 1./pivot;
}
// In einem weiteren Durchlauf Zeilen/Spalten vertauschen
for(k = dim - 2; k >= 0; --k)
{
i = pivot_j[k]; // mit Pivot-Spalte korrespondierende Zeilen
if(i != k) // wenn sich die Zeilen unterscheiden...
for(j = 0; j < dim; ++j)
{
tmp = MAT_AT(A, dim, k, j);
MAT_AT(A, dim, k, j) = -MAT_AT(A, dim, i, j);
MAT_AT(A, dim, i, j) = tmp;
}
j = pivot_i[k]; // mit Pivot-Zeile korrespondierende Spalten
if(j != k) // wenn sich die Zeilen unterscheiden...
for(i = 0; i < dim; ++i)
{
tmp = MAT_AT(A, dim, i, k);
MAT_AT(A, dim, i, k) = -MAT_AT(A, dim, i, j);
MAT_AT(A, dim, i, j) = tmp;
}
}
free(pivot_i);
free(pivot_j);
return det; // Determinante zurückgeben
};