Matrizeninversion mit dem Gauss-Jordan-Algorithmus

• Aufrufbeispiel:

#include "matinvert.h"
.
.
.
void __fastcall TForm1::Button1Click(TObject *Sender)
{
    int iSize = 3;
    double Determinante;

    double Matrix[iSize][iSize]={{5,5,4},
                                 {9,3,2},
                                 {3,8,7}};

    // invertierte Matrix steht nach der Inversion wieder in 'Matrix'
    // zusätzlich wird die Determinante der Matrix zurückgegeben

    Determinante = invert_matrix(*Matrix, iSize);
}

• matinvert.h

#ifndef MATINVERT_H
#define MATINVERT_H

#define MAT_AT(A, COLW, J, I) (A[(J)*(COLW)+(I)])
double invert_matrix(double *A, int dim);

#endif

• matinvert.cpp

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "matinvert.h"

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 spaeter 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 */
}