3
0
mirror of https://github.com/triqs/dft_tools synced 2024-12-25 13:53:40 +01:00

array & gf transpose

- implement transposed_view for arrays.
- .transpose method for gf
- wrapped to python
- add call op. for GfImTime, using C++
- Added ChangeLog
- rm matrix_stack
- start cleaning old code
This commit is contained in:
Olivier Parcollet 2014-05-22 14:51:26 +02:00
parent 5c8270c437
commit cfe3532c94
32 changed files with 109 additions and 1258 deletions

10
ChangeLog Normal file
View File

@ -0,0 +1,10 @@
From version 1.0 to 1.2
Python
----------
gf:
* transpose method now returns a NEW green function, like conjugate, etc...

View File

@ -1,31 +0,0 @@
############### Density #########################
cdef extern from "triqs/gfs/local/functions.hpp":
matrix_view density(gf_imfreq &)
matrix_view density(gf_legendre &)
void enforce_discontinuity(gf_legendre &, array_view[double,TWO])
############### Fourier #########################
cdef extern from "triqs/gfs/local/fourier_matsubara.hpp" :
gf_imfreq fourier (gf_imtime & )
gf_imtime inverse_fourier (gf_imfreq & )
cdef extern from "triqs/gfs/local/fourier_real.hpp" :
gf_refreq fourier (gf_retime & )
gf_retime inverse_fourier (gf_refreq & )
#gf_refreq fourier (gf_retime & ) except +
#gf_retime inverse_fourier (gf_refreq & ) except +
############### Legendre #########################
cdef extern from "triqs/gfs/local/legendre_matsubara.hpp" :
gf_imfreq legendre_to_imfreq (gf_legendre &)
gf_imtime legendre_to_imtime (gf_legendre &)
gf_legendre imfreq_to_legendre (gf_imfreq &)
gf_legendre imtime_to_legendre (gf_imtime &)
############### Pade #########################
cdef extern from "triqs/gfs/local/pade.hpp":
void pade(gf_refreq &, gf_imfreq &, int, double) except +

View File

@ -1,62 +0,0 @@
from dcomplex cimport *
from arrays cimport *
from libcpp.vector cimport vector
from libcpp.string cimport string as std_string
from shared_ptr cimport *
from extractor cimport *
from h5 cimport *
# -------------------- Some generic tools -------------------------------
cdef extern from "triqs/gfs/tools.hpp" namespace "triqs::gfs" :
cdef enum statistic_enum "triqs::gfs::statistic_enum" :
Boson,Fermion
cdef cppclass nothing :
nothing ()
cdef cppclass freq_infty:
freq_infty ()
cdef cppclass mesh_pt_generator "triqs::gfs::mesh_pt_generator" [MeshType] :
mesh_pt_generator( MeshType * )
mesh_pt_generator()
complex to_point()
mesh_pt_generator operator++()
bint at_end()
void increment()
cdef cppclass indices_2_t :
indices_2_t( )
indices_2_t( vector[vector[std_string]] &)
vector[std_string] & operator[](int)
vector[vector[std_string]] & operator()()
bint same()
cdef extern from "triqs/gfs/meshes/linear.hpp" namespace "triqs::gfs" :
cdef enum mesh_enum "triqs::gfs::mesh_kind":
half_bins, full_bins, without_last
cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
cdef cppclass discrete_domain :
discrete_domain ()
discrete_domain (vector[std_string] &)
vector[std_string] & names()
cdef cppclass discrete_mesh :
discrete_mesh ()
discrete_mesh (discrete_mesh &)
discrete_domain & domain()
long size()
bint operator ==( discrete_mesh &)
include "tail.pxd"
include "imfreq.pxd"
include "imtime.pxd"
include "refreq.pxd"
include "retime.pxd"
include "legendre.pxd"
#include "two_real_times.pxd"

View File

@ -1,26 +0,0 @@
#!python
#cython: embedsignature=True
from cython.operator cimport dereference as deref, preincrement as inc #dereference and increment operators
cimport cython
import numpy
import string
import warnings
from block_gf import BlockGf
from math import pi
from h5 cimport *
include "mesh_imfreq.pyx"
include "mesh_imtime.pyx"
include "mesh_refreq.pyx"
include "mesh_retime.pyx"
#include "mesh_two_real_times.pyx"
include "mesh_legendre.pyx"
include "imfreq.pyx"
include "imtime.pyx"
include "refreq.pyx"
include "retime.pyx"
#include "two_real_times.pyx"
include "legendre.pyx"
include "tail.pyx"
include "functions.pxd"

View File

@ -1,57 +0,0 @@
cdef extern from "triqs/gfs/imfreq.hpp" namespace "triqs::gfs" :
cdef cppclass imfreq_domain :
double beta
statistic_enum statistic
imfreq_domain ()
cdef cppclass mesh_imfreq "triqs::gfs::gf_mesh<triqs::gfs::imfreq>" :
mesh_imfreq ()
mesh_imfreq (mesh_imfreq &)
imfreq_domain & domain()
double x_min()
long size()
bint positive_only()
bint operator ==( mesh_imfreq &)
cdef mesh_imfreq make_mesh_imfreq "triqs::gfs::gf_mesh<triqs::gfs::imfreq>" (double beta, statistic_enum S, size_t n_max, bint positive_only)
#cdef mesh_imfreq make_mesh_imfreq "triqs::gfs::gf_implementation::gf_factories<triqs::gfs::imfreq,triqs::gfs::matrix>::make_mesh" (double beta, statistic_enum S, size_t n_max)
cdef cppclass gf_imfreq "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::imfreq> >" :
gf_imfreq()
gf_imfreq(gf_imfreq &)
# The constructor must be no_except, or the cython code won't be correct...
gf_imfreq(mesh_imfreq, array_view[dcomplex, THREE], tail, nothing) #except +
void operator << (gf_imfreq &)
mesh_imfreq mesh()
array_view[dcomplex, THREE] data()
tail singularity()
cdef extern from "triqs/gfs/imfreq.hpp" :
cdef void h5_write (h5_group, char *, gf_imfreq &)
cdef extern from "triqs/utility/serialization.hpp" :
cdef std_string boost_serialize "triqs::serialize" (gf_imfreq &)
cdef void boost_unserialize_into "triqs::deserialize_into_view" (std_string, gf_imfreq &)
# Python -> C
cdef gf_imfreq as_gf_imfreq (g) except +
# C -> Python
cdef make_GfImFreq (gf_imfreq x, indices_pack=*, name=*)
############### Blocks of Im Freq #########################
cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
cdef cppclass gf_block_imfreq "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::imfreq>>>" :
gf_block_imfreq()
gf_imfreq & operator [](int)
discrete_mesh & mesh()
cdef gf_block_imfreq make_gf_block_imfreq "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::imfreq>>" (vector[gf_imfreq] &)
#cdef gf_block_imfreq make_gf_block_imfreq "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::imfreq>>" ( vector[gf_imfreq] &)
cdef gf_block_imfreq as_gf_block_imfreq (G) except +
cdef make_BlockGfImFreq (gf_block_imfreq G, block_indices_pack=*, name=*)

View File

