Decomposição LU

[Arthur Santos]

Então, pessoal, tentei fazer o código da Fatoração LU, um método numérico utilizado para resolver  sistemas de equações lineares. O princípio do método está em dividir a matriz dos coeficientes em duas, obtendo uma matriz triangular superior e outra inferior. Até aí tudo bem, acontece que meu programa volta valores absurdos. Abaixo o código:
 

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

 

int main()

{

    int i,j,k,r;

    float a[5][5];

    float b[5];

    float c[5];

    float m,aux,soma;

    float y[5];

    float x[5];

    int p[5];

    float pv;

 

/* Receber valores da matriz*/

    for(i=0;i<5;i++){

      for(j=0;j<5;j++){

        a[j]=0;

       }

    }

    a[0][0]=3;

    a[0][2]=-1;

    a[1][0]=-1;

    a[1][1]=1;

    a[2][1]=-2;

    a[2][2]=9;

    a[3][1]=1;

    a[3][2]=6;

    a[3][3]=-9;

    a[3][4]=2;

    a[4][0]=5;

    a[4][1]=1;

    a[4][4]=-6;

 

    b[0]=40;

    b[1]=0;

    b[2]=350;

    b[3]=0;

    b[4]=0;

    /* COMEÇO DO ALGORITMO */

 

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

    {

        p=i;

    }

 

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

    {

        pv=abs(a[k][k]);

        r=k;

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

        {

            if(abs(a[k])>pv)

            {

                pv=abs(a[k]);

                r=i;

            }

        }

        if(pv==0)

        {   printf("A matriz é singular");

            return 0;

        }

        if(r=!k)

        {

            aux=p[k];

            p[k]=p[r];

            p[r]=aux;

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

            {

                aux=a[k][j];

                a[k][j]=a[r][j];

                a[r][j]=aux;

            }

        }

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

        {

            m=(a[k])/a[k][k];

            a[k]=m;

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

            {

                a[j]=a[j]-m*(a[k][j]);

            }

        }

    }

 

    /* RESOLUÇÃO DOS SISTEMAS TRIANGULARES */

    for(i=0;i<5;i++)    // c= Pb //

    {

        r=p;

        c=b[r];

    }

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

    {

        soma=0;

        for(j=0;j<(i-1);j++)

        {

            soma=soma+(a[j])*(y[j]);

            y=c-soma;

        }

    }

    for(i=4;i>=0;i--)

    {

        soma=0;

        for(j=i+1;j<5;j++)

        {

            soma=soma+(a[j])*x[j];

            x=(y-soma)/a;

        }

    }

 

    /* IMPRIME OS VALORES X */

 

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

    {

        printf("\n %5.2f \n",x);

    }

return 0;

 

}