printf("\nStart LU-decomposition\n");
for(k=0;k<N;k++){
if(a[k][k]==0){
printf("a[k][k]==0:\n");
return 1;
}
for(i=k+1;i<N;i++){
fct=a[i][k]/a[k][k];
for(j=k+1;j<N;j++){
a[i][j]-=fct*a[k][j];
}
}
}
printf("\ny[] vector:\n");
for(i=0;i<N;i++){
sum=0.0;
for(j=0;j<i;j++){
sum+=a[i][j]*y[j];
}
y[i]=b[i]-sum;
printf("y[%d]=%9f\n",i,y[i]);
}
printf("\nx[] vector:\n");
for(i=N-1;i>=0;i--){
sum=0.0;
for(j=i+1;j<N;j++){
sum+=a[i][j]*x[j];
}
x[i]=(y[i]-sum)/a[i][j];
printf(" x[%d]=%9f\n",N);
}