@ -1,58 +0,0 @@
cdef extern from "triqs/gfs/imtime.hpp" namespace "triqs::gfs" :
cdef cppclass imtime_domain :
double beta
statistic_enum statistic
imtime_domain ()
cdef cppclass mesh_imtime "triqs::gfs::gf_mesh<triqs::gfs::imtime>" :
mesh_imtime ()
mesh_imtime (mesh_imtime &)
imtime_domain & domain()
double x_min()
long size()
long kind()
bint operator ==( mesh_imtime &)
cdef mesh_imtime make_mesh_imtime "triqs::gfs::gf_mesh<triqs::gfs::imtime>" (double beta, statistic_enum S, size_t n_time_slices, mesh_enum mk)
#cdef mesh_imtime make_mesh_imtime "triqs::gfs::gf_factories<triqs::gfs::imtime>::make_mesh" (double beta, statistic_enum S, size_t n_time_slices, mesh_enum mk)
cdef cppclass gf_imtime "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::imtime>>" :
gf_imtime()
gf_imtime(gf_imtime &)
# The constructor must be no_except, or the cython code won't be correct...
gf_imtime(mesh_imtime, array_view[double, THREE], tail, nothing) #except +
void operator << (gf_imtime &)
mesh_imtime mesh()
array_view[double, THREE] data()
matrix_view[double] operator()(double)
tail singularity()
cdef extern from "triqs/gfs/imtime.hpp" :
cdef void h5_write (h5_group, char *, gf_imtime &)
cdef extern from "triqs/utility/serialization.hpp" :
cdef std_string boost_serialize "triqs::serialize" (gf_imtime &)
cdef void boost_unserialize_into "triqs::deserialize_into_view" (std_string, gf_imtime &)
# Python -> C
cdef gf_imtime as_gf_imtime (g) except +
# C -> Python
cdef make_GfImTime (gf_imtime x, indices_pack=*, name=*)
############### Blocks of Im Time #########################
cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
cdef cppclass gf_block_imtime "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::imtime>>>" :
gf_block_imtime()
gf_imtime & operator [](int)
discrete_mesh & mesh()
cdef gf_block_imtime make_gf_block_imtime "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::imtime>>" (vector[gf_imtime] &)
#cdef gf_block_imtime make_gf_block_imtime "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::imtime>>" ( vector[gf_imtime] &)
cdef gf_block_imtime as_gf_block_imtime (G) except +
cdef make_BlockGfImTime (gf_block_imtime G, block_indices_pack=*, name=*)

View File

@ -1,54 +0,0 @@
cdef extern from "triqs/gfs/legendre.hpp" namespace "triqs::gfs" :
cdef cppclass legendre_domain :
size_t Nmax
double beta
statistic_enum statistic
legendre_domain ()
cdef cppclass mesh_legendre "triqs::gfs::gf_mesh<triqs::gfs::legendre>" :
mesh_legendre ()
mesh_legendre (mesh_legendre &)
legendre_domain & domain()
long size()
bint operator == (mesh_legendre &)
cdef mesh_legendre make_mesh_legendre "triqs::gfs::gf_mesh<triqs::gfs::legendre>" (double beta, statistic_enum S, size_t n_leg)
#cdef mesh_legendre make_mesh_legendre "triqs::gfs::gf_factories<triqs::gfs::legendre>::make_mesh" (double beta, statistic_enum S, size_t n_leg)
cdef cppclass gf_legendre "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::legendre>>" :
gf_legendre()
gf_legendre(gf_legendre &)
gf_legendre(mesh_legendre, array_view[double, THREE], nothing, nothing) #except +
void operator << (gf_legendre &)
mesh_legendre mesh()
array_view[double, THREE] data()
cdef extern from "triqs/gfs/legendre.hpp" :
cdef void h5_write (h5_group, char *, gf_legendre &)
cdef extern from "triqs/utility/serialization.hpp" :
cdef std_string boost_serialize "triqs::serialize" (gf_legendre &)
cdef void boost_unserialize_into "triqs::deserialize_into_view" (std_string, gf_legendre &)
# Python -> C
cdef gf_legendre as_gf_legendre(g) except +
# C -> Python
cdef make_GfLegendre(gf_legendre x, indices_pack=*, name=*)
############### Blocks of Im Time #########################
cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
cdef cppclass gf_block_legendre "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::legendre>>>" :
gf_block_legendre()
gf_legendre & operator [](int)
discrete_mesh & mesh()
cdef gf_block_legendre make_gf_block_legendre "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::legendre>>" (vector[gf_legendre] &)
#cdef gf_block_legendre make_gf_block_legendre "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::legendre>>" ( vector[gf_legendre] &)
cdef gf_block_legendre as_gf_block_legendre (G) except +
cdef make_BlockGfLegendre (gf_block_legendre G, block_indices_pack=*, name=*)

View File

@ -1,17 +0,0 @@
include "arrays.pxd"
include "dcomplex.pxd"
cdef extern from "<triqs/arrays/matrix_stack_view.hpp>" namespace "triqs::arrays":
cdef cppclass matrix_stack "triqs::arrays::matrix_stack_view" [T]:
matrix_stack(object)
void iadd "operator +=" (matrix_view[T] &) except +
void isub "operator -=" (matrix_view[T] &) except +
void imul "operator *=" (matrix_view[T] &) except +
void idiv "operator /=" (matrix_view[T] &) except +
void invert() except +
void onsite_matmul_L_R(matrix_view[T] &, matrix_stack &, matrix_view[T] &) except +

View File

@ -1,34 +0,0 @@
import numpy
class MatrixStack:
def __init__(self, A):
"""
"""
self.data = A
self.is_double = A.dtype in [numpy.float, numpy.float32, numpy.float64, numpy.float_]
def add(self, M):
if self.is_double:
matrix_stack[double](self.data).iadd(matrix_view[double](M))
else:
matrix_stack[dcomplex](self.data).iadd(matrix_view[dcomplex](M.astype(numpy.complex)))
def sub(self, M):
if self.is_double:
matrix_stack[double](self.data).isub(matrix_view[double](M))
else:
matrix_stack[dcomplex](self.data).isub(matrix_view[dcomplex](M.astype(numpy.complex)))
def matmul_L_R(self, L, D, R):
if self.is_double:
matrix_stack[double](self.data).onsite_matmul_L_R(matrix_view[double](L), matrix_stack[double](D), matrix_view[double](R))
else:
matrix_stack[dcomplex](self.data).onsite_matmul_L_R(
matrix_view[dcomplex](L.astype(numpy.complex)), matrix_stack[dcomplex](D), matrix_view[dcomplex](R.astype(numpy.complex)))
def invert(self):
if self.is_double:
matrix_stack[double](self.data).invert()
else:
matrix_stack[dcomplex](self.data).invert()

View File

@ -1,36 +0,0 @@
cdef class MeshImFreq:
cdef mesh_imfreq _c
def __init__(self, beta, stat, int n_pts, bint positive_only=True):
self._c = make_mesh_imfreq(beta,{'F' :Fermion, 'B' : Boson}[stat], n_pts, positive_only)
def __len__ (self) : return self._c.size()
property beta :
"""Inverse temperature"""
def __get__(self): return self._c.domain().beta
property statistic :
def __get__(self): return 'F' if self._c.domain().statistic==Fermion else 'B'
property positive_only :
def __get__(self): return self._c.positive_only()
def __iter__(self) : # I use the C++ generator !
cdef mesh_pt_generator[mesh_imfreq ] g = mesh_pt_generator[mesh_imfreq ](&self._c)
while not g.at_end() :
yield g.to_point()
g.increment()
def __richcmp__(MeshImFreq self, MeshImFreq other,int op) :
if op ==2 : # ==
return self._c == other._c
def __reduce__(self):
return self.__class__, (self.beta, self.statistic, len(self))
# C -> Python
cdef inline make_MeshImFreq (mesh_imfreq x) :
return MeshImFreq( x.domain().beta, 'F' if x.domain().statistic==Fermion else 'B', x.size(), x.positive_only() )

