Lineare Gleichungssysteme Lösen

Thorsten90

Hi,

ich versuche, das lösen linearer Gleichungssysteme der Form A*x=b zu programmieren. Dabei ist A eine Matrix, b und x Vektoren. der vektor x soll berechnet werden. Gauß'sches Eliminationsverfahren habe ich bereits erfolgreich programmiert, jetzt versuche ich Gauss-Seidel- und Jakobi-Verfahren zur iterativen Lösung zu schreiben. Ich habe beides als Funktion ausgelagert. Das erste Gleichungssystem, das ich zum testen erstellt habe, wird korrekt berechnet, aber danach hängt sich das Programm auf, als würde es eine Endlosschleife beginnen. Ich hoffe, jemand hier kann mir bei der Suche nach der Lösung helfen.

double* matrix_solve_gauss_seidel(double** matrix, double* vector_result, int lines, double residual)
{
    int i, j;
    double* vector_solution={NULL};
    vector_solution=(double*)calloc(lines, sizeof(double));
    double residual_tmp_1=1;
    double residual_tmp_2=1;
    while (residual_tmp_2>residual)
    {
        for(i=0; i<lines; i++)
        {
            vector_solution[i]=vector_result[i]/matrix[i][i];
            for(j=0; j<lines; j++)
            {
                if(j!=i) vector_solution[i]=vector_solution[i]-matrix[i][j]*vector_solution[j]/matrix[i][i];
            }
        }
        residual_tmp_2=0;
        for(i=0; i<lines; i++)
        {
            residual_tmp_1=vector_result[i];
            for(j=0; j<lines; j++)
            {
                residual_tmp_1=residual_tmp_1-vector_solution[j]*matrix[i][j];
            }
            if(fabs(residual_tmp_1)>residual_tmp_2) residual_tmp_2=fabs(residual_tmp_1);
        }
    }
    return vector_solution;
}
double* matrix_solve_jacobi(double** matrix, double* vector_result, int lines, double residual)
{
    int i, j;
    double* vector_solution={NULL};
    vector_solution=(double*)calloc(lines, sizeof(double));
    double* vector_solution_tmp={NULL};
    vector_solution_tmp=(double*)calloc(lines, sizeof(double));
    double residual_tmp_1=1;
    double residual_tmp_2=1;
    while (residual_tmp_2>residual)
    {
        for(i=0; i<lines; i++)
        {
            vector_solution_tmp[i]=vector_result[i]/matrix[i][i];
            for(j=0; j<lines; j++)
            {
                if(j!=i) vector_solution_tmp[i]=vector_solution_tmp[i]-matrix[i][j]*vector_solution[j]/matrix[i][i];
            }
        }
        for(i=0; i<lines; i++) vector_solution[i]=vector_solution_tmp[i];
        residual_tmp_2=0;
        for(i=0; i<lines; i++)
        {
            residual_tmp_1=vector_result[i];
            for(j=0; j<lines; j++)
            {
                residual_tmp_1=residual_tmp_1-vector_solution[j]*matrix[i][j];
            }
            if(fabs(residual_tmp_1)>residual_tmp_2) residual_tmp_2=fabs(residual_tmp_1);
        }

    }
    free (vector_solution_tmp);
    return vector_solution;
}

Und wie und mit welchen Werten rufst du die Funktionen auf?

Thorsten90

Der Code sieht im Prinziep so aus, Die Funktionen stehen in seperaten header und cpp Files.

