Gauss Jordan sustitution

#include <iostream>

#define SIZE 3
using namespace std;

void print_matrix(int M[SIZE][SIZE + 1]) {
  // La stampiamo
  for (int i = 0; i < SIZE; i++) {
    for (int j = 0; j < SIZE + 1; j++) {
      cout << M[i][j] << " ";
    }
    cout << endl;
  }
}

double solve_matrix(int M[SIZE][SIZE + 1]) {
  double
      solutions[SIZE]; // l'array delle soluzioni che ho trovato alle incognite

  for (int i = 0; i < SIZE; i++) {
    cout << SIZE - i << " , " << SIZE - i + 1 << " : ";
    cout << M[SIZE - i - 1][SIZE - i] << endl;
    for (int j = 0; j < i; j++) {
      cout << M[(SIZE - i)][SIZE - i + 1] << " * " << solutions[j] << " = "
           << M[SIZE - i][SIZE - i + 1] * solutions[j] << endl;
      M[SIZE - i - 1][SIZE - i - 1] =
          M[SIZE - i - 1][SIZE - i - 1] * solutions[j];
    }
    solutions[i] = (double)M[SIZE - i - 1][SIZE - i] /
                   (double)M[SIZE - i - 1][SIZE - i - 1];
  }

  cout << "Soluzione della Z: " << solutions[0] << endl;
  return 0;
}

int main(void) {
  // Dichiriamo una SIZE x SIZE
  int M[SIZE][SIZE + 1];

  // La carichiamo
  for (int i = 0; i < SIZE; i++) {
    for (int j = 0; j < SIZE + 1; j++) {
      //   cout << "M"
      //        << "[" << i << "]"
      //        << "[" << j << "]: ";
      cin >> M[i][j];
    }
  }

  print_matrix(M);

  cout << endl << endl;
  solve_matrix(M);

  cout << endl;
  print_matrix(M);
  return 0;
}

/**
 *     [0][1][2][3]
 *      #  #  #  S          [0]
 *      0  #  #  S          [1]
 *      0  0  #  S          [2]
 */