View File

@ -1,37 +0,0 @@
cdef class MeshImTime:
cdef mesh_imtime _c
def __init__(self, beta, stat, int n_max, kind):
self._c = make_mesh_imtime(beta, {'F' :Fermion, 'B' : Boson}[stat], n_max, {'H':half_bins, 'F':full_bins, 'W':without_last}[kind])
def __len__ (self) : return self._c.size()
property beta :
"""Inverse temperature"""
def __get__(self): return self._c.domain().beta
property kind :
"""Inverse temperature"""
def __get__(self): return {half_bins: 'H', full_bins: 'F', without_last: 'W'}[self._c.kind()]
property statistic :
def __get__(self): return 'F' if self._c.domain().statistic==Fermion else 'B'
def __iter__(self) : # I use the C++ generator !
cdef mesh_pt_generator[mesh_imtime ] g = mesh_pt_generator[mesh_imtime ](&self._c)
while not g.at_end() :
yield g.to_point()
g.increment()
def __richcmp__(MeshImTime self, MeshImTime other,int op) :
if op ==2 : # ==
return self._c == other._c
def __reduce__(self):
return self.__class__, (self.beta, self.statistic, len(self), self.kind)
# C -> Python
cdef inline make_MeshImTime ( mesh_imtime x) :
return MeshImTime( x.domain().beta, 'F' if x.domain().statistic==Fermion else 'B',
x.size(), {half_bins: 'H', full_bins: 'F', without_last: 'W'}[x.kind()] )

View File

@ -1,34 +0,0 @@
cdef class MeshReFreq:
cdef mesh_refreq _c
def __init__(self, double omega_min, double omega_max, int n_max, kind):
self._c = make_mesh_refreq(omega_min, omega_max, n_max, {'H':half_bins, 'F':full_bins, 'W':without_last}[kind])
property omega_min:
def __get__(self): return self._c.x_min()
property omega_max:
def __get__(self): return self._c.x_max()
property kind:
def __get__(self): return {half_bins: 'H', full_bins: 'F', without_last: 'W'}[self._c.kind()]
def __len__ (self) : return self._c.size()
def __iter__(self) : # I use the C++ generator !
cdef mesh_pt_generator[mesh_refreq] g = mesh_pt_generator[mesh_refreq](&self._c)
while not g.at_end() :
yield g.to_point()
g.increment()
def __richcmp__(MeshReFreq self, MeshReFreq other, int op) :
if op ==2 : # ==
return self._c == other._c
def __reduce__(self):
return self.__class__, (self.omega_min, self.omega_max, len(self), self.kind)
# C -> Python
cdef inline make_MeshReFreq ( mesh_refreq x) :
return MeshReFreq( x.x_min(), x.x_max(), x.size(), {half_bins: 'H', full_bins: 'F', without_last: 'W'}[x.kind()] )

View File

@ -1,34 +0,0 @@
cdef class MeshReTime:
cdef mesh_retime _c
def __init__(self, double t_min, double t_max, int n_max, kind):
self._c = make_mesh_retime(t_min, t_max, n_max, {'H':half_bins, 'F':full_bins, 'W':without_last}[kind])
property t_min:
def __get__(self): return self._c.x_min()
property t_max:
def __get__(self): return self._c.x_max()
property kind:
def __get__(self): return {half_bins: 'H', full_bins: 'F', without_last: 'W'}[self._c.kind()]
def __len__ (self) : return self._c.size()
def __iter__(self) : # I use the C++ generator !
cdef mesh_pt_generator[mesh_retime] g = mesh_pt_generator[mesh_retime](&self._c)
while not g.at_end() :
yield g.to_point()
g.increment()
def __richcmp__(MeshReTime self, MeshReTime other, int op) :
if op ==2 : # ==
return self._c == other._c
def __reduce__(self):
return self.__class__, (self.t_min, self.t_max, len(self), self.kind)
# C -> Python
cdef inline make_MeshReTime ( mesh_retime x) :
return MeshReTime( x.x_min(), x.x_max(), x.size(), {half_bins: 'H', full_bins: 'F', without_last: 'W'}[x.kind()] )

View File

@ -1,43 +0,0 @@
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011-2012 by M. Ferrero, O. Parcollet
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
class Nothing (object):
def __init__(self): pass
def __getitem__(self, key): return Nothing()
def __setitem__(self, key, val): pass
def __iadd__(self, arg): return Nothing()
def __add__(self, y): return Nothing()
def __radd__(self, y): return Nothing()
def __isub__(self, arg): return Nothing()
def __sub__(self, y): return Nothing()
def __rsub__(self, y): return Nothing()
def __imul__(self, arg): return Nothing()
def __mul__(self, arg): return Nothing()
def __rmul__(self, x): return Nothing()
def __idiv__(self, arg): return Nothing()
def __div__(self, arg): return Nothing()
def _make_slice(self, sl1, sl2): return Nothing()
def invert(self): pass
def transpose(self): return Nothing()
def conjugate(self): return Nothing()
def copy(self): return Nothing()

View File

@ -1,56 +0,0 @@
cdef extern from "triqs/gfs/refreq.hpp" namespace "triqs::gfs" :
cdef cppclass refreq_domain :
refreq_domain()
cdef cppclass mesh_refreq "triqs::gfs::gf_mesh<triqs::gfs::refreq>" :
mesh_refreq ()
mesh_refreq (mesh_refreq &)
refreq_domain & domain()
double x_min()
double x_max()
long size()
long kind()
bint operator ==( mesh_refreq &)
cdef mesh_refreq make_mesh_refreq "triqs::gfs::gf_mesh<triqs::gfs::refreq>" (double omega_min, double omega_max, size_t n_freq, mesh_enum mk)
#cdef mesh_refreq make_mesh_refreq "triqs::gfs::gf_factories<triqs::gfs::refreq>::make_mesh" (double omega_min, double omega_max, size_t n_freq, mesh_enum mk)
cdef cppclass gf_refreq "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::refreq>>" :
gf_refreq()
gf_refreq(gf_refreq &)
# The constructor must be no_except, or the cython code won't be correct...
gf_refreq(mesh_refreq, array_view[dcomplex, THREE], tail, nothing) #except +
void operator << (gf_refreq &)
mesh_refreq mesh()
array_view[dcomplex, THREE] data()
tail singularity()
cdef extern from "triqs/gfs/refreq.hpp" :
cdef void h5_write (h5_group, char *, gf_refreq &)
cdef extern from "triqs/utility/serialization.hpp" :
cdef std_string boost_serialize "triqs::serialize" (gf_refreq &)
cdef void boost_unserialize_into "triqs::deserialize_into_view" (std_string, gf_refreq &)
# Python -> C
cdef gf_refreq as_gf_refreq (g) except +
# C -> Python
cdef make_GfReFreq (gf_refreq x, indices_pack=*, name=*)
############### Blocks of Im Time #########################
cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
cdef cppclass gf_block_refreq "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::refreq>>>" :
gf_block_refreq()
gf_refreq & operator [](int)
discrete_mesh & mesh()
cdef gf_block_refreq make_gf_block_refreq "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::refreq>>" (vector[gf_refreq] &)
#cdef gf_block_refreq make_gf_block_refreq "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::refreq>>" ( vector[gf_refreq] &)
cdef gf_block_refreq as_gf_block_refreq (G) except +
cdef make_BlockGfReFreq (gf_block_refreq G, block_indices_pack=*, name=*)

