Ajuda no código de decomposição matricial

Boa noite,

Estou trabalhando num código (elementar) que calcula a inversa de uma matriz usando a chamada decomposição LU.
Não sou especialista em C, mas a saída gerada está errada, e eu não consigo encontrar o erro.
Agradeço de antemão qualquer ajuda

Fiquem com Deus

#include <stdio.h>

#include <stdlib.h>

#include <math.h>



int main(void)

{

    int i,j,k,l,r,s; // contadores

    int n; // dimensão das matrizes que serão alocadas dinâmicamente



    // Entrada da dimensão da matriz

    printf("Por favor, entre com a dimensão da matriz : &quot;

    scanf("%i",&n);

    



    

    // Alocando vetores auxiliares

   float *b = malloc( n * sizeof( float ) );

   float *x = malloc( n * sizeof( float ) );

   float *y = malloc( n * sizeof( float ) );



    // Alocando as matrizes

    

    float **M  = malloc( n * sizeof( float * ) );

    float **I  = malloc( n * sizeof( float * ) );

    

    for( i = 0; i < n; ++i )

    {

        M[i] = malloc( n * sizeof( float ) );

        I[i] = malloc( n * sizeof( float ) );

        

    }    

    // Carregando a matriz via teclado

    for( i = 0; i < n; ++i )

        for( j = 0; j < n; j++ )

        {

            printf("Entrada M%i%i da matriz ? ",i + 1,j + 1);

            scanf("%f",&M[i][j]);

            printf("\n&quot;

        }



    

    // Loop que gera a decomposição LU

    

    for (k = 0; k < n - 1; k++)

    {

        for ( i = k + 1; i < n;i++)

        {

            M[i][k] = M[i][k]/M[k][k];

            for (j = k + 1; j < n;j++)

                M[i][j] = M[i][j] - M[i][k]*M[k][j];

        }

    

    }

    

    

    

    // Loop para geração da inversa usando a decomposição LU

    

    for (i = 0; i < n; i++)

    {

        

        // Carga do i-ésimo vetor da base canônica

        for (j = 0; j < n; j++)

            if (j == i)

                b[j] = (float) 1;

            else

                b[j] = (float) 0;

        

        

        // Cálculo do vetor y 

        

        y[0] = b[0];

        

        for (k = 1;k < n;k++)

        {

            y[k] = b[k];

                

            for (l = 0;l < k;l++)

                y[k] = y[k] - M[k][l]*y[l];    

        }    

            

        // Gerando a matriz inversa    

        I[n - 1][i] = y[n - 1]/M[n - 1][n - 1];

        x[n - 1] = y[n - 1]/M[n - 1][n - 1];

        for (r = n - 2; r >= 0; r--)

        {

            x[r] = y[r];

            for(s = n - 1; s > r;s--)

                x[r] = x[r] - M[r][s]*x[s];

                

            I[r][i] = x[r]/M[r][r];

            

        }

    }

    

    

    // Saída do processamento - Matriz de saída

    

    printf("Decomposição LU \n\n&quot;

    for( i = 0; i < n; i++ )

    {

        for( j = 0; j < n; j++ )

            printf("%f\t",M[i][j]);

            

        printf("\n&quot;

        

    }

    printf("\n\n&quot; 

    

    printf("A Matriz Inversa\n\n&quot;

    for( i = 0; i < n; i++ )

    {

        for( j = 0; j < n; j++ )

            printf("%f\t",I[i][j]);

            

        printf("\n&quot;

        

    }

    printf("\n\n&quot;

    

    

    // Liberando a área de memória alocada para as matrizes 

    for( i = 0; i < n; i++ )

    {

        free(M[i]);    

        free(I[i]);

    }

    free(M);

    free(I);

    

    // Liberando a área de memória alocada para o vetores 

    free(b);

    free(x);

    free(y);

        

        

    return(0);

}