Merge pull request #221 from gbalduzz/extra_symmetrization

      for (int b0 = 0; b0 < BDmn::dmn_size(); ++b0) {

        for (int b1 = 0; b1 < BDmn::dmn_size(); ++b1) {

          scalartype tmp_0 = f(b0, b1, c_ind, w_ind);

          scalartype tmp_1 = f(b1, b0, new_c_idx, w_0 - w_ind);

          constexpr bool real_hamiltonian = true;  // TODO : replace with actual check.

          if constexpr (real_hamiltonian) {

            const auto tmp1 = f(b0, b1, c_ind, w_ind);

            const auto tmp2 = f(b1, b0, new_c_idx, w_0 - w_ind);  // F(w) = conj(F^t(-w))

            const auto tmp3 = f(b0, b1, c_ind, w_0 - w_ind);      // F(w) = conj(F(-w))

            const auto tmp4 = f(b1, b0, new_c_idx, w_ind);        // F(w) = F^t(w)

          scalartype tmp = (tmp_0 + std::conj(tmp_1)) / 2.;

            const auto tmp = (tmp1 + std::conj(tmp2) + std::conj(tmp3) + tmp4) / 4.;

            f_new(b0, b1, c_ind, w_ind) = tmp;

            f_new(b1, b0, new_c_idx, w_0 - w_ind) = std::conj(tmp);

            f_new(b0, b1, c_ind, w_0 - w_ind) = std::conj(tmp);

            f_new(b1, b0, new_c_idx, w_ind) = tmp;

          }

          else {  // Hamiltonian is complex.

            const auto tmp1 = f(b0, b1, c_ind, w_ind);

            const auto tmp2 = f(b1, b0, new_c_idx, w_0 - w_ind);  // F(w) = conj(F^t(-w))

            const auto tmp = (tmp1 + std::conj(tmp2)) / 2.;

            f_new(b0, b1, c_ind, w_ind) = tmp;

            f_new(b1, b0, new_c_idx, w_0 - w_ind) = std::conj(tmp);

          }

        }

      }

    }

  }

  if (do_diff) {

    double max = 0;

  for (int ind = 0; ind < f.size(); ++ind) {

    for (std::size_t ind = 0; ind < f.size(); ++ind) {

      max = std::max(max, abs(f(ind) - f_new(ind)));

    }

    difference(max, f.get_name(), "WDmn-domain of the function : " + f.get_name() + "\n");

  }

  f = std::move(f_new);

  if (do_diff)

    difference(max, f.get_name(), "WDmn-domain of the function : " + f.get_name() + "\n");

}

template <typename scalartype>