View File

@ -1,79 +0,0 @@
from gf_refreq import GfReFreq
cdef class GfReFreq_cython:
cdef gf_refreq _c
def __init__(self, MeshReFreq mesh, data, TailGf tail):
self._c = gf_refreq(mesh._c, array_view[dcomplex,THREE](data), tail._c, nothing())
def __write_hdf5_cython__ (self, gr , char * key) :
h5_write (make_h5_group(gr), key, self._c)
def set_from_fourier(self, GfReTime_cython gt) :
"""Fills self with the Fourier transform of gt"""
self._c << fourier( gt._c )
# put if back with make_gf_from_fourier when approved
#def inverse_fourier(self):
# return make_GfReTime(inverse_fourier(self._c))
def set_from_pade(self, GfImFreq_cython gw, n_points = 100, freq_offset = 0.0) :
pade(self._c, gw._c, n_points, freq_offset)
def __dealloc__ (self):
pass
#---------------- Reading from h5 ---------------------------------------
def h5_read_GfReFreq(gr, key):
try:
indicesL = gr[key]['indices']['left']
indicesR = gr[key]['indices']['right']
pack = [indicesL, indicesR]
except:
pack = []
try:
name = gr[key]['name']
except:
name = key
return make_GfReFreq(h5_extractor[gf_refreq]()(make_h5_group(gr),key), pack, name)
from pytriqs.archive.hdf_archive_schemes import register_class
register_class (GfReFreq, read_fun = h5_read_GfReFreq)
#---------------- Convertions functions ---------------------------------------
# Python -> C
cdef gf_refreq as_gf_refreq (g) except +:
return (<GfReFreq_cython?>g)._c
# C -> Python. Do NOT add except +
cdef make_GfReFreq (gf_refreq x, indices_pack = [], name = "g"):
data = x.data().to_python()
if indices_pack == []:
indices_pack = [range(data.shape[1]), range(data.shape[2])]
return GfReFreq(
mesh = make_MeshReFreq (x.mesh()),
data = data,
tail = make_TailGf (x.singularity()),
indices_pack = indices_pack,
name = name)
# Python -> C for blocks
cdef gf_block_refreq as_gf_block_refreq (G) except +:
cdef vector[gf_refreq] v_c
for n,g in G:
v_c.push_back(as_gf_refreq(g))
return make_gf_block_refreq (v_c)
# C -> Python for block
cdef make_BlockGfReFreq (gf_block_refreq G, block_indices_pack = [], name = "G"):
gl = []
name_list = G.mesh().domain().names()
if block_indices_pack == []:
for i,n in enumerate(name_list):
sha = G[i].data().to_python().shape[1:3]
block_indices_pack.append( [range(sha[0]), range(sha[1])] )
for i,n in enumerate(name_list):
gl.append( make_GfReFreq(G[i], block_indices_pack[i]) )
return BlockGf( name_list = name_list, block_list = gl, name = name )

View File

@ -1,56 +0,0 @@
cdef extern from "triqs/gfs/retime.hpp" namespace "triqs::gfs" :
cdef cppclass retime_domain :
retime_domain()
cdef cppclass mesh_retime "triqs::gfs::gf_mesh<triqs::gfs::retime>" :
mesh_retime ()
mesh_retime (mesh_retime &)
retime_domain & domain()
double x_min()
double x_max()
long size()
long kind()
bint operator ==( mesh_retime &)
cdef mesh_retime make_mesh_retime "triqs::gfs::gf_mesh<triqs::gfs::retime>" (double t_min, double t_max, size_t n_freq, mesh_enum mk)
#cdef mesh_retime make_mesh_retime "triqs::gfs::gf_factories<triqs::gfs::retime>::make_mesh" (double t_min, double t_max, size_t n_freq, mesh_enum mk)
cdef cppclass gf_retime "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::retime>>" :
gf_retime()
gf_retime(gf_retime &)
# The constructor must be no_except, or the cython code won't be correct...
gf_retime(mesh_retime, array_view[dcomplex, THREE], tail, nothing) #except +
void operator << (gf_retime &)
mesh_retime mesh()
array_view[dcomplex, THREE] data()
tail singularity()
cdef extern from "triqs/gfs/retime.hpp" :
cdef void h5_write (h5_group, char *, gf_retime &)
cdef extern from "triqs/utility/serialization.hpp" :
cdef std_string boost_serialize "triqs::serialize" (gf_retime &)
cdef void boost_unserialize_into "triqs::deserialize_into_view" (std_string, gf_retime &)
# Python -> C
cdef gf_retime as_gf_retime (g) except +
# C -> Python
cdef make_GfReTime (gf_retime x, indices_pack=*, name=*)
############### Blocks of Im Time #########################
cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
cdef cppclass gf_block_retime "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::retime>>>" :
gf_block_retime()
gf_retime & operator [](int)
discrete_mesh & mesh()
cdef gf_block_retime make_gf_block_retime "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::retime>>" (vector[gf_retime] &)
#cdef gf_block_retime make_gf_block_retime "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::retime>>" ( vector[gf_retime] &)
cdef gf_block_retime as_gf_block_retime (G) except +
cdef make_BlockGfReTime (gf_block_retime G, block_indices_pack=*, name=*)

View File

