C++ gf_view类代码示例

 void inverse_fourier_impl  (gf_view<imtime,matrix_valued> gt , gf_view<imfreq,matrix_valued> const gw, matrix_valued){
  impl_worker w;
  for (size_t n1=0; n1<gw.data().shape()[1];n1++)
   for (size_t n2=0; n2<gw.data().shape()[2];n2++){
    auto gt_sl=slice_target_to_scalar(gt, n1, n2);
    auto gw_sl=slice_target_to_scalar(gw, n1, n2);
    w.inverse( gt_sl, gw_sl);
   }
 }

 void fit_tail(gf_view<imfreq> gf, tail_view known_moments, int n_moments, int n_min, int n_max,
   bool replace_by_fit ) {
  if (get_target_shape(gf) != known_moments.shape()) TRIQS_RUNTIME_ERROR << "shape of tail does not match shape of gf";
  gf.singularity() = fit_tail_impl(gf, known_moments, n_moments, n_min, n_max);
  if (replace_by_fit) { // replace data in the fitting range by the values from the fitted tail
   int i = 0;
   for (auto iw : gf.mesh()) { // (arrays::range(n_min,n_max+1)) {
    if (i >= n_min) gf[iw] = evaluate(gf.singularity(),iw);
    i++;
   }
   }
  }

 // compute a tail from the Legendre GF
 // this is Eq. 8 of our paper
 local::tail_view get_tail(gf_view<legendre> const & gl, int size = 10, int omin = -1) {

   auto sh = gl.data().shape().front_pop();
   local::tail t(sh, size, omin);
   t.data() = 0.0;

   for (int p=1; p<=t.order_max(); p++)
     for (auto l : gl.mesh())
       t(p) += (triqs::utility::legendre_t(l.index(),p)/pow(gl.domain().beta,p)) * gl[l];

   return t;

 }

 void legendre_matsubara_direct(gf_view<imfreq> gw, gf_const_view<legendre> gl) {

  gw() = 0.0;
  triqs::arrays::range R;

  // Use the transformation matrix
  for (auto om : gw.mesh()) {
   for (auto l : gl.mesh()) {
    gw[om] += legendre_T(om.index(), l.index()) * gl[l];
   }
  }

  gw.singularity() = get_tail(gl, gw.singularity().size(), gw.singularity().order_min());
 }

 tqa::matrix<double> density( gf_view<legendre> const & gl) { 

   auto sh = gl.data().shape().front_pop();
   tqa::matrix<double> res(sh);
   res() = 0.0;

   for (auto l : gl.mesh()) {
     res -= sqrt(2*l.index()+1) * gl[l];
   }
   res /= gl.domain().beta;

   return res;

 }

 // Impose a discontinuity G(\tau=0)-G(\tau=\beta)
 void enforce_discontinuity(gf_view<legendre> & gl, arrays::array_view<double,2> disc) {

   double norm = 0.0;
   arrays::vector<double> t(gl.data().shape()[0]);
   for (int i=0; i<t.size(); ++i) {
     t(i) = triqs::utility::legendre_t(i,1) / gl.domain().beta;
     norm += t(i)*t(i);
   }

   arrays::array<dcomplex, 2> corr(disc.shape());
   corr() = 0;
   for (auto const &l : gl.mesh()) corr += t(l.index()) * gl[l];

   auto _ = arrays::range{};
   for (auto const& l : gl.mesh()) gl.data()(l.index(), _, _) += (disc - corr) * t(l.index()) / norm;
 }

 auto curry_impl(gf_view<cartesian_product<Ms...>, Target, Singularity, Evaluator, IsConst> g) {
  // pick up the meshed corresponding to the curryed variables
  auto meshes_tuple = triqs::tuple::filter<pos...>(g.mesh().components());
  using var_t = cart_prod<triqs::tuple::filter_t<std::tuple<Ms...>, pos...>>;
  auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple);
  auto l = [g](auto&&... x) { return partial_eval_linear_index<pos...>(g, std::make_tuple(x...)); };
  return make_gf_view_lambda_valued<var_t>(m, l);
 };

 static auto invoke(gf_view<cartesian_product<Ms...>, Target, Singularity, Evaluator, IsConst> g, XTuple const& x_tuple) {
  using var_t = cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>;
  // meshes of the returned gf_view : just drop the mesh of the evaluated variables
  auto meshes_tuple_partial = triqs::tuple::filter_out<pos...>(g.mesh().components());
  // The mesh of the resulting function
  auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple_partial);
  // rebuild a tuple of the size sizeof...(Ms), containing the linear indices and range at the position of evaluated variables.
  auto arr_args = triqs::tuple::inverse_filter<sizeof...(Ms), pos...>(x_tuple, arrays::range());
  // from it, we make a slice of the array of g, corresponding to the data of the returned gf_view
  auto arr2 = triqs::tuple::apply(g.data(), std::tuple_cat(arr_args, std::make_tuple(arrays::ellipsis{})));
  // We now also partial_eval the singularity
  auto singv = partial_eval_linear_index<pos...>(g.singularity(), x_tuple);
  using r_sing_t = typename decltype(singv)::regular_type;
  // finally, we build the view on this data.
  using r_t = gf_view<var_t, Target, r_sing_t, void, IsConst>;
  return r_t{m, arr2, singv, {}, {}};
 }

 void legendre_matsubara_inverse(gf_view<legendre> gl, gf_const_view<imtime> gt) {

  gl() = 0.0;
  legendre_generator L;
  auto N = gt.mesh().size() - 1;
  double coef;

  // Do the integral over imaginary time
  for (auto t : gt.mesh()) {
   if (t.index()==0 || t.index()==N) coef = 0.5;
   else coef = 1.0;
   L.reset(2 * t / gt.domain().beta - 1);
   for (auto l : gl.mesh()) {
    gl[l] += coef * sqrt(2 * l.index() + 1) * L.next() * gt[t];
   }
  }
  gl.data() *= gt.mesh().delta();
 }

 // Impose a discontinuity G(\tau=0)-G(\tau=\beta)
 void enforce_discontinuity(gf_view<legendre> & gl, tqa::array_view<double,2> disc) {

   double norm = 0.0;
   tqa::vector<double> t(gl.data().shape()[0]);
   for (int i=0; i<t.size(); ++i) {
     t(i) = triqs::utility::legendre_t(i,1) / gl.domain().beta;
     norm += t(i)*t(i);
   }

   tqa::array<double,2> corr(disc.shape()); corr() = 0;
   for (auto l : gl.mesh()) {
     corr += t(l.index()) * gl[l];
   }

   tqa::range R;
   for (auto l : gl.mesh()) {
     gl.data()(l.index(),R,R) += (disc - corr) * t(l.index()) / norm;
   }

 }

 void pade (gf_view<refreq> &gr, gf_view<imfreq> const &gw, int n_points, double freq_offset) {

  // make sure the GFs have the same structure
  //assert(gw.shape() == gr.shape());

  // copy the tail. it doesn't need to conform to the pade approximant
  gr.singularity() = gw.singularity();

  auto sh = gw.data().shape().front_pop();
  int N1 = sh[0], N2 = sh[1];
  for (int n1=0; n1<N1; n1++) {
    for (int n2=0; n2<N2; n2++) {

      arrays::vector<dcomplex> z_in(n_points); // complex points
      arrays::vector<dcomplex> u_in(n_points); // values at these points
      arrays::vector<dcomplex> a(n_points);    // corresponding Pade coefficients

      for (int i=0; i < n_points; ++i) z_in(i) = gw.mesh()[i];
      for (int i=0; i < n_points; ++i) u_in(i) = gw.on_mesh(i)(n1,n2);

      triqs::utility::pade_approximant PA(z_in,u_in);

      gr() = 0.0;
      for (auto om : gr.mesh()) {
        dcomplex e = om + dcomplex(0.0,1.0)*freq_offset;
        gr[om](n1,n2) = PA(e);
      }

    }
  }

 }

 template <typename M, typename T, typename S, typename E> gf<M, real_target_t<T>, S> real(gf_view<M, T, S, E> g) {
  return {g.mesh(), real(g.data()), g.singularity(), g.symmetry(), {}, {}};
 }

