Algoritmo di Gauss senza riordinamento pivotale

Devo risolvere un sistema lineare con il metodo di Gauss senza il riordinamento, ma l'output della matrice a scala per righe dà i valori sbagliati nella diagonale principale e io non riesco a capire cosa sto sbagliando.
Questo è il codice completo
C++:

/*Programma che stampa una matrice data
e la riduce a scala per righe con l'algoritmo
di Gauss senza riordinamento pivotale*/
#include <stdio.h>
#include <math.h>
#define N 5
void matrice (double [N][N]);
void sostituzione_indietro (double [N][N], int, double[], double[]);
int main ()

{double a[N][N];
matrice (a);
double b[N]={2,1,1,1,2};
printf ("c="); for (int i=0; i<N; i++) printf ("\t%.0lf\n", b[i]);
double m;
int k, i, j;
for (int k=0; k<N;k++)
    {for (int i=k+1; i<N; i++)
            {m=a[i][k]/a[k][k]; a[i][k]=0;
                for(j=k+1; j<N; j++) { a[i][j]=a[i][j]-m*a[k][j];}
                b[i]=b[i]-m*b[i];}
}

printf ("\nMatrice a scala="); for (int i=0; i<N; i++) {printf ("\n");
for (int j=0; j<N; j++) printf ("\t%.2lf", a[i][j]);}
return 0;}


void matrice (double a[N][N])
{int i,j;
for (i=0; i<N; i++) {for (j=0; j<N; j++) {int s=i-j;
    switch (s){
        case 0: a[i][j]=3;
        break;
        case 1:
        case -1: a[i][j]=-1;
        break;
        default: a[i][j]=0;}
}
}
printf ("B="); for (int k=0; k<N; k++) {for (int h=0; h<N; h++) {printf ("\t%.0lf\t", a[k][h]);} printf ("\n");}
return;}

void sostituzione_indietro (double L[N][N], int n, double b[], double x[])
{x[n]=b[n]/L[n][n];
for (int i=n-1; i>1; i--) {x[i]=b[i]; for (int j=i; j>0; j++) {x[i]=x[i]-L[i][j]*x[j];} x[i]=x[i]/L[i][i];}
for (int k=0; k<n; k++) {printf ("\nx%d=%.2lf\n", k+1, x[k]); }