int main()
{
    double** M1={NULL};
    M1=matrix_create(M1, 3, 3);
    M1[0][0]=5;
    M1[0][1]=2;
    M1[0][2]=1;
    M1[1][0]=2;
    M1[1][1]=4;
    M1[1][2]=2;
    M1[2][0]=1;
    M1[2][1]=1;
    M1[2][2]=4;
    double* V1={NULL};
    V1=(double*)calloc(3, sizeof(double));
    V1[0]=1;
    V1[1]=6;
    V1[2]=3;
    double* V11={NULL};
    double* V12={NULL};
    double residual=0.000001;

    V11=matrix_solve_gauss_seidel(M1, V1, 3, residual);
    V12=matrix_solve_jacobi(M1, V1, 3, residual);
    V1=matrix_solve_gauss_elimination(M1, V1, 3);

Danach werden die Ergebnisse noch in der konsole ausgegeben. Gauss-Eliminationsverfahren funktioniert, daher habe ich die Funktion nicht mit angegeben. Übergeben werden Matrix und Vektor, die das lineare Gleichungssystem bilden. zurückgegeben wird der Zeiger für den Vektor mit den Resultaten(V11 und V12). Beim Gauss-Eliminationsverfahren geht die matrix verloren und der Vektor(V1) wird zum Ergebnis. Das Programm enthält 5 Gleichungssysteme zum testen der Funktionen, aber seid ich einige Fehler korrigiert habe (vorher gab es für Gauss-Seidel und Jakobi verfahren nur falsche Ergebnisse) löst das Programm nurnoch das erste Gleichungssystem und scheint in einer Endlosschleife hängen zu bleiben, wenn es das zweite lösen soll. Wenn ich die Funktionen getrennt rausnehme hängt es ebenfalls oder bricht sogar ab.

Bitsy

Bau das Programm mal so um, dass Speicherreservierungen und -Freigaben innerhalb einer Routine passieren.
Dein matrix_solve_gauss_seidel z.B. fordert mit calloc an, gibt aber nichts mehr frei, die Außenwelt muß sich also um die Freigabe kümmern - was nicht gut ist.

Interessant wäre auch die Frage, woher jene Funktion wissen soll, wie groß die
zweite Dimension von matrix ist, wenn es auf matrix[...][i] zugreift... das steht nämlich in einem ** nicht mit drin.

Thorsten90

Hi,

die übergebenen Matrizen sind quadratisch, also muss nur eine Zahl übergeben werden für Zeilen und Spalten.

Den angeforderten Speicher freizugeben macht keinen Sinn, die Funktion erstellt ein Array mit den Lösungen des Gleichungssystems und giebt dieses zurück. Gebe ich den Speicher in der Funktion frei lösche ich damit meine Ergebnisse. Außerdem sollen sie außerhalb der Funktion zur verfügung stehen.

Dann übergib den Funktionen genug Platz für das Ergebnis. Die rufende Funktion soll sich um den Speicher kümmern.

double* matrix_solve_gauss_seidel(double* vector_solution, double** matrix, double* vector_result, int lines, double residual);

Du kannst auch weiterhin den Zeiger von vector_solution zurückgeben.

Gibt es das Problem auch schon in dem kurzen Code von der main, die du gepostet hast?
Du sprichst von Wiederholungen. Die hat deine main aber gar nicht.

Poste ein Minimalbeispiel, bei dem der Fehler noch auftritt.

Thorsten90

Habe die Funktion so geändert, das vector_solution als Parameter mit übergeben wird.

Das ist mein vollständiges Programm

int main()
{
    double** M1={NULL};
    M1=matrix_create(M1, 3, 3);
    M1[0][0]=5;    M1[0][1]=2;    M1[0][2]=1;    M1[1][0]=2;    M1[1][1]=4;
    M1[1][2]=2;    M1[2][0]=1;    M1[2][1]=1;    M1[2][2]=4;    
    double** M2={NULL};
    M2=matrix_create(M2, 2, 2);
    M2[0][0]=0.01;    M2[0][1]=1;    M2[1][0]=1;    M2[1][1]=1;
    double** M3={NULL};
    M3=matrix_create(M3, 5, 5);
    M3[0][0]=5;    M3[0][1]=4;    M3[0][2]=3;    M3[0][3]=2;    M3[0][4]=1;
    M3[1][0]=3;    M3[1][1]=2;    M3[1][2]=1;    M3[1][3]=0;    M3[1][4]=4;
    M3[2][0]=1;    M3[2][1]=0;    M3[2][2]=1;    M3[2][3]=0;    M3[2][4]=1;
    M3[3][0]=0;    M3[3][1]=1;    M3[3][2]=1;    M3[3][3]=1;    M3[3][4]=4;
    M3[4][0]=2;    M3[4][1]=1;    M3[4][2]=0;    M3[4][3]=2;    M3[4][4]=1;
    double** M4={NULL};
    M4=matrix_create(M4, 6, 6);
    M4[0][0]=1;    M4[0][1]=1;    M4[0][2]=1;    M4[0][3]=1;    M4[0][4]=1;
    M4[0][5]=1;    M4[1][0]=1;    M4[1][1]=2;    M4[1][2]=1;    M4[1][3]=0;
    M4[1][4]=0;    M4[1][5]=0;    M4[2][0]=0;    M4[2][1]=1;    M4[2][2]=2;
    M4[2][3]=0;    M4[2][4]=1;    M4[2][5]=1;    M4[3][0]=1;    M4[3][1]=1;
    M4[3][2]=0;    M4[3][3]=0;    M4[3][4]=0;    M4[3][5]=1;    M4[4][0]=0;
    M4[4][1]=2;    M4[4][2]=0;    M4[4][3]=2;    M4[4][4]=0;    M4[4][5]=2;
    M4[5][0]=1;    M4[5][1]=2;    M4[5][2]=3;    M4[5][3]=4;    M4[5][4]=5;
    M4[5][5]=6;
    double** M5={NULL};
    M5=matrix_create(M5, 3, 3);
    M5[0][0]=4;    M5[0][1]=1;    M5[0][2]=1;    M5[1][0]=2;    M5[1][1]=8;
    M5[1][2]=2;    M5[2][0]=1;    M5[2][1]=1;    M5[2][2]=4;
    double* V1={NULL};
    V1=(double*)calloc(3, sizeof(double));
    V1[0]=1;    V1[1]=6;    V1[2]=3;
    double* V2={NULL};
    V2=(double*)calloc(2, sizeof(double));
    V2[0]=0.5;    V2[1]=1;
    double* V3={NULL};
    V3=(double*)calloc(5, sizeof(double));
    V3[0]=15;    V3[1]=10;    V3[2]=3;    V3[3]=7;    V3[4]=6;
    double* V4={NULL};
    V4=(double*)calloc(6, sizeof(double));
    V4[0]=9;    V4[1]=6;    V4[2]=7;    V4[3]=5;    V4[4]=12;    V4[5]=33;
    double* V5={NULL};
    V5=(double*)calloc(3, sizeof(double));
    V5[0]=9;    V5[1]=24;    V5[2]=15;

    double* V11={NULL};
    double* V12={NULL};
    double* V21={NULL};
    double* V22={NULL};
    double* V31={NULL};
    double* V32={NULL};
    double* V41={NULL};
    double* V42={NULL};
    double* V51={NULL};
    double* V52={NULL};
    double residual=0.000001;

    V11=matrix_solve_gauss_seidel(M1, V1, V11, 3, residual);
    V12=matrix_solve_jacobi(M1, V1,V12, 3, residual);
    V1=matrix_solve_gauss_elimination(M1, V1, 3);
    matrix_output_display_solution(V1, 3);
    matrix_output_display_solution(V11, 3);
    matrix_output_display_solution(V12, 3);
    free(V11);
    free(V12);
    free(V1);
    matrix_free(M1, 3);
    printf("\n\n");
    V21=matrix_solve_gauss_seidel(M2, V2,V21, 2, residual);
    V22=matrix_solve_jacobi(M2, V2,V22, 2, residual);
    V2=matrix_solve_gauss_elimination(M2, V2, 2);
    matrix_output_display_solution(V2, 2);
    matrix_output_display_solution(V21, 2);
    matrix_output_display_solution(V22, 2);
    free(V21);
    free(V22);
    free(V2);
    matrix_free(M2, 2);
    printf("\n\n");
    V31=matrix_solve_gauss_seidel(M3, V3,V31, 5, residual);
    V32=matrix_solve_jacobi(M3, V3,V32, 5, residual);
    V3=matrix_solve_gauss_elimination(M3, V3, 5);
    matrix_output_display_solution(V3, 5);
    matrix_output_display_solution(V31, 5);
    matrix_output_display_solution(V32, 5);
    free(V31);
    free(V32);
    free(V3);
    matrix_free(M3, 5);
    printf("\n\n");
    V41=matrix_solve_gauss_seidel(M4, V4,V41, 6, residual);
    V42=matrix_solve_jacobi(M4, V4,V42, 6, residual);
    V4=matrix_solve_gauss_elimination(M4, V4, 6);
    matrix_output_display_solution(V4, 6);
    matrix_output_display_solution(V41, 6);
    matrix_output_display_solution(V42, 6);
    free(V41);
    free(V42);
    free(V4);
    matrix_free(M4, 6);
    printf("\n\n");
    V51=matrix_solve_gauss_seidel(M5, V5,V51, 3, residual);
    V52=matrix_solve_jacobi(M5, V5,V52, 3, residual);
    V5=matrix_solve_gauss_elimination(M5, V5, 3);
    matrix_output_display_solution(V5, 3);
    matrix_output_display_solution(V51, 3);
    matrix_output_display_solution(V52, 3);
    free(V51);
    free(V52);
    free(V5);
    matrix_free(M5, 3);
    printf("\n\n");
}

Am Anfang stehen fünf Matrizen und die Vektoren für die Lösungen. Das erste Gleichungssystem wird fehlerfrei gelöst, danach hängt sich das Programm auf als würde eine Endlosschleife gestartet werden. Wenn ich das zweite Gleichungssystem auskommentiere läuft das Programm zwar durch, aber nur für das erste und letzte Gleichungssystem gibt es korrekte Ergebnisse. Für die anderen beiden sind die Ergebnisse völlig sinnlos.

Bitsy

Huhu - hast meine Antwort nicht verstanden!

** und [][] haben nicht unbedingt was miteinander zu tun.

Benutze ich [i][j], dann liegt mein Element an Position i*YDIM+j im Speicher.
- und YDIM kennt die Funktion nicht!

matrix_create wäre mal interessant - dann läßt es sich besser erklären.

Bitsy schrieb:

** und [][] haben nicht unbedingt was miteinander zu tun.

Doch, haben sie. z.B. bei einem *double[] (einem Array aus Zeigern auf double)

Bitsy schrieb:

Benutze ich [i][j], dann liegt mein Element an Position i*YDIM+j im Speicher.

Das gilt für echte 2D-Arrays.

Arrays sind während ihrer Lebenszeit nicht mehr in der Größe und im Speicherposition änderbar.

Bitsy

Eben, und deshalb ist matrix_create interessant.

Thorsten90

double** matrix_create (double** matrix, int lines, int columns)
{
    int i;
    matrix=(double**)calloc(lines, sizeof(*matrix));
    for (i=0; i<lines; i++) matrix[i]=(double*)calloc(columns, sizeof(**matrix));
    return (matrix);
}

So sieht matrix_create aus.

Bitsy

Jo, das ist dann okay.
Ist dann wohl kein Speicherproblem, sondern was im mathematischen Verfahren.

Ich hab's mal nachgestellt und er bleibt bei mir dann ebenfalls hängen.
Da stimmt noch was nicht mit den residual-Werten, die steigen auf unendlich, und dann kann er die Abbruchbedingung nicht erfüllen.
(Im Jacobi bleibt er hängen, im Gauss-Seidel liefert er INF/-INF für x und y, also auch falsch. Und ja, die 1.Matrix macht er korrekt)

Wutz

Bitsy schrieb:

Jo, das ist dann okay.
Ist dann wohl kein Speicherproblem, sondern was im mathematischen Verfahren.

Bitsy schrieb:

Huhu - hast meine Antwort nicht verstanden!

** und [][] haben nicht unbedingt was miteinander zu tun.

?? Du weißt aber auch nicht so recht, was du erzählst?
Die Dereferenzierung von double**matrix mittels matrix[i][i] usw. wird knallen, wie nur sonst irgendwas, und zwar aus dem von dir selbst zuvor genanntem Grund.

Bitsy

Tut sie nicht, und zwar aus dem von DirkB genannten Grund.
Deswegen wollte ich ja matrix_create sehen.

Wenn Du die Sourcen jetzt nimmst, wie sie derzeit da stehen, mußt Du noch vector_solution als Parameter zu den Matrix-Lösefunktionen hinzufügen und V11...V52 mit calloc vorbelegen, aber dann liegen die Daten stabil wo sie hingehören und er kann bequem mit [][] darauf zugreifen. (Eben weil es über ein Pointer-Array läuft - ist schon etwas verwirrend). Das Problem liegt im mathematischen Bereich. Leider kenne ich diese Näherungsverfahren nicht.

Da, guck, hab's extra nichtquadratisch gemacht:

double** matrix_create(int lines, int columns)
{
  int i;
  double** m = (double**)calloc(lines, sizeof(double*));
  for (i=0; i<lines; i++)
    m[i] = (double*)calloc(columns, sizeof(double));
  return m;
}

void matrix_out(double** m, int d1, int d2)
{
  for (int i=0; i<d1; i++)
  {
    for (int j=0; j<d2; j++)
      printf("%f ", m[i][j]);
    printf("\n");
  }
  // aber der rauscht jetzt ab, weil die zweite Dimension
  // überschritten ist
  printf("%f\n", m[1][3]);
}

int main()
{
  double** M=matrix_create(3, 2);
  M[0][0]=5; M[0][1]=2;
  M[1][0]=2; M[1][1]=4;
  M[2][0]=1; M[2][1]=1;
  matrix_out(M, 3, 2);
  // free...
  return 0;
}

Das hier wäre fatal gewesen, ohne dass er den cast angemeckert hätte:

double** matrix_create(int lines, int columns)
{
  return (double**)calloc(lines*columns, sizeof(double));
}

Thorsten90

Also was das Rechenverfahren angeht (falls es hilft):

Gelöst werden soll A*x=b, A ist eine Matrix, x und b sind Vektoren. Wenn man das ganze als einzelne gleichungen schreibt sind es ja im Grunde alles Gleichungen, z.b. sowas wie 2x1+3x2-x3+2x4=5. Sagen wir mal das wäre die erste Zeile. umgestellt ist dann x1=5/2-3/2+2+1/2x3-x4. So kann jede Zeile umgestellt werden. Es ist also x(i)=1/a(ii)*( b(i)- summe (j=1 bis n, j!=i) x(j)*a(ij) ). dabei ist x in meinem programm vector_solution[i], A ist matrix[i][j] und b ist vector_result[j]. Mit jedem Iterationsschritt wird dabei jedes x genauer und nährt sich so immer weiter der Lösung. Der Unterschied zwischen gauß-Seidel und Jakobi ist, das bei Gauß-Seidel immer das x aus vector_solution benutzt wird, während bei Jakobi die Ergebnisse x erstmal gespeichert werden (vector_solution_tmp) und anschließend erst dem eigentlichen x (vector_solution) übergeben werden bevor eine weitere iteration startet.

Was den Fehler (residual) angeht, der ist eigentlich die selbe Gleichung, nur das x(i) auch auf die andere seite gebracht wurde, zu 0=b-A*x. Das Programm sucht den größten Einzelfehler aus allen Zeilen, der wird dann mit dem übergebenen Maximalwert verglichen.

Ich hoffe, ich habe es einigermaßen verständlich erklährt, vlt hilft es ja bei der Suche nach dem Fehler.

Wo ich aber eure Beiträge so lese hätte ich noch ein paar Frage:

Die Matrizen und Vektoren in meinem Programm sind nur zum testen, es sollen hinterher Gleichungssysteme eingelesen werden. Aber sie machen das Programm irgendwie unübersichtlich. Gibt es für die dynamischen arrays auch eine kürzere möglichkeit, sie zu initialisieren, so wie bei statischen mit a[n]={...}?

Kann ich einer solchen Funktion auch ein statisches Array übergeben? Hatte es mal mit (double**) array als Übergabeparameter versucht, aber das hat nicht funktioniert.

Ihr diskutiert so viel über meine Variante, gibt es vlt eine bessere möglichkeit, dynamische Arrays an Funktionen zu übergeben bzw. als Rückgabewert zu setzen?

Bitsy

Ich versuche mich mal reinzulesen. Das kriegen wir schon gebacken.

Zu den Fragen nach den Arrays:
Create und Free hast Du ja jetzt schon, und viel komfortabler wirst Du es in C nicht hinbekommen. Statisch würde ich vermeiden, Du siehst ja, Deine Beispiele haben schon unterschiedliche Dimensionen. Es wäre sinnvoller sich zumindest eine kleine Fileread-Routine zu basteln, dann kloppst Du die Daten einfach mit einem Texteditor in eine Datei, machst einen Doppelklick auf's Programm und fertig.

Was Du jetzt mit diesem Problem aber auch schon mal machen könntest:
Schnupper mal in den Debugger Deines Compilers rein.
Den wirst Du früher oder später brauchen - je eher, desto besser.
Da kannst Du Zwischenwerte auf festgelegten Breakpoints anschauen und kommst der Problemlösung um einiges schneller auf die Spur.

Wutz

Thorsten90 schrieb:

Die Matrizen und Vektoren in meinem Programm sind nur zum testen, es sollen hinterher Gleichungssysteme eingelesen werden. Aber sie machen das Programm irgendwie unübersichtlich.

Richtig erkannt.

Thorsten90 schrieb:

Gibt es für die dynamischen arrays auch eine kürzere möglichkeit, sie zu initialisieren,

Es gibt keine dynamischen Arrays. Es gibt nur Arrays. Und die haben immer eine feste Größe.

Thorsten90 schrieb:

so wie bei statischen mit a[n]={...}?

Wenn du damit VLA meinst, deren Größe nicht zur Compilezeit bekannt ist: nein, VLA lassen sich nicht initialisieren.

Thorsten90 schrieb:

Kann ich einer solchen Funktion auch ein statisches Array übergeben?

Ja, natürlich.

Thorsten90 schrieb:

Hatte es mal mit (double**) array als Übergabeparameter versucht, aber das hat nicht funktioniert.

Da musst du unterscheiden zwischen Matrix (zweidimensionale lückenlose Folge von Elementen gleichen Typs) und Array (eindimensionale lückenlose Folge von Elementen gleichen Typs).
Erstere lässt sich nicht via ** übertragen, letzte mittels .
Ein Array von Zeigern auf double lässt sich via double* oder double*[] übergeben, ist aber dann keine Matrix mehr, sondern wie gesagt ein Array.

Thorsten90 schrieb:

Ihr diskutiert so viel über meine Variante, gibt es vlt eine bessere möglichkeit, dynamische Arrays an Funktionen zu übergeben bzw. als Rückgabewert zu setzen?

Definiere dir ein Array mit fester, maximal vorkommender Größe und verzichte auf malloc/calloc/free, und übergebe die (max.) untere Dimensionsgröße im Parametertyp beim Prototyp:

enum {MAXDIM=20};
void fkt(double m[][MAXDIM],double v[],double bla,int n)
{
  int i;
  assert(n<MAXDIM);
  for(i=0;i<n;++i)
    m[i][i]=...
}

Bitsy

Also, ich habe das jetzt nochmal richtig umgebaut und erhalte den gleichen Effekt.
Wenn's geht, schreib doch mal für die zweite Matrix
0.01 1
1 1
mit dem Zielvektor (0.5 | 1) die jeweils ersten beiden Iterationsschritte manuell auf. Das kann nur ein Start-Problem sein.

Meine Routinen sehen so aus:

double matrix_sum(double** matrix, int index, double lines, double* vin)
{
  double vsum = 0;
  for (int j=0; j<lines; j++)
    if (index != j)    
      vsum += vin[j] * matrix[index][j];
  return vsum;
}

void matrix_solve_gauss_seidel(double** matrix, double* vector_result, double* vector_solution, int lines, double residual)
{
  int i;
  double residual_tmp_1;
  double residual_tmp_2 = 1;

  double* vector_old=(double*)malloc(lines * sizeof(double));
  for (int i=0; i<lines; i++)
    vector_old[i] = 0;

  while (residual_tmp_2>residual)
  {
    residual_tmp_2 = 0;
    for(i=0; i<lines; i++)
    {
      vector_solution[i] =
         (vector_result[i] - matrix_sum(matrix, i, lines, vector_old)) / matrix[i][i];
      residual_tmp_1 = vector_solution[i] - vector_old[i];
      if(fabs(residual_tmp_1)>residual_tmp_2)
        residual_tmp_2=fabs(residual_tmp_1);
    }
    for(i=0; i<lines; i++)
      vector_old[i] = vector_solution[i];
  }
  free(vector_old);
}

Beachte: vector_result bleibt hier stabil der Zielvektor und vector_old wird auf 0|0 gestartet. Endet die Routine, liegt das Resultat in vector_solution.
Irgendwo mu0 hier der Denkfehler sein.
Damit iteriert er sich bei der ersten 3*3-Matrix wunderschön ran, aber bei der
kleinen haut er ab!
In anderen Quellen habe ich was davon gesehen, dass die Matrix zu Beginn irgendwie transformiert werden muss - ist da was dran?

Bitsy

Also, ich habe das jetzt nochmal richtig umgebaut und erhalte den gleichen Effekt.

Meine Routinen sehen so aus: (das scheint nun eher Jacobi zu sein als Gauss-Seidel)

double matrix_sum(double** matrix, int index, double lines, double* vin)
{
  double vsum = 0;
  for (int j=0; j<lines; j++)
    if (index != j)    
      vsum += vin[j] * matrix[index][j];
  return vsum;
}

void matrix_solve_gauss_seidel(double** matrix, double* vector_result, double* vector_solution, int lines, double residual)
{
  int i;
  double residual_tmp_1;
  double residual_tmp_2 = 1;

  double* vector_old=(double*)malloc(lines * sizeof(double));
  for (int i=0; i<lines; i++)
    vector_old[i] = 0;

  while (residual_tmp_2>residual)
  {
    residual_tmp_2 = 0;
    for(i=0; i<lines; i++)
    {
      vector_solution[i] =
         (vector_result[i] - matrix_sum(matrix, i, lines, vector_old)) / matrix[i][i];
      residual_tmp_1 = vector_solution[i] - vector_old[i];
      if(fabs(residual_tmp_1)>residual_tmp_2)
        residual_tmp_2=fabs(residual_tmp_1);
    }
    for(i=0; i<lines; i++)
      vector_old[i] = vector_solution[i];
  }
  free(vector_old);
}

Beachte: vector_result bleibt hier stabil der Zielvektor und vector_old wird auf 0|0 gestartet. Endet die Routine, liegt das Resultat in vector_solution.
Damit iteriert er sich bei der ersten 3*3-Matrix wunderschön ran, aber bei der
kleinen haut er ab! Und dann lese ich das:

http://www2.hs-esslingen.de/~koch/numerik/numerik_skript.pdf,
pdf Seite 14, programmiere diesen Mini und verfolge die Entwicklung

int main(void)
{
  double xalt[2] = { 0.0, 0.0 };
  double xneu[2];

  for (int i=0; i<20; i++)
  {
    //        ziel - sum matrix*alt ohne diag  diag
    xneu[0] = (0.5 -(   0       + 1*xalt[1]))/ 0.01;
    xneu[1] = (1   -( 1*xalt[0] +   0      ))/ 1;

    xalt[0]=xneu[0];
    xalt[1]=xneu[1];
    printf("step %d   %f %f\n", i, xneu[0], xneu[1]);
  }
  return 0;
}

Aha! Es gibt da eine schöne Nebenbedingung für die Lösbarkeit:
siehe hier http://de.wikipedia.org/wiki/Diagonaldominante_Matrix
und - genau das ist diese Matrix nicht! Da ist der Fehler