From 3a9f986461ea64ce5098edc150a0fbf694336ce4 Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Sat, 10 May 2014 21:39:11 +0200 Subject: [PATCH] gf: add a few functions in C++ --- triqs/arrays/python/numpy_extractor.hpp | 5 +- triqs/gfs/gf.hpp | 94 ++++++++++++++++++++++++- triqs/gfs/local/tail.hpp | 37 ++++++---- 3 files changed, 118 insertions(+), 18 deletions(-) diff --git a/triqs/arrays/python/numpy_extractor.hpp b/triqs/arrays/python/numpy_extractor.hpp index 1d617961..c1d2925b 100644 --- a/triqs/arrays/python/numpy_extractor.hpp +++ b/triqs/arrays/python/numpy_extractor.hpp @@ -106,7 +106,10 @@ namespace triqs { namespace arrays { namespace numpy_interface { // do several checks if (!numpy_obj) {// The convertion of X to a numpy has failed ! - if (PyErr_Occurred()) {PyErr_Print();PyErr_Clear();} + if (PyErr_Occurred()) { + //PyErr_Print(); + PyErr_Clear(); + } TRIQS_RUNTIME_ERROR<<"numpy interface : the python object is not convertible to a numpy. "; } assert (PyArray_Check(numpy_obj)); assert((numpy_obj->ob_refcnt==1) || ((numpy_obj ==X))); diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp index e134086b..47d55288 100644 --- a/triqs/gfs/gf.hpp +++ b/triqs/gfs/gf.hpp @@ -28,11 +28,16 @@ #include #include "./tools.hpp" #include "./data_proxies.hpp" +#include "./local/tail.hpp" namespace triqs { namespace gfs { using utility::factory; using arrays::make_shape; + using arrays::array; + using arrays::array_view; + using arrays::matrix; + using arrays::matrix_view; using triqs::make_clone; // the gf mesh @@ -424,8 +429,8 @@ namespace gfs { using target_shape_t = typename factory::target_shape_t; - gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}, typename B::indices_t const &ind = typename B::indices_t{}) - : B(std::move(m), factory::make_data(m, shape), factory::make_singularity(m, shape), typename B::symmetry_t{}, ind, {}, // clean unncessary types + gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}, typename B::indices_t const &ind = typename B::indices_t{}, std::string name = "") + : B(std::move(m), factory::make_data(m, shape), factory::make_singularity(m, shape), typename B::symmetry_t{}, ind, name, // clean unncessary types typename B::evaluator_t{}) {} friend void swap(gf &a, gf &b) noexcept { a.swap_impl(b); } @@ -647,6 +652,91 @@ namespace gfs { g.singularity(), g.symmetry()); }*/ + // ---------------------------------- some functions : invert, conjugate, transpose, ... ------------------------------------ + + // ---- inversion + // auxiliary function : invert the data : one function for all matrix valued gf (save code). + template void _gf_invert_data_in_place(A3 && a) { + for (int i = 0; i < first_dim(a); ++i) {// Rely on the ordering + auto v = a(i, arrays::range(), arrays::range()); + v = inverse(v); + } + } + + template + void invert_in_place(gf_view g) { + _gf_invert_data_in_place(g.data()); + g.singularity() = inverse(g.singularity()); + } + + template gf inverse(gf_view g) { + auto res = gf(g); + invert_in_place(res); + return res; + } + + // ---- transpose : a new view +/* + template + gf_view invert_in_place(gf_view g) { + return { g.mesh(), transpose( g.data(), xxxx), transpose(g.singularity()), g.symmetry(), g.indices() XXX, g.name}; + } +*/ + + // ---- conjugate : always a new function -> changelog + + template gf conj(gf_view g) { + return {g.mesh(), conj(g.data()), conj(g.singularity()), g.symmetry(), g.indices(), g.name}; + } + + // ---- matrix mult R and L + + template void _gf_data_mul_R(A3 &&a, matrix const &r) { + for (int i = 0; i < first_dim(a); ++i) { // Rely on the ordering + matrix_view v = a(i, arrays::range(), arrays::range()); + v = v * r; + } + } + + template void _gf_data_mul_L(matrix const &l, A3 &&a) { + for (int i = 0; i < first_dim(a); ++i) { // Rely on the ordering + matrix_view v = a(i, arrays::range(), arrays::range()); + v = l * v; + } + } + + template typename A3::regular_type _gf_data_mul_LR(matrix const &l, A3 &a, matrix const &r) { + auto res = typename A3::regular_type(first_dim(a), first_dim(l), second_dim(r)); + for (int i = 0; i < first_dim(a); ++i) { // Rely on the ordering + auto _ = arrays::range{}; + res(i, _, _) = l * make_matrix_view(a(i, _, _))* r; + } + return res; + } + + template + gf operator*(gf g, matrix r) { + _gf_data_mul_R(g.data(), r); + g.singularity() = g.singularity() * r; + return g; + } + + template + gf operator*(matrix l, gf g) { + _gf_data_mul_L(l, g.data()); + g.singularity() = l * g.singularity(); + return g; + } + + template + gf L_G_R(matrix l, gf g, matrix r) { + auto res = gf{g.mesh(), {first_dim(l), second_dim(r)}}; + res.data() = _gf_data_mul_LR(l, g.data(), r); + res.singularity() = l * g.singularity() * r; + return res; + } + + namespace gfs_implementation { // implement some default traits // ------------------------- default factories --------------------- diff --git a/triqs/gfs/local/tail.hpp b/triqs/gfs/local/tail.hpp index 6a62af0a..d017e578 100644 --- a/triqs/gfs/local/tail.hpp +++ b/triqs/gfs/local/tail.hpp @@ -31,7 +31,9 @@ namespace triqs { namespace gfs { namespace local { static constexpr double small = 1.e-10; } - namespace tqa= triqs::arrays; namespace tql= triqs::clef; namespace mpl= boost::mpl; + using arrays::matrix_view; + using arrays::matrix; + namespace tqa= triqs::arrays; //namespace tql= triqs::clef; namespace mpl= boost::mpl; typedef std::complex dcomplex; class tail; // the value class @@ -40,7 +42,6 @@ namespace triqs { namespace gfs { namespace local { template struct LocalTail : mpl::false_{}; // a boolean trait to identify the objects modelling the concept LocalTail template<> struct LocalTail : mpl::true_{}; template<> struct LocalTail: mpl::true_{}; - template<> struct LocalTail>: mpl::true_{}; // a trait to find the scalar of the algebra i.e. the true scalar and the matrix ... template struct is_scalar_or_element : mpl::or_< tqa::ImmutableMatrix, utility::is_in_ZRC > {}; @@ -80,8 +81,6 @@ namespace triqs { namespace gfs { namespace local { typedef tqa::mini_vector shape_type; shape_type shape() const { return shape_type(_data.shape()[1], _data.shape()[2]);} - size_t shape(int i) const { return _data.shape()[i];} - bool is_decreasing_at_infinity() const { return (smallest_nonzero() >=1);} protected: @@ -231,6 +230,7 @@ namespace triqs { namespace gfs { namespace local { typedef tqa::mini_vector shape_type; tail(size_t N1, size_t N2, size_t size_=10, long order_min=-1): B(N1,N2,size_,order_min) {} tail(shape_type const & sh, size_t size_=10, long order_min=-1): B(sh[0],sh[1],size_,order_min) {} + tail(B::data_type const &d, B::mask_type const &m, long order_min): B(d, m, order_min) {} tail(tail const & g): B(g) {} tail(tail_view const & g): B(g) {} tail(tail &&) = default; @@ -279,6 +279,8 @@ namespace triqs { namespace gfs { namespace local { return *this; } + inline tail conj(tail_view t) { return {conj(t.data()), t.mask_view(),t.order_min()};} + /// Slice in orbital space //template tail_view slice_target(tail_impl const & t, tqa::range R1, tqa::range R2) { inline tail_view slice_target(tail_view t, tqa::range R1, tqa::range R2) { @@ -326,19 +328,24 @@ namespace triqs { namespace gfs { namespace local { return res; } - template - TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator* (T1 const & a, T2 const & b) { return mult_impl(a,b); } + template + TYPE_ENABLE_IF(tail, mpl::and_, LocalTail>) operator*(T1 const &a, T2 const &b) { + return mult_impl(a, b); + } - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator* (T1 const & a, T2 const & b) { - tail res(b); for (long n=res.order_min(); n<=res.order_max(); ++n) res(n)=a*res(n); return res; - } + template + TYPE_ENABLE_IF(tail, mpl::and_, LocalTail>) operator*(T1 const &a, T2 const &b) { + auto res = tail{first_dim(a), b.shape()[1], b.size(), b.order_min()}; + for (int n = res.order_min(); n <= res.order_max(); ++n) res(n) = a * b(n); + return res; + } - template TYPE_ENABLE_IF(tail,mpl::and_, tqa::ImmutableMatrix>) - operator* (T1 const & a, T2 const & b) { - tail res(a); for (long n=res.order_min(); n<=res.order_max(); ++n) res(n)=res(n)*b; return res; - } + template + TYPE_ENABLE_IF(tail, mpl::and_, tqa::ImmutableMatrix>) operator*(T1 const &a, T2 const &b) { + auto res = tail{a.shape()[0], second_dim(b), a.size(), a.order_min()}; + for (int n = res.order_min(); n <= res.order_max(); ++n) res(n) = a(n) * b; + return res; + } inline tail operator * (dcomplex a, tail_view const & r) { tail res(r); res.data()*=a; return res;} inline tail operator * (tail_view const & r, dcomplex a) { return a*r; }