Decomposiçao LU em C/C++

[RSABINO]

Aoo galera do Clube do Hardware!

Tenho um trabalho para a faculdade e nao consigo conclui-lo. Preciso fazer uma programaçao em C/C++ para resolver um sistema linear, mas em DECOMPOSIÇAO LU (uma matriz A = duas matrizes L e U, triangular inferior e superior respectivamente).

Como faço para por em rotina de C/C++ essa decomposiçao??? Alguem já fez isso em algum lugar e tem isso pronto??

Tenho uma rotina semi-pronta. Se alguem puder me dizer porque o programa acusa erro e fecha, eu agradeço.

#include <stdio.h>

#include <conio.h>

#include <math.h>

#define NMAX 50

int main()

{

int i, j, nmais, a[NMAX][NMAX+1], U[NMAX][NMAX], L[NMAX][NMAX], var;

int aux, aux2, aux3, aux4, n, k, coef;

printf("Numero de equacoes: ");

scanf("%d", &n);

printf("\nCoeficientes do sistema de equacoes lineares\n");

for (i = 1; i <= n; i++) { //contagem para linhas

nmais = n +1;

printf("");

for (j = 1; j <= nmais; j++) { //contagem para colunas ( coloca-se uma coluna extra pois é a coluna dos coeficientes indep.)

printf(" A( %d , %d ) = ", i, j);

scanf("%d", &coef);

a[j] = coef;

} //fecha for

} //fecha for

/*i = 1;

j = 1;

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

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

printf("/n A[ %d ][ %d ] = %d", i, j, coef);

}

}*/

i = 1;

j = 1;

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

U[1][j] = a[1][j];

L[1] = a[1]/U[1][1];

i = i+1;

} //fecha for

aux = 0;

aux2 = 0;

i = 2;

j = 2;

var = 0;

while (var <= n-1){

if (i <= n){

for (k = 1; k <= j-1; k++){ //loop para somatoria (matriz L)

aux = L[k]*U[k][j];

aux2 = aux2 + aux;

L[j] = (a[j] - aux2);

L[j] = L[j]/U[j][j];

i = i+1;

} //fecha for

} //fecha while

else{

if (j <= n){

for (k=1; k <= i-1; k++){

aux3 = L[k]*U[k][j];

aux4 = aux4 + aux3;

U[j] = a[j] - aux4;

j = j+1;

} //fecha for

} //fecha if

else

break;

} //fecha else

var = var+1;

} //fecha while

getch();

}