Der folgende Source berechnet nxn-Matrizen nach der Gauss-Elimination mit Pivot-Suche. Dadurch ist dieser Algorithmus sehr stabil.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <exec/types.h>

/* Alles von 0 weg! n = 1..max */

void pmat(double **mat, int n)
{
  int i,j;
  for (j=0; j<n;j++)
  {
    for (i=0; i<=n;i++)
    {
      printf("%f ",mat[j][i]);
    }
    printf("\n");
  }
    printf("\n");
}

int pivot(double **mat, int n, int zeile) /* Ab Zeile wird das Pivot-Element auf 0 Kontrolliert und ev. Zeile getauscht */
{
  int i,j;
  double help;

  if (mat[zeile][zeile] != 0.0)
    return(0);                  /* alles OK */

  for (i=zeile+1;i<n;i++)         /* untere Zeile anschauen */
  {
    if (mat[i][zeile] != 0.0)
    {
      for (j=0;j<=n;j++)   // mit Lösungsvektor kopieren
      {
	help = mat[zeile][j];
	mat[zeile][j] = mat[i][j];
	mat[i][j] = help;
      }
      return(1);                  /* Zeile mit einer unteren getauscht */
    }

  }
  return(-1);
}

int gauss(double **mat, double *erg, int n)
{
  int i,j,k;
  double q;

  for (j=0;j<n;j++)
  {
    pivot(mat,n,j);  /* Ist erste Zeile ok? */
    // Neue Eliminationszeile berechnen */
    for (i=j+1;i<n;i++)
    {
      q = -1 * (mat[i][j] / mat[j][j]);   // Durch Linkes oberes dividieren */
      for (k=j;k<=n;k++)  // Zeilenweise neue Elemente berechnen, incl. Lösungsvektor
      {
	mat[i][k]= mat[i][k] + (mat[j][k] * q);
      }
    }

  }

  for (i=n-1;i>=0;i--)
  {
    for (j=i+1;j<n;j++)  // Von links nach Rechts korrigieren
    {
      mat[i][n] = mat[i][n] - (mat[i][j] * erg[j]);
    }
    erg[i] = mat[i][n] / mat[i][i]; // Ergebnis eintragen
  }
  return(0);
}

#define DIM 5

int main(int argc, char *argv[])
{
  int i;
  double *mat[DIM];  /* Zeile,Spalte, Lösungsvektor (Quadratisch) */
  double *erg;   /* Ergebnis der Berechnung */

  for (i=0;i<DIM;i++)
  {
    mat[i] = (double *)calloc(DIM,sizeof(double));
  }
  erg = (double *)calloc(DIM-1,sizeof(double));



   mat[0][0] = 1.0;
   mat[0][1] = 2.0;
   mat[0][2] = -0.7;
   mat[0][3] = 21.0;

   mat[1][0] = 3.0;
   mat[1][1] = 0.2;
   mat[1][2] = -1.0;
   mat[1][3] = 24.0;

   mat[2][0] = 0.9;
   mat[2][1] = 7.0;
   mat[2][2] = -2.0;
   mat[2][3] = 27.0;

   gauss((double **)mat,erg,3);

   for (i=0;i<3;i++)
   {
     printf("%f\n",erg[i]);

   }
   return(0);

© 1995 by Thomas Dorn
Design by comdes