#include<stdio.h>
#include<math.h>

#define MAX 100
typedef double matrix[MAX][MAX];

void print_matrix(matrix A, int w, int k)
{
   printf("\nMacierz %dx%x\n", w, k);

   for(int i=0; i< w;i++)
   {
      for(int j=0; j< k; j++)
         printf(" %5.1lf", A[i][j]);
      printf("\n");
   }
}

void rozklad_lu(const matrix A, matrix L, matrix U, int n)
{
   for (int i = 0; i < n; i++)
   {
      for (int j = 0; j < i; j++)
      {
         U[i][j] = 0;
         L[j][i] = 0;
      }

      L[i][i] = 1.0;

      for (int j = i; j < n; j++)
      {
         U[i][j] = A[i][j]; 
         for (int k = 0; k <= i-1; k++)
         {
            U[i][j] -= L[i][k] * U[k][j]; 
         }         
      }

      for (int j = i+1; j < n; j++)
      {
         L[j][i] = A[j][i]; 
         for (int k = 0; k <= i-1; k++)
         {
            L[j][i] -= L[j][k] * U[k][i]; 
         }         
         L[j][i] /= U[i][i];
      }
   }
}

void mul(const matrix A, const matrix B, matrix C, int nwa, int nka, int awb, int nkb)
{
   for (int i = 0; i < nwa; i++)
   {
      for (int j  = 0; j < nkb; j++)
      {
         C[i][j] = 0.0;
         for (int k = 0; k < nka; k++)
         {
            C[i][j] += A[i][k] * B[k][j];
         }
      }
   }
}

int main()
{
   matrix A = {{1,2,3}, {4,5,6}, {7,8,9}};
   
   matrix L = {{0,0,0}, {0,0,0}, {0,0,0}};
   matrix U = {{0,0,0}, {0,0,0}, {0,0,0}};

   matrix C;

   printf("Macierz A\n");
   print_matrix(A, 3, 3);
   
   rozklad_lu(A, L, U, 3);
   
   printf("Macierz L\n");
   print_matrix(L, 3, 3);

   printf("Macierz U\n");
   print_matrix(U, 3, 3);

   mul(L, U, C, 3, 3, 3, 3);
   printf("Sprawdam A = LU\n");
   print_matrix(C, 3, 3);

   return 0;
}