Untitled

std::map<rational, vectorspace> math::find_eigenspaces(matrix mat, std::vector<rational> eigenvalues)
{
    std::map<rational, vectorspace> eigenspaces{};
    rational coef{};
    if (mat.get_size()[0] != mat.get_size()[1])
        throw std::exception("Ошибка при нахождении собственных подпространств: матрица должна быть квадратной");
    for (int i = 0; i < static_cast<int>(eigenvalues.size()); i++)
    {
        matrix diag_mat = matrix::create_diag_mat(mat.get_size()[0], eigenvalues[i]);
        matrix mat_to_solve = (mat - diag_mat).concat_columns(matrix(mat.get_size()[0], 1));
        matrix gaussian_elim_mat = gaussian_elimination(mat_to_solve, coef).del_empty_rows();

        int m = gaussian_elim_mat.get_size()[0];
        int n = gaussian_elim_mat.get_size()[1];

        std::vector<int> basis_var_indexes = {};
        std::vector<int> free_var_indexes = {};
        for (int j = 0; j < m; j++)
            for (int k = j; k < n - 1; k++)
                if (gaussian_elim_mat[j][k] != 0)
                {
                    basis_var_indexes.push_back(k);
                    break;
                }
        for (int j = 0; j < n - 1; j++)
            if (std::find(basis_var_indexes.begin(), basis_var_indexes.end(), j) == basis_var_indexes.end())
                free_var_indexes.push_back(j);

        std::vector<int> row_indexes_arr(m);
        std::iota(std::begin(row_indexes_arr), std::end(row_indexes_arr), 0);
        matrix sliced_mat = gaussian_elim_mat.intersection_slice(row_indexes_arr, basis_var_indexes);

        std::vector<point> basis_vectors{};
        for (auto& f : free_var_indexes) 
        {
            std::vector<rational> solution(n - 1);
            solution[f] = 1;
            matrix mat_b(m, 1);

            for (int j = 0; j < m; j++)
                mat_b[j][0] = gaussian_elim_mat[j][f] * -1;

            matrix expand_sliced_mat = sliced_mat.concat_columns(mat_b);
            std::vector<rational> gauss_solution = solve_square_matrix_by_gauss(expand_sliced_mat);

            for (int j = 0; j < static_cast<int>(gauss_solution.size()); j++)
                solution[basis_var_indexes[j]] = gauss_solution[j];
            basis_vectors.push_back(solution);
        }

        eigenspaces[eigenvalues[i]] = vectorspace(basis_vectors);
    }

    return eigenspaces;