void LU(float A[MAX][MAX], float U[MAX][MAX], float L[MAX][MAX], const int tam) {
  float sumatoria;
  for (int i = 0; i < tam; i++)
    L[i][i] = 1;

  for (int k = 0; k < tam; k++) {
    for (int j = 0; j < tam; j++) {
      sumatoria = 0;
      for (int s = 0; s < tam; s++) {
        sumatoria += L[k][s] * U[s][j];
        if (k==tam-1 && j==1) {
          cout << "ERROR  " << sumatoria << endl;
          cout << "L[" << k << "," << s << "] X U[" << s << "," << j << "] = " << L[k][s]*U[s][j] << endl;
        }
      }
      U[k][j] = (A[k][j] - sumatoria) / L[k][k];
    }
    for (int i = k+1; i < tam; i++) {
      sumatoria = 0;
      for (int s = 0; s < tam; s++) {
        sumatoria += L[i][s] * U[s][k];
      }
      L[i][k] = (A[i][k] - sumatoria) / U[k][k];
    }
  }
}