@ -1,76 +0,0 @@
from gf_retime import GfReTime
cdef class GfReTime_cython:
cdef gf_retime _c
def __init__(self, MeshReTime mesh, data, TailGf tail):
self._c = gf_retime(mesh._c, array_view[dcomplex,THREE](data), tail._c, nothing())
def __write_hdf5_cython__ (self, gr , char * key) :
h5_write (make_h5_group(gr), key, self._c)
# Put it back ?
#def fourier(self):
# return make_GfReFreq(fourier(self._c))
def set_from_inverse_fourier(self, GfReFreq_cython gw) :
"""Fills self with the Inverse Fourier transform of gw"""
self._c << inverse_fourier( gw._c)
def __dealloc__ (self):
pass
#---------------- Reading from h5 ---------------------------------------
def h5_read_GfReTime(gr, key):
try:
indicesL = gr[key]['indices']['left']
indicesR = gr[key]['indices']['right']
pack = [indicesL, indicesR]
except:
pack = []
try:
name = gr[key]['name']
except:
name = key
return make_GfReTime(h5_extractor[gf_retime]()(make_h5_group(gr),key), pack, name)
from pytriqs.archive.hdf_archive_schemes import register_class
register_class (GfReTime, read_fun = h5_read_GfReTime)
#---------------- Convertions functions ---------------------------------------
# Python -> C
cdef gf_retime as_gf_retime (g) except +:
return (<GfReTime_cython?>g)._c
# C -> Python. Do NOT add except +
cdef make_GfReTime (gf_retime x, indices_pack = [], name = "g"):
data = x.data().to_python()
if indices_pack == []:
indices_pack = [range(data.shape[1]), range(data.shape[2])]
return GfReTime(
mesh = make_MeshReTime (x.mesh()),
data = data,
tail = make_TailGf (x.singularity()),
indices_pack = indices_pack,
name = name)
# Python -> C for blocks
cdef gf_block_retime as_gf_block_retime (G) except +:
cdef vector[gf_retime] v_c
for n,g in G:
v_c.push_back(as_gf_retime(g))
return make_gf_block_retime (v_c)
# C -> Python for block
cdef make_BlockGfReTime (gf_block_retime G, block_indices_pack = [], name = "G"):
gl = []
name_list = G.mesh().domain().names()
if block_indices_pack == []:
for i,n in enumerate(name_list):
sha = G[i].data().to_python().shape[1:3]
block_indices_pack.append( [range(sha[0]), range(sha[1])] )
for i,n in enumerate(name_list):
gl.append( make_GfReTime(G[i], block_indices_pack[i]) )
return BlockGf( name_list = name_list, block_list = gl, name = name )

View File

@ -1,37 +0,0 @@
from dcomplex cimport *
from arrays cimport *
cdef extern from "triqs/gfs/local/tail.hpp" :
cdef cppclass tail "triqs::python_tools::cython_proxy<triqs::gfs::local::tail_view>" :
tail()
tail(array_view[dcomplex,THREE], array_view[long,TWO], long) except +
matrix_view[dcomplex] operator()(int) except +
array_view[dcomplex,THREE] data()
array_view[long,TWO] mask_view() except +
void operator << (tail &) except +
void operator << (dcomplex &) except +
long order_min() except +
long order_max() except +
size_t size() except +
size_t shape(int) except +
cdef tail operator +( tail &, tail &) except +
cdef tail operator -( tail &, tail &) except +
cdef tail operator *( tail&, tail &) except +
cdef tail operator *( double, tail &) except +
cdef tail operator *( tail &, double) except +
cdef tail operator /( double, tail &) except +
cdef tail operator /( tail &, double) except +
cdef tail operator *( dcomplex, tail &) except +
cdef tail operator *( tail &, dcomplex) except +
cdef tail operator /( dcomplex, tail &) except +
cdef tail operator /( tail &, dcomplex) except +
cdef tail inverse_c "inverse" ( tail &) except +
cdef void h5_write (h5_group, char *, tail &)
cdef extern from "triqs/utility/serialization.hpp" :
cdef std_string boost_serialize "triqs::serialize" (tail &)
cdef void boost_unserialize_into "triqs::deserialize_into_view" (std_string, tail &)

View File