void fill_data(gf_view<legendre> g_l, int i, vector<double> const& data) {
 g_l.data()(range(),i,i) = data;
}

template <typename T> size_t n_blocks(gf_view<block_index, T> const &g) { return g.mesh().size(); }

 tail fit_tail_impl(gf_view<imfreq> gf, const tail_view known_moments, int n_moments, int n_min, int n_max) {

  // precondition : check that n_max is not too large
  n_max = std::min(n_max, int(gf.mesh().size()-1));

  tail res(get_target_shape(gf));
  if (known_moments.size())
   for (int i = known_moments.order_min(); i <= known_moments.order_max(); i++) res(i) = known_moments(i);

  // if known_moments.size()==0, the lowest order to be obtained from the fit is determined by order_min in known_moments
  // if known_moments.size()==0, the lowest order is the one following order_max in known_moments

  int n_unknown_moments = n_moments - known_moments.size();
  if (n_unknown_moments < 1) return known_moments;

  // get the number of even unknown moments: it is n_unknown_moments/2+1 if the first
  // moment is even and n_moments is odd; n_unknown_moments/2 otherwise
  int omin = known_moments.size() == 0 ? known_moments.order_min() : known_moments.order_max() + 1; // smallest unknown moment
  int omin_even = omin % 2 == 0 ? omin : omin + 1;
  int omin_odd = omin % 2 != 0 ? omin : omin + 1;
  int size_even = n_unknown_moments / 2;
  if (n_unknown_moments % 2 != 0 && omin % 2 == 0) size_even += 1;
  int size_odd = n_unknown_moments - size_even;

  int size1 = n_max - n_min + 1;
  // size2 is the number of moments

  arrays::matrix<double> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT);
  arrays::matrix<double> B(size1, 1, FORTRAN_LAYOUT);
  arrays::vector<double> S(std::max(size_even, size_odd));
  const double rcond = 0.0;
  int rank;

  for (int i = 0; i < get_target_shape(gf)[0]; i++) {
   for (int j = 0; j < get_target_shape(gf)[1]; j++) {

    // fit the odd moments
    S.resize(size_odd);
    A.resize(size1,size_odd); //when resizing, gelss segfaults
    for (int k = 0; k < size1; k++) {
     auto n = n_min + k;
     auto iw = std::complex<double>(gf.mesh().index_to_point(n));

     B(k, 0) = imag(gf.data()(gf.mesh().index_to_linear(n), i, j));
     // subtract known tail if present
     if (known_moments.size() > 0)
      B(k, 0) -= imag(evaluate(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)), iw)(0, 0));

     for (int l = 0; l < size_odd; l++) {
      int order = omin_odd + 2 * l;
      A(k, l) = imag(pow(iw, -1.0 * order)); // set design matrix for odd moments
     }
    }

    arrays::lapack::gelss(A, B, S, rcond, rank);
    for (int m = 0; m < size_odd; m++) {
     res(omin_odd + 2 * m)(i, j) = B(m, 0);
    }

    // fit the even moments
    S.resize(size_even);
    A.resize(size1,size_even); //when resizing, gelss segfaults
    for (int k = 0; k < size1; k++) {
     auto n = n_min + k;
     auto iw = std::complex<double>(gf.mesh().index_to_point(n));

     B(k, 0) = real(gf.data()(gf.mesh().index_to_linear(n), i, j));
     // subtract known tail if present
     if (known_moments.size() > 0)
      B(k, 0) -= real(evaluate(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)), iw)(0, 0));

     for (int l = 0; l < size_even; l++) {
      int order = omin_even + 2 * l;
      A(k, l) = real(pow(iw, -1.0 * order)); // set design matrix for odd moments
     }
    }

    arrays::lapack::gelss(A, B, S, rcond, rank);
    for (int m = 0; m < size_even; m++) {
     res(omin_even + 2 * m)(i, j) = B(m, 0);
    }
   }
  }
  res.mask()()=n_moments;
  return res; // return tail
 }