@ -8,19 +8,7 @@ cdef class TailGf:
TailGf ( shape )
TailGf ( data, mask, order_min )
"""
c_obj = d.pop('encapsulated_c_object', None)
if c_obj :
assert d == {}
self._c = extractor [tail] (c_obj) ()
return
bss = d.pop('boost_serialization_string', None)
if bss :
assert d == {}, "Internal error : boost_serialization_string must be the only argument"
boost_unserialize_into(<std_string>bss,self._c)
return
# default values
# default values
omin = -1
omax = 8
@ -39,50 +27,12 @@ cdef class TailGf:
self._c = tail(array_view[dcomplex,THREE](a), array_view[long,TWO](m), o)
#-------------- Reduction -------------------------------
def __reduce__(self):
return py_deserialize, (self.__class__,boost_serialize(self._c),)
#-------------- Properties -------------------------------
property data :
"""Access to the data array"""
def __get__(self) :
return self._c.data().to_python()
property mask:
"""Access to the mask"""
def __get__(self) :
return self._c.mask_view().to_python()
property shape :
def __get__(self) : return self.data.shape[1:3]
property order_min :
"""Min order of the expansion"""
def __get__(self) : return self._c.order_min()
property order_max :
"""Max order of the expansion"""
def __get__(self) : return self._c.order_max()
property N1 :
def __get__(self): return self.shape[0]
property N2 :
def __get__(self): return self.shape[1]
property size :
"""Length of the expansion"""
def __get__(self) : return self._c.size()
def copy(self) :
return self.__class__(data = self.data.copy(), mask = self.mask.copy())
def copy_from(self, TailGf T) :
self._c << T._c
def _make_slice(self, sl1, sl2):
return self.__class__(data = self.data[:,sl1,sl2], mask = self.mask[sl1,sl2])
@ -95,85 +45,13 @@ cdef class TailGf:
else:
return string.join([ "%s"%self[r]+ (" /" if r>0 else "") + " Om^%s"%(abs(r)) for r in range(omin, self.order_max+1) ] , " + ")
def __getitem__(self,i) :
"""Returns the i-th coefficient of the expansion, or order Om^i"""
if not self.has_coef(i) : raise IndexError, "Index %s is out of range"%i
return self.data[i-self.order_min,:,:]
def __setitem__(self,i, val) :
"""Sets the i-th coefficient of the expansion, or order Om^i"""
if not self.has_coef(i) : raise IndexError, "Index %s is out of range"%i
self.data[i-self.order_min,:,:] = val
def has_coef(self, i):
return (i >= self.order_min) and (i <= self.order_max)
def __call__(self, x) :
val = 0.0
for n in range(self.order_min, self.order_max+1):
val += self[n] * x**(-n)
return val
def invert(self) :
self._c << inverse_c (self._c)
######### arithmetic operations #################
def __iadd__(self, TailGf arg):
self._c << self._c + arg._c
return self
def __add__(self, TailGf y):
cdef tail t = (<TailGf?>self)._c + y._c
return make_TailGf(t)
#c = self.copy()
#c += y
#return c
def __isub__(self, TailGf arg):
self._c << self._c - arg._c
return self
def __sub__(self,TailGf y):
cdef tail t = (<TailGf?>self)._c - y._c
return make_TailGf(t)
#c = self.copy()
#c -= y
#return c
def __imul__(self,arg):
n = type(arg).__name__
if n == 'TailGf' :
self._c << self._c * (<TailGf?>arg)._c
elif descriptors.is_scalar(arg):
self._c << as_dcomplex(arg)* self._c
else :
raise RuntimeError, " argument type not recognized in imul for %s"%arg
return self
def __mul__(self,arg):
c = self.copy()
c *= arg
return c
def __idiv__(self,arg):
if descriptors.is_scalar(arg):
self._c << self._c / as_dcomplex(arg)
else:
raise RuntimeError, "rhs must be a scalar"
return self
def __div__(self,arg):
c = self.copy()
c /= arg
return c
#---- other operations ----
def zero(self) :
"""Sets the expansion to 0"""
self._c << as_dcomplex(0.0)
def transpose (self) :
"""Transpose the array : new view as in numpy"""
return TailGf(data=self.data.transpose(), mask=self.mask.transpose())
@ -182,21 +60,5 @@ cdef class TailGf:
"""Transpose the array : new view as in numpy"""
return TailGf(data=self.data.conjugate(), mask=self.mask)
def __write_hdf5__ (self, gr , char * key) :
h5_write (make_h5_group(gr), key, self._c)
#---------------- Reading from h5 ---------------------------------------
def h5_read_TailGf( gr, std_string key) :
return make_TailGf( h5_extractor[tail]()(make_h5_group(gr),key))
from pytriqs.archive.hdf_archive_schemes import register_class
register_class (TailGf, read_fun = h5_read_TailGf)
#-----------------------------------------------------
# C -> Python
#-----------------------------------------------------
cdef inline make_TailGf ( tail x) :
return TailGf(encapsulated_c_object = encapsulate (&x))

View File

@ -1,59 +0,0 @@
cdef extern from "triqs/gfs.hpp" namespace "triqs::gfs" :
cdef cppclass two_real_times_domain :
two_real_times_domain()
cdef cppclass mesh_two_real_times "triqs::gfs::gf_mesh<triqs::gfs::two_real_times>" :
mesh_two_real_times ()
mesh_two_real_times (mesh_two_real_times &)
two_real_times_domain & domain()
double x_min()
double x_max()
long size()
double kind()
bint operator ==( mesh_two_real_times &)
cdef mesh_two_real_times make_mesh_two_real_times "triqs::gfs::gf_mesh<triqs::gfs::two_real_times>" (double t_max, double n_time_slices)
cdef mesh_retime & get_1d_mesh_from_2times_mesh "triqs::gfs::get_1d_mesh_from_2times_mesh" (mesh_two_real_times &)
cdef cppclass gf_two_real_times "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::two_real_times>>" :
gf_two_real_times()
gf_two_real_times(gf_two_real_times &)
# The constructor must be no_except, or the cython code won't be correct...
gf_two_real_times(mesh_two_real_times, array_view[dcomplex, THREE], nothing, nothing) #except +
void operator << (gf_two_real_times &)
mesh_two_real_times mesh()
array_view[dcomplex, THREE] data()
tail singularity()
cdef gf_retime slice1d "triqs::gfs::slice" (gf_two_real_times &, double t) except +
cdef extern from "triqs/gfs/two_real_times.hpp" :
cdef void h5_write (h5_group, char *, gf_two_real_times &)
cdef extern from "triqs/utility/serialization.hpp" :
cdef std_string boost_serialize "triqs::serialize" (gf_two_real_times &)
cdef void boost_unserialize_into "triqs::deserialize_into_view" (std_string, gf_two_real_times &)
# Python -> C
cdef gf_two_real_times as_gf_two_real_times (g) except +
# C -> Python
cdef make_GfTwoRealTime (gf_two_real_times x, indices_pack=*, name=*)
############### Blocks of Im Time #########################
cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
cdef cppclass gf_block_two_real_times "triqs::python_tools::cython_proxy<triqs::gfs::gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::two_real_times>>>" :
gf_block_two_real_times()
gf_two_real_times & operator [](int)
discrete_mesh & mesh()
cdef gf_block_two_real_times make_gf_block_two_real_times "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::two_real_times>>" ( vector[gf_two_real_times] &)
#cdef gf_block_two_real_times make_gf_block_two_real_times "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::two_real_times>>" ( vector[gf_two_real_times] &)
cdef gf_block_two_real_times as_gf_block_two_real_times (G) except +
cdef make_BlockGfTwoRealTime (gf_block_two_real_times G, block_indices_pack=*, name=*)

View File

@ -122,14 +122,10 @@ def make_mesh( py_type, c_tag, has_kind=True, is_im=False) :
getter = cfunction(calling_pattern="mesh_kind result = self_c.kind()", signature = "mesh_kind()"),
doc = "")
#def __richcmp__(MeshImFreq self, MeshImFreq other,int op) :
# if op ==2 : # ==
# return self._c == other._c
#def __reduce__(self):
# return self.__class__, (self.beta, self.statistic, len(self))
return m
########################
@ -291,6 +287,11 @@ def make_gf( py_type, c_tag, is_complex_data = True, is_im = False) :
g.add_method(py_name = "invert", calling_pattern = "invert_in_place(self_c)" , signature = "void()", doc = "Invert (in place)")
g.add_method(py_name = "transpose",
calling_pattern = "auto result = transpose(self_c)",
signature = "gf<%s>()"%c_tag,
doc = "Returns a NEW gf, with transposed data...")
if c_tag != "imtime" :
g.add_method(py_name = "conjugate", calling_pattern = "auto result = conj(self_c)" , signature = "gf<%s>()"%c_tag, doc = "Return a new function, conjugate of self.")
@ -353,6 +354,9 @@ g.add_method(py_name = "set_from_inverse_fourier",
calling_pattern = "self_c = inverse_fourier(*gw)",
doc = """Fills self with the Inverse Fourier transform of gw""")
# add the call operator using the interpolation
g.add_call(signature = "matrix<double>(double tau)", doc = "G(tau) using interpolation")
module.add_class(g)
########################

View File

@ -1,89 +0,0 @@
#include "./common.hpp"
#include <triqs/arrays/matrix_stack_view.hpp>
using namespace triqs::arrays;
using namespace triqs;
int main() {
clef::placeholder<0> i_;
clef::placeholder<1> j_;
clef::placeholder<2> k_;
{
std::cerr<< "------------- Testing inversion ---------------"<< std::endl;
array <double,3> A(4,2,2);
A(k_,i_,j_) << clef::if_else(i_ == j_, i_ + k_+1, 0); // + 2.3* j_;
matrix_stack_view<double> S (std::move(A));
std::cout << S(0) << S(1)<< S(2) << std::endl ;
S.invert();
std::cout << S(0) << S(1)<< S(2)<< std::endl ;
}
{
std::cerr<< "------------- Testing += matrix ---------------"<< std::endl;
array <double,3> A(4,2,3);
matrix<double> M(2,3); M()= 1;
A(k_,i_,j_) << i_ + 2*j_ + k_/2.0;
matrix_stack_view<double> S (std::move(A));
std::cout << S(0) << S(1)<< S(2) << std::endl ;
S += M;
std::cout << S(0) << S(1)<< S(2)<< std::endl ;
}
{
std::cerr<< "------------- Testing += me---------------"<< std::endl;
array <double,3> A(4,2,3);
matrix<double> M(2,3); M()= 1;
A(k_,i_,j_) << i_ + 2*j_ + k_/2.0;
matrix_stack_view<double> S (std::move(A));
std::cout << S(0) << S(1)<< S(2) << std::endl ;
S += S;
std::cout << S(0) << S(1)<< S(2)<< std::endl ;
}
{
std::cerr<< "------------- Testing *= 2---------------"<< std::endl;
array <double,3> A(4,2,3);
matrix<double> M(2,3); M()= 1;
A(k_,i_,j_) << i_ + 2*j_ + k_/2.0;
matrix_stack_view<double> S (std::move(A));
std::cout << S(0) << S(1)<< S(2) << std::endl ;
S *= 2;
std::cout << S(0) << S(1)<< S(2)<< std::endl ;
}
{
std::cerr<< "------------- Testing R, L mul ---------------"<< std::endl;
array <double,3> A(4,3,3);
matrix<double> L(2,3); L()= 2;
matrix<double> R(3,2); R()= 3;
A(k_,i_,j_) << i_ + 2*j_ + k_/2.0;
matrix_stack_view<double> S (std::move(A));
std::cout << S(0) << S(1)<< S(2) << std::endl ;
auto Sb = matmul_L_R(L,S,R);
std::cout << Sb(0) << Sb(1)<< Sb(2)<< std::endl ;
std::cout << S(0) << std::endl ;
}
}

View File

@ -96,7 +96,10 @@ namespace triqs { namespace arrays {
// to forbid serialization of views...
//template<class Archive> void serialize(Archive & ar, const unsigned int version) = delete;
};
template <typename... INT> friend array_view transposed_view(array_view const& a, INT... is) {
return array_view{transpose(a.indexmap_, is...), a.storage_};
}
};
#undef IMPL_TYPE
@ -213,6 +216,13 @@ namespace triqs { namespace arrays {
TRIQS_DEFINE_COMPOUND_OPERATORS(array);
template <typename... INT> friend const_view_type transposed_view(array const& a, INT... is) {
return view_type{transpose(a.indexmap_, is...), a.storage_};
};
template <typename... INT> friend view_type transposed_view(array& a, INT... is) {
return view_type{transpose(a.indexmap_, is...), a.storage_};
};
};//array class
#undef IMPL_TYPE

View File

@ -48,10 +48,9 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid {
domain_t () =default;
domain_t (const domain_t & C) = default;
domain_t (domain_t && C) { *this = std::move(C);}
friend void swap( domain_t & a, domain_t & b) { swap(a.lengths_,b.lengths_);}
domain_t (domain_t && C) = default;
domain_t & operator =( domain_t const &) = default;
domain_t & operator =( domain_t && x) { swap(*this,x); return *this;}
domain_t & operator =( domain_t && x)= default;
domain_t (n_uple lengths):lengths_(std::move(lengths)) {}
domain_t (mini_vector<int,Rank> const & lengths):lengths_(lengths) {}

View File

@ -60,31 +60,36 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid {
typedef domain_t<Rank> domain_type;
domain_type const & domain() const { return mydomain;}
// basic construction
map (memory_layout<Rank> const & ml = memory_layout<Rank>(traversal_order)):mydomain(), start_shift_(0), memory_order_(ml) {}
map(domain_type const & C): mydomain(C), start_shift_(0), memory_order_(traversal_order) {compute_stride_compact();}
map(domain_type const & C, memory_layout<Rank> ml): mydomain(C), start_shift_(0), memory_order_(ml) {compute_stride_compact();}
/// Construction from the length, the stride, start_shift
map(lengths_type const & Lengths, strides_type const & strides, std::ptrdiff_t start_shift ):
mydomain(Lengths), strides_(strides), start_shift_(start_shift),
memory_order_ (memory_layout_from_strides(strides_)) {}
/// Construction from the length, the stride, start_shift
map(lengths_type && Lengths, strides_type && strides, std::ptrdiff_t start_shift ):
mydomain(std::move(Lengths)), strides_(std::move(strides)), start_shift_(start_shift),
memory_order_ (memory_layout_from_strides(strides_)) {}
/// Construction from another map with the same order (used in grouping indices)
template<ull_t Opt2, ull_t To2> map (map<Rank,Opt2,To2> const & C):
mydomain(C.domain()), strides_(C.strides()), start_shift_(C.start_shift()), memory_order_ (C.memory_indices_layout()) {}
// regular type
// semi-regular type
map (): start_shift_(0), memory_order_(memory_layout<Rank>(traversal_order)) {}
map (map const & C) = default;
map (map && C) = default;
map & operator = (map const & m) = default;
map & operator = (map && m) = default;
// basic construction
map(memory_layout<Rank> const & ml):mydomain(), start_shift_(0), memory_order_(ml) {}
map(domain_type const & C): mydomain(C), start_shift_(0), memory_order_(traversal_order) {compute_stride_compact();}
map(domain_type const & C, memory_layout<Rank> ml): mydomain(C), start_shift_(0), memory_order_(ml) {compute_stride_compact();}
/// Construction from the length, the stride, start_shift
map(lengths_type Lengths, strides_type strides, std::ptrdiff_t start_shift ):
mydomain(std::move(Lengths)), strides_(std::move(strides)), start_shift_(start_shift),
memory_order_ (memory_layout_from_strides(strides_)) {}
/// Construction from the length, the stride, start_shift, ml
map(lengths_type Lengths, strides_type strides, std::ptrdiff_t start_shift, memory_layout<Rank> const & ml ):
mydomain(std::move(Lengths)), strides_(std::move(strides)), start_shift_(start_shift), memory_order_ (ml) {}
/// Construction from another map with the same order (used in grouping indices)
template<ull_t Opt2, ull_t To2> map (map<Rank,Opt2,To2> const & C):
mydomain(C.domain()), strides_(C.strides()), start_shift_(C.start_shift()), memory_order_ (C.memory_indices_layout()) {}
// transposition
template <typename... INT> friend map transpose(map const& m, INT... is) {
return map{{m.domain().lengths()[is]...}, {m.strides_[is]...}, m.start_shift_, transpose(m.memory_order_, is...)};
}
/// Returns the shift in position of the element key.
template <typename KeyType>
size_t operator[] (KeyType const & key ) const {
@ -192,6 +197,9 @@ namespace triqs { namespace arrays { namespace indexmaps { namespace cuboid {
};
}; //------------- end class ---------------------
}//namespace cuboid
template<int R1, int R2, ull_t OptFlags1, ull_t OptFlags2, ull_t To1, ull_t To2>

View File

@ -65,8 +65,6 @@ namespace triqs { namespace arrays {
template< int rank, ull_t fl, ull_t to> struct get_traversal_order { static constexpr ull_t value = _get_traversal_order (rank,fl,to); };
}}
struct memory_layout_fortran {};
struct memory_layout_c {};
@ -90,18 +88,23 @@ namespace triqs { namespace arrays {
static_assert( sizeof...(in)==Rank-2, "Error");
}
memory_layout (const memory_layout & C) = default;
memory_layout (memory_layout && C) { *this = std::move(C);}
friend void swap( memory_layout & a, memory_layout & b){ std::swap(a.value,b.value);}
memory_layout (memory_layout && C) = default;
memory_layout & operator =( memory_layout const &) = default;
memory_layout & operator =( memory_layout && x) { swap(*this,x); return *this;}
memory_layout & operator =( memory_layout && x) = default;
bool operator ==( memory_layout const & ml) const { return value == ml.value;}
bool operator !=( memory_layout const & ml) const { return value != ml.value;}
friend std::ostream & operator <<( std::ostream & out, memory_layout const & s) { permutations::print(out,s.value); return out;}
friend std::ostream &operator<<(std::ostream &out, memory_layout const &s) {
permutations::print(out, s.value);
return out;
}
};
//template<int Rank> void swap( memory_layout<Rank> & a, memory_layout<Rank> & b) { std::swap(a.value,b.value);}
template <int R, typename... INT> memory_layout<R> transpose(memory_layout<R> ml, INT... is) {
static_assert(sizeof...(INT)==R, "!");
return memory_layout<R>{permutations::compose(ml.value, permutations::inverse(permutations::permutation(is...)))};
}
}}//namespace triqs::arrays
#endif

View File

@ -1,86 +0,0 @@
/*******************************************************************************
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2013 by O. Parcollet
*
* TRIQS is free software: you can redistribute it and/or modify it under the
* terms of the GNU General Public License as published by the Free Software
* Foundation, either version 3 of the License, or (at your option) any later
* version.
*
* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
* details.
*
* You should have received a copy of the GNU General Public License along with
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef TRIQS_ARRAYS_MATRIX_STACK_VIEW_H
#define TRIQS_ARRAYS_MATRIX_STACK_VIEW_H
#include "./array.hpp"
#include "./matrix.hpp"
#include "./matrix_tensor_proxy.hpp"
#include <triqs/arrays/linalg/det_and_inverse.hpp>
namespace triqs { namespace arrays {
template<typename T> class matrix_stack_view {
array_view<T,3> a;
public:
typedef array_view<T,3> array_view_t;
matrix_stack_view( typename array_view_t::view_type const & a_):a(a_) {}
matrix_stack_view( array_view_t && a_):a(std::move(a_)) {}
#ifdef TRIQS_WITH_PYTHON_SUPPORT
/// Build from a numpy.array (X is a borrowed reference) : throws if X is not a numpy.array
explicit matrix_stack_view (PyObject * X):a(typename array_view_t::view_type (X)){}
#endif
//matrix_view_proxy <array_view_t,0> operator()(size_t i) { return matrix_view_proxy <array_view_t,0>(a,i);}
//const_matrix_view_proxy<array_view_t,0> operator()(size_t i) const { return const_matrix_view_proxy <array_view_t,0>(a,i);}
auto operator()(long i) DECL_AND_RETURN(make_matrix_proxy(this->a, i));
auto operator()(long i) const DECL_AND_RETURN(make_const_matrix_proxy(this->a, i));
matrix_view<T> view(size_t i) const { return a(i,range(),range());}
size_t size() const { return a.shape(0);}
// BE CAREFUL to the shift : it is 1, 2, not 0,1, because of the stack !
friend size_t first_dim (matrix_stack_view const & m) { return second_dim(m.a);}
friend size_t second_dim (matrix_stack_view const & m) { return third_dim(m.a);}
matrix_stack_view & operator +=(matrix_stack_view const & arg) { a += arg.a; return *this; }
matrix_stack_view & operator -=(matrix_stack_view const & arg) { a -= arg.a; return *this; }
template<typename RHS>
typename std::enable_if<RHS::rank ==2, matrix_stack_view &>::type
operator +=(RHS const & arg) { for (size_t i=0; i<size(); ++i) view(i) +=arg; return *this; }
template<typename RHS>
typename std::enable_if<RHS::rank ==2, matrix_stack_view &>::type
operator -=(RHS const & arg) { for (size_t i=0; i<size(); ++i) view(i) -=arg; return *this; }
template<typename RHS> matrix_stack_view & operator *=(RHS const & arg) { a*= arg; return *this; }
template<typename RHS> matrix_stack_view & operator /=(RHS const & arg) { a/= arg; return *this; }
void invert() {for (size_t i=0; i<size(); ++i) { auto v = view(i); v = inverse(v);} }
friend matrix_stack_view matmul_L_R ( matrix_view<T> const & L, matrix_stack_view const & M, matrix_view<T> const & R) {
matrix_stack_view res (typename array_view_t::regular_type (M.size(), first_dim(L), second_dim(R)));
for (size_t i=0; i<M.size(); ++i) { res.view(i) = L * M.view(i) * R; }
return res;
}
void onsite_matmul_L_R ( matrix_view<T> const & L, matrix_stack_view const & M, matrix_view<T> const & R) {
if ((first_dim(*this) != first_dim(L)) || (second_dim(*this) != second_dim(R)) || (second_dim(L) != first_dim(R)))
TRIQS_RUNTIME_ERROR << "dimensions do not match!";
for (size_t i=0; i<M.size(); ++i) { view(i) = L * M.view(i) * R; }
}
};
}}
#endif

View File

@ -676,13 +676,12 @@ namespace gfs {
return res;
}
// ---- transpose : a new view
/*
// ---- transpose : a new gf
template <typename Variable, typename Opt>
gf_view<Variable, matrix_valued, Opt> invert_in_place(gf_view<Variable, matrix_valued, Opt> g) {
return { g.mesh(), transpose( g.data(), xxxx), transpose(g.singularity()), g.symmetry(), g.indices() XXX, g.name};
gf<Variable, matrix_valued, Opt> transpose(gf_view<Variable, matrix_valued, Opt> g) {
return {g.mesh(), transposed_view(g.data(), 0, 2, 1), transpose(g.singularity()), g.symmetry(), transpose(g.indices()), g.name};
}
*/
// ---- conjugate : always a new function -> changelog

View File

@ -281,6 +281,8 @@ namespace triqs { namespace gfs { namespace local {
inline tail conj(tail_view t) { return {conj(t.data()), t.mask_view(),t.order_min()};}
inline tail transpose(tail_view t) { return {transposed_view(t.data(),0,2,1), transposed_view(t.mask_view(),1,0),t.order_min()};}
/// Slice in orbital space
//template<bool V> tail_view slice_target(tail_impl<V> const & t, tqa::range R1, tqa::range R2) {
inline tail_view slice_target(tail_view t, tqa::range R1, tqa::range R2) {

View File

@ -142,6 +142,16 @@ namespace gfs {
TRIQS_RUNTIME_ERROR << "Not implemented : slice of string indices";
}
inline indices_2 transpose(indices_2 const &x) {
auto ind2 = x.ind;
if (x.ind.size() > 0) {
int im = x.ind.size(), jm = x.ind[0].size();
for (int i = 0; i < im; ++i)
for (int j = 0; j < jm; ++j) ind2[i][j] = x.ind[j][i];
}
return indices_2(std::move(ind2));
}
//------------------------------------------------------
// A simple replacement of tail when there is none to maintain generic code simple...
@ -162,6 +172,7 @@ namespace gfs {
template<typename A> bool check_size(A) {return true;}
};
inline nothing transpose(nothing) { return {};}
template <typename... T> nothing slice_target(nothing, T...) { return nothing(); }
template <typename T> nothing operator+(nothing, T const &) { return nothing(); }
template <typename T> nothing operator-(nothing, T const &) { return nothing(); }

View File

@ -1,4 +1,5 @@
#pragma once
#pragma once
#include <Python.h>
#include "structmember.h"
#include <string>
@ -198,11 +199,14 @@ template <> struct py_converter<double> {
template <> struct py_converter<std::complex<double>> {
static PyObject *c2py(std::complex<double> x) { return PyComplex_FromDoubles(x.real(), x.imag()); }
static std::complex<double> py2c(PyObject *ob) {
auto r = PyComplex_AsCComplex(ob);
return {r.real, r.imag};
if (PyComplex_Check(ob)) {
auto r = PyComplex_AsCComplex(ob);
return {r.real, r.imag};
}
return PyFloat_AsDouble(ob);
}
static bool is_convertible(PyObject *ob, bool raise_exception) {
if (PyComplex_Check(ob)) return true;
if (PyComplex_Check(ob) || PyFloat_Check(ob) || PyInt_Check(ob)) return true;
if (raise_exception) { PyErr_SetString(PyExc_TypeError, "Can not convert to complex");}
return false;
}