From f0dfabff3882562fb27e2bcab4df980aa3f5b9ed Mon Sep 17 00:00:00 2001 From: Michel Ferrero Date: Wed, 11 Sep 2013 18:23:38 +0200 Subject: [PATCH] Change tail implementation with fixed array size Now the tail have a fixed size. It actually makes everything simpler. I took order_min = -1 and order_max = 8. This makes the tails compatible with the previous implementation. However we might want to change this to something like -10, 10 so that they are self-contained. This commit should also fix issue #11. --- doc/changelog.rst | 11 +- pytriqs/gf/local/descriptors.py | 47 +-- pytriqs/gf/local/gf_generic.py | 19 +- pytriqs/gf/local/gf_imfreq.py | 4 +- pytriqs/gf/local/gf_imtime.py | 2 +- pytriqs/gf/local/gf_refreq.py | 2 +- pytriqs/gf/local/gf_retime.py | 2 +- pytriqs/gf/local/tail.pxd | 2 +- pytriqs/gf/local/tail.pyx | 40 +- test/pytriqs/base/gf_init.output.h5 | Bin 45824 -> 47848 bytes test/triqs/gfs/test_gf_triqs.cpp | 4 +- triqs/gfs/local/functions.cpp | 2 +- triqs/gfs/local/tail.hpp | 593 +++++++++++++--------------- 13 files changed, 359 insertions(+), 369 deletions(-) diff --git a/doc/changelog.rst b/doc/changelog.rst index 4fc2484c..1a455f22 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -6,10 +6,15 @@ Changelog This document describes the main changes in TRIQS. -From TRIQS 0.x to TRIQS 1.0 ---------------------------- +master (latest commit on github) +-------------------------------- -There have been changes from versions 0.x to 1.0 that will most likely have +* The tails now have fixed size avoid mpi problems + +version 1.0.0 +------------- + +There have been changes from versions 0.x to 1.0.0 that will most likely have consequences for your scripts and archives. Python classes diff --git a/pytriqs/gf/local/descriptors.py b/pytriqs/gf/local/descriptors.py index dd94a740..ac463a60 100644 --- a/pytriqs/gf/local/descriptors.py +++ b/pytriqs/gf/local/descriptors.py @@ -23,7 +23,7 @@ r""" """ import numpy from math import * -from gf import MeshImFreq, TailGf, MeshReFreq +from gf import MeshImFreq, MeshReFreq from lazy_expressions import LazyExprTerminal, LazyExpr, transform def is_lazy(y): @@ -68,7 +68,7 @@ class Function (Base): :param function: the function:math:`\omega \rightarrow function(\omega)` :param tail: The tail. Use None if you don't use any tail (will be put to 0) """ - Base.__init__(self,function=function, tail=tail) + Base.__init__(self, function=function, tail=tail) def __call__(self,G): if not(callable(self.function)): raise RuntimeError, "GFInitializer.Function: f must be callable" @@ -97,8 +97,7 @@ class Const(Base): C = C*numpy.identity(G.N1) if C.shape !=(G.N1,G.N2): raise RuntimeError, "Size of constant incorrect" - t = G.tail - G.tail = TailGf(shape = t.shape, size = t.size, order_min=0) + G.tail.zero() G.tail[0][:,:] = C Function(lambda om: C, None)(G) @@ -116,9 +115,12 @@ class Omega_(Base): Id = numpy.identity(G.N1) G.tail.zero() G.tail[-1][:,:] = Id + for n,om in enumerate(G.mesh): G.data[n,:,:] = om*Id return G +########################################################################## + Omega = Omega_() iOmega_n = Omega_() @@ -139,8 +141,7 @@ class A_Omega_Plus_B(Base): if A.shape !=(G.N1,G.N2): raise RuntimeError, "Size of A incorrect" if B.shape !=(G.N1,G.N2): raise RuntimeError, "Size of B incorrect" - t,Id = G.tail, numpy.identity(G.N1) - G.tail = TailGf(shape = t.shape, size = t.size, order_min=-1) + G.tail.zero() G.tail[-1][:,:] = A G.tail[0][:,:] = B @@ -160,17 +161,17 @@ class OneFermionInTime(Base): if G.mesh.TypeGF not in [GF_Type.Imaginary_Time]: raise TypeError, "This initializer is only correct in frequency" - t,Id = G.tail, numpy.identity(G.N1) - G.tail = TailGf(shape = t.shape, size = 3, order_min=1) - t[1][:,:] = 1 - t[2][:,:] = L - t[3][:,:] = L*L + Id = numpy.identity(G.N1) + G.tail.zero() + G.tail[1][:,:] = 1*Id + G.tail[2][:,:] = L*Id + G.tail[3][:,:] = L*L*Id + G.tail.mask.fill(3) fact = -1/(1+exp(-L*G.beta)) Function(lambda t: fact* exp(-L*t) *Id, None)(G) return G - ################################################## def _SemiCircularDOS(half_bandwidth): @@ -223,11 +224,12 @@ class SemiCircular (Base): raise TypeError, "This initializer is only correct in frequency" # Let's create a new tail - G.tail = TailGf(shape = G.tail.shape, size=5, order_min=1) - for i in range(G.N1): - G.tail[1][i,i] = 1.0 - G.tail[3][i,i] = D**2/4.0 - G.tail[5][i,i] = D**4/8.0 + Id = numpy.identity(G.N1) + G.tail.zero() + G.tail[1][:,:] = 1.0*Id + G.tail[3][:,:] = D**2/4.0*Id + G.tail[5][:,:] = D**4/8.0*Id + G.tail.mask.fill(6) Function(f,None)(G) return G @@ -267,11 +269,12 @@ class Wilson (Base): raise TypeError, "This initializer is only correct in frequency" # Let's create a new tail - G.tail = TailGf(shape = G.tail.shape, size=5, order_min=1) - for i in range(G.N1): - G.tail[1][i,i] = 1.0 - G.tail[3][i,i] = D**2/3.0 - G.tail[5][i,i] = D**4/5.0 + Id = numpy.identity(G.N1) + G.tail.zero() + G.tail[1][:,:] = 1.0*Id + G.tail[3][:,:] = D**2/3.0*Id + G.tail[5][:,:] = D**4/5.0*Id + G.tail.mask.fill(6) Function(f,None)(G) return G diff --git a/pytriqs/gf/local/gf_generic.py b/pytriqs/gf/local/gf_generic.py index f098df44..893fc0e1 100644 --- a/pytriqs/gf/local/gf_generic.py +++ b/pytriqs/gf/local/gf_generic.py @@ -238,15 +238,10 @@ class GfGeneric: else: raise RuntimeError, " argument type not recognized in += for %s"%arg if rhs !=None : - new_tail = TailGf(shape=lhs.tail.shape, size=lhs.tail.size, order_min=min(0,lhs.tail.order_min)) + new_tail = TailGf(shape=lhs.tail.shape) new_tail[0][:,:] = rhs - # if it is add, then we CAN change the shape of the tail, and - # reassign it since it is a new object, just create (then use the - # _singularity object. - # otherwise we can not, since it could be view, so we use the tail - # and if shape is not correct, = i.e. copy_from will raise an error - if is_add : lhs._singularity = lhs.tail + new_tail - else : lhs.tail = lhs.tail + new_tail + if is_add : lhs._singularity = lhs.tail + new_tail + else : lhs.tail = lhs.tail + new_tail return lhs def __iadd__(self, arg): @@ -275,10 +270,10 @@ class GfGeneric: else: raise RuntimeError, " argument type not recognized in -= for %s"%arg if rhs !=None : - new_tail = TailGf(shape=lhs.tail.shape, size=lhs.tail.size, order_min=min(0,lhs.tail.order_min)) + new_tail = TailGf(shape=lhs.tail.shape) new_tail[0][:,:] = rhs - if is_sub : lhs._singularity = lhs.tail - new_tail - else : lhs.tail = lhs.tail - new_tail + if is_sub : lhs._singularity = lhs.tail - new_tail + else : lhs.tail = lhs.tail - new_tail return lhs def __isub__(self, arg): @@ -336,7 +331,7 @@ class GfGeneric: MatrixStack(self.data).matmul_L_R(L, G.data, R) # this might be a bit slow - t = TailGf(shape=(N1,N2), size=G.tail.order_max-G.tail.order_min+1, order_min=G.tail.order_min) + t = TailGf(shape=(N1,N2)) for o in range(t.order_min, t.order_max+1): t[o] = numpy.dot(L, numpy.dot(G.tail[o], R)) self.tail = t diff --git a/pytriqs/gf/local/gf_imfreq.py b/pytriqs/gf/local/gf_imfreq.py index 924401c7..fded7a3b 100644 --- a/pytriqs/gf/local/gf_imfreq.py +++ b/pytriqs/gf/local/gf_imfreq.py @@ -49,7 +49,7 @@ class GfImFreq ( GfGeneric, GfImFreq_cython ) : indicesL, indicesR = indices_pack N1, N2 = len(indicesL),len(indicesR) data = d.pop('data') if 'data' in d else numpy.zeros((len(mesh),N1,N2), self.dtype ) - tail = d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2), size=10, order_min=-1) + tail = d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2)) symmetry = d.pop('symmetry', Nothing()) name = d.pop('name','g') assert len(d) ==0, "Unknown parameters in GFBloc constructions %s"%d.keys() @@ -93,7 +93,7 @@ class GfImFreq ( GfGeneric, GfImFreq_cython ) : # Change the order_max # It is assumed that any known_coef will start at order -1 - self.tail = TailGf(shape = (self.N1,self.N2), size = order_max+2, order_min = -1) + self.tail = TailGf(shape = (self.N1,self.N2)) # Fill up two arrays with the frequencies and values over the range of interest ninit, nstop = 0, -1 diff --git a/pytriqs/gf/local/gf_imtime.py b/pytriqs/gf/local/gf_imtime.py index 34ecbd8a..c87ce0fd 100644 --- a/pytriqs/gf/local/gf_imtime.py +++ b/pytriqs/gf/local/gf_imtime.py @@ -49,7 +49,7 @@ class GfImTime ( GfGeneric, GfImTime_cython ) : indicesL, indicesR = indices_pack N1, N2 = len(indicesL),len(indicesR) data = d.pop('data') if 'data' in d else numpy.zeros((len(mesh),N1,N2), self.dtype ) - tail = d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2), size=10, order_min=-1) + tail = d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2)) symmetry = d.pop('symmetry', Nothing()) name = d.pop('name','g') assert len(d) ==0, "Unknown parameters in GFBloc constructions %s"%d.keys() diff --git a/pytriqs/gf/local/gf_refreq.py b/pytriqs/gf/local/gf_refreq.py index 15cbfd77..1b5f4e67 100644 --- a/pytriqs/gf/local/gf_refreq.py +++ b/pytriqs/gf/local/gf_refreq.py @@ -47,7 +47,7 @@ class GfReFreq ( GfGeneric, GfReFreq_cython ) : indicesL, indicesR = indices_pack N1, N2 = len(indicesL),len(indicesR) data = d.pop('data') if 'data' in d else numpy.zeros((len(mesh),N1,N2), self.dtype ) - tail= d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2), size=10, order_min=-1) + tail= d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2)) symmetry = d.pop('symmetry',None) name = d.pop('name','g') assert len(d) ==0, "Unknown parameters in GFBloc constructions %s"%d.keys() diff --git a/pytriqs/gf/local/gf_retime.py b/pytriqs/gf/local/gf_retime.py index b6e5c798..65fee26c 100644 --- a/pytriqs/gf/local/gf_retime.py +++ b/pytriqs/gf/local/gf_retime.py @@ -47,7 +47,7 @@ class GfReTime ( GfGeneric, GfReTime_cython ) : indicesL, indicesR = indices_pack N1, N2 = len(indicesL),len(indicesR) data = d.pop('data') if 'data' in d else numpy.zeros((len(mesh),N1,N2), self.dtype ) - tail= d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2), size=10, order_min=-1) + tail= d.pop('tail') if 'tail' in d else TailGf(shape = (N1,N2)) symmetry = d.pop('symmetry',None) name = d.pop('name','g') assert len(d) ==0, "Unknown parameters in GFBloc constructions %s"%d.keys() diff --git a/pytriqs/gf/local/tail.pxd b/pytriqs/gf/local/tail.pxd index 24cfa97f..1143871a 100644 --- a/pytriqs/gf/local/tail.pxd +++ b/pytriqs/gf/local/tail.pxd @@ -3,7 +3,7 @@ from arrays cimport * cdef extern from "triqs/gfs/local/tail.hpp" : cdef cppclass tail "triqs::python_tools::cython_proxy" : tail() - tail(array_view[dcomplex,THREE], int, array_view[long,TWO]) except + + 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 + diff --git a/pytriqs/gf/local/tail.pyx b/pytriqs/gf/local/tail.pyx index 12effa39..c5f9e05f 100644 --- a/pytriqs/gf/local/tail.pyx +++ b/pytriqs/gf/local/tail.pyx @@ -5,8 +5,8 @@ cdef class TailGf: cdef tail _c def __init__(self, **d): """ - TailGf ( shape, size, order_min ) - TailGf ( data, order_min ) + TailGf ( shape ) + TailGf ( data, mask, order_min ) """ c_obj = d.pop('encapsulated_c_object', None) if c_obj : @@ -18,19 +18,26 @@ cdef class TailGf: if bss : assert d == {}, "Internal error : boost_serialization_string must be the only argument" boost_unserialize_into(bss,self._c) - return + return + + # default values + omin = -1 + omax = 8 - omin = d.pop('order_min') a = d.pop('data',None) if a==None : - (N1, N2), s = d.pop('shape'), d.pop('size') - a = numpy.zeros((s,N1,N2), numpy.complex) + (N1, N2) = d.pop('shape') + a = numpy.zeros((omax-omin+1,N1,N2), numpy.complex) m = d.pop('mask',None) if m==None : m = numpy.zeros(a.shape[1:3], int) - m.fill(omin+a.shape[0]-1) + m.fill(omax) + o = d.pop('order_min',None) + if o==None: o = omin + assert len(d) == 0, "Unknown parameters in TailGf constructions %s"%d.keys() - self._c = tail(array_view[dcomplex,THREE](a), omin, array_view[long,TWO](m)) + + self._c = tail(array_view[dcomplex,THREE](a), array_view[long,TWO](m), o) #-------------- Reduction ------------------------------- @@ -71,17 +78,22 @@ cdef class TailGf: def __get__(self) : return self._c.size() def copy(self) : - return self.__class__(data = self.data.copy(), order_min = self.order_min, mask = self.mask.copy()) + return self.__class__(data = self.data.copy(), mask = self.mask.copy()) def copy_from(self, TailGf T) : - assert self.order_min <= T.order_min, "Copy_from error " self._c << T._c def _make_slice(self, sl1, sl2): - return self.__class__(data = self.data[:,sl1,sl2], order_min = self.order_min, mask = self.mask[sl1,sl2]) + return self.__class__(data = self.data[:,sl1,sl2], mask = self.mask[sl1,sl2]) def __repr__ (self) : - return string.join([ "%s"%self[r]+ (" /" if r>0 else "") + " Om^%s"%(abs(r)) for r in range(self.order_min, self.order_max+1) ] , " + ") + omin = self.order_min + while ((omin <= self.order_max) and (numpy.max(numpy.abs(self.data[omin-self.order_min,:,:])) < 1e-8)): + omin = omin+1 + if omin == self.order_max+1: + return "%s"%numpy.zeros(self.shape) + 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""" @@ -164,11 +176,11 @@ cdef class TailGf: def transpose (self) : """Transpose the array : new view as in numpy""" - return TailGf(data=self.data.transpose(), order_min=self.order_min, mask=self.mask.transpose()) + return TailGf(data=self.data.transpose(), mask=self.mask.transpose()) def conjugate(self) : """Transpose the array : new view as in numpy""" - return TailGf(data=self.data.conjugate(), order_min=self.order_min, mask=self.mask) + 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) diff --git a/test/pytriqs/base/gf_init.output.h5 b/test/pytriqs/base/gf_init.output.h5 index e1f2997f04ac4723e2ffc2a66cd41f8e077080a3..5b945d4cfc8569ef98f91cef2df30c2f35f4ca5c 100644 GIT binary patch delta 4020 zcmbVPT})g>6rOwe*~M17P?g2RdZGL*$UoRFRLd3%rCO1oX-Mlsv8>DTYk^ikYq!xD zd@$L|qz!GvauwpH#&lPj)S9R~8T#hW%T|ryso5vPixC?Wow?`C-Fx>3+TH}t%$zxA z&d;20&IaztvA^V?OOAb`R=BGqA`wEI?2muB#G%s}{XZcg83wO*x3v)>(ZZ7aZG^04 z!KyQ&{VmbpZby@X`&V&4PKmD$M!tHa4M6#s{Bx^XO^rk{EJirc2$u~u^#CDOBZ78V zb_gR5*W$hzG!LhwuXm!CkH`HZzSBd!WB!RLoi*$o@8<*mNT2_tcU-rD&t0naRNgC{ zFCXNF;VDGPP_&W0r2871tz-Ey>o{$7bb28!p5;E4c{al+>J<({Iu%->0UecV6o<%2_y_Ow4ymOuObW8oIx@%~>*TgyHAK}?VrlGj(~_0_}wbVYIs zof1AKp645d7kAn1$ig-tJ*~x zs5~YG=x>#UIx@{x1yR3j=#MgRjw4>}zKHgS#KelW+4Vx$$<#^o+sY?%W-$8_-O*TJ z<3LBqoT|<1XMAWrLtjxXE9;SRn!B%Dx@4?Uh9C?RD(R8Er8yev&g8lQL>W6R@`XI$ z?gu_0L*6$g_;}3k>*oV$mhMr*Ioy3(QCur8RngaqlK?Ls zJJ0jx*v++x*y+Mwq!f<71-ELsut@CsS#$EVdb8r&b*lDV!mU$Y1SbIAIXX|(Z;sBB zFPU}rA7L{cogFWj(fK^nl7P-3UkY^Il@8GW$Cf(}$pticp;Qekw7EK$4%Rr>-wyi0 zh1y8?SLvC};qK2N@Gdpn9o3&JLW~3(J;6nk1L^ov_CRoVbO6+x8KPIuuE@^!Q@~x) z3*SD>qb#zza-Gif4;2ZB1!m^pT9@9>!@-IZx`?7~xgrr>EmzKu{P--_2j zDy-``V^*>-De5%|>wdyVC>7SpGo41Eu|pVG9G9#s2zK&ZeA$9_Kcn`0L_6oZ;HNtP z)QWa1=V8^ziGg-KxF7giquoO1=4clVn4@QZPee~gJJ*0kn6*92+2`)2Z(cp~VAhE% zK&?*27GO2Gc@5xx;FGP}1NPHHZ(ctcHL}K&k<`3;J~Mk)`zf30@EE;fhDYa85 zQov){HpuxUSU$N3G_eg5$T5d};paUztpxi$_4Bf#M^`Ri;@?iA+vJUgx!dFo5^bz* zU*2=|vMu3(=?KHo@avED5%T(?5y5zIWfBEV3@@7l7>tMjy07tdQO65kG1TJ?g1li- z$8w~12r^+2HvPyjS_JzNoO6#w6T1%E?Kgm0?^p0^SWPY(7w!lC)CQ`?EU6BXYK(!f2lT~0?m?3ojp4Drm_96#7ZYMk^q%{D=k{FarJa}1b9>Hr zzVrS4{Z8|LNh^OyYaU71NjCdMJMm7PaOSE6{TYNo7#->F*Ir4k_2l+;)vr1`#@60M2`XV~SYBe^gPK4rnXrL|_3VTE0v+yYjk5k=G$}ONH z0D~|vu=87ubh5fcOt)S@c^vqd66LiK$!{kQ`!NoV#*{e%^;Qv>YXp6MlI(w5>M`r+ zwP7|i|0-ufR}%IYL7MZsYC#Q065|*!G_Hg%c@N6r!29B3vCj9vD$5OTOEeUZoWtB< z&YWxQAU}voOW3e?8s%`{Z74H#GMeTm+e=eoD|y&};m>xjC9^i?>k93y@H>%T8vB}* zXHGe9{cvKdcQUzslYF}=YyW{ATg#FtR*~PEogJ)a{>6u{9k9_--!O@VSGKd?i?>eM z$$qo(-%(5#b@cXU^m$4Fy-#N+;YBd@re)3i+6Lk>2@a{!y*BdW5IVj;w0g*dAWLT* zrl-(44qXc^t!OYlz2O6;%w#^@lfgtVmw|PDem1Jc_R00F=abp)GIN%H8Pw2v_i0%1 zQ^T{;jGZiZp$F-%^h^8^mqUr*3U(}4^)>NCFrmg1>cn0<`Lv;(toJt7rdmOFr9g(g zE^|r^X15y4bdY3&Og`#qlD^?IXkPfZb1J7buyelxn)vy+x)jQzeWCz-Br70SW?!yh z+ruWoAP34vS|kfaLyKYsd5MegV{tVWzB&_JP!rcY-DH0Fx^%CVEH`0(Ef~VG)>mHz zL!j2%WK*p-vz*-TLn*YHx04q4qC#pfj7g(%0V58Fl)B^Yh5h_?TI0KD1f%0TlH(X2 zR|4wBQ5>aBA?;%>L_UD!BK;|a&mN$I1%eKYgM5LA1Lunnv4(>?XD+U|$eD{v<`%Ui z<*^CS%z&Q(d>T=z%t8`S4Jb1gQY-U%KLX1IuU+~CN9d_>a5`VaN2k|=5c;Ci z>t(Z!EO#M`Y(5w52+-Yru{D#mmaXdqgG2ILvo!)2#H z1FKPk2IjDsO4EQ_WnrKrt4Ys1%mNK0%>EfP5LrMe7!K2bYfh(uWsBXJ_%XZH*x0;5 z12;GghcS!NzypjVUni6XQfMEef!zCWVwrWKU>Tx!rdmBh1Ihr1ahc(1U@xxIz$fzm zrvd*pjF*1psq+=5fyE`1LaX_hP#W+o>NL>eJDvs}pzr*4N&^uzg3&+sx27fh+dl`541gHW1lHL} zFqR;_t=k~RKKc8-0srYVEFS~^gDD;U+q9fLI{atufPVOO)oIpo1pXEDoGW91{}j58 zzKB9C{L=pdCdzE7!8GnS)W8e9Ss6E(Oa>-d^tt?3Quxo^Ayd6$wF=(W4`G}+#g+)} zzYB)|!ky#VK|^o@;TAVg3Wmc7=lN1cxFL((yPL=C)_Uo#fpDXohN1|U22Fe@D8gmX zK8A3YFeYWDBY=@c?@aA@1mX6vce&Ew5iY-}Bi#O1MG#KAfzh9vJF(5LtCLwbaa (beta, Fermion, make_shape(1,1),100,full_bins, t); diff --git a/triqs/gfs/local/functions.cpp b/triqs/gfs/local/functions.cpp index f72982b6..f12b86bc 100644 --- a/triqs/gfs/local/functions.cpp +++ b/triqs/gfs/local/functions.cpp @@ -90,7 +90,7 @@ namespace triqs { namespace gfs { local::tail_view get_tail(gf_view const & gl, int size = 10, int omin = -1) { auto sh = gl.data().shape().front_pop(); - local::tail t(sh, size, omin); + local::tail t(sh); t.data() = 0.0; for (int p=1; p<=t.order_max(); p++) diff --git a/triqs/gfs/local/tail.hpp b/triqs/gfs/local/tail.hpp index 8cdbebcd..730c6cc8 100644 --- a/triqs/gfs/local/tail.hpp +++ b/triqs/gfs/local/tail.hpp @@ -27,375 +27,352 @@ namespace triqs { namespace gfs { namespace local { - namespace details { - static constexpr double small = 1.e-10; - } + namespace details { + static constexpr double small = 1.e-10; + } - namespace tqa= triqs::arrays; namespace tql= triqs::clef; namespace mpl= boost::mpl; - typedef std::complex dcomplex; + namespace tqa= triqs::arrays; namespace tql= triqs::clef; namespace mpl= boost::mpl; + typedef std::complex dcomplex; - class tail; // the value class - class tail_view; // the view class + class tail; // the value class + class tail_view; // the view class - 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_{}; + 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 > {}; + // 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 > {}; - // ---------------------- implementation -------------------------------- + // ---------------------- implementation -------------------------------- - /// A common implementation class. Idiom : ValueView - template class tail_impl { - public : - typedef tail_view view_type; - typedef tail regular_type; + /// A common implementation class. Idiom : ValueView + template class tail_impl { + public : + typedef tail_view view_type; + typedef tail regular_type; - typedef arrays::array data_regular_type; - typedef arrays::array_view data_view_type; - typedef typename mpl::if_c::type data_type; + typedef arrays::array data_regular_type; + typedef arrays::array_view data_view_type; + typedef typename mpl::if_c::type data_type; - typedef arrays::array mask_regular_type; - typedef arrays::array_view mask_view_type; - typedef typename mpl::if_c::type mask_type; + typedef arrays::array mask_regular_type; + typedef arrays::array_view mask_view_type; + typedef typename mpl::if_c::type mask_type; - typedef arrays::matrix_view mv_type; - typedef arrays::matrix_view const_mv_type; - //typedef arrays::matrix_view const_mv_type; + typedef arrays::matrix_view mv_type; + typedef arrays::matrix_view const_mv_type; - data_view_type data() { return _data;} - const data_view_type data() const { return _data;} - mask_view_type mask_view() { return mask;} - const mask_view_type mask_view() const { return mask;} + data_view_type data() { return _data;} + const data_view_type data() const { return _data;} + mask_view_type mask_view() { return mask;} + const mask_view_type mask_view() const { return mask;} - long order_min() const {return omin;} - long order_max() const {return min_element(mask);} - size_t size() const {return _data.shape()[0];} - long smallest_nonzero() const { - long om = omin; - while ((om < this->order_max()) && (max_element(abs(_data(om-omin,tqa::range(),tqa::range()))) < details::small)) om++; - return om; - } + long order_min() const {return omin;} + long order_max() const {return min_element(mask);} + size_t size() const {return _data.shape()[0];} + long smallest_nonzero() const { + long om = omin; + while ((om < this->order_max()) && (max_element(abs(_data(om-omin,tqa::range(),tqa::range()))) < details::small)) om++; + return om; + } - 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];} + 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);} + bool is_decreasing_at_infinity() const { return (smallest_nonzero() >=1);} - protected: + protected: - long omin; - mask_type mask; - data_type _data; + long omin; + mask_type mask; + data_type _data; - // All constructors - tail_impl(): omin(0), mask(), _data() {} // all arrays of zero size (empty) - tail_impl(size_t N1, size_t N2, size_t size_, long order_min): - omin(order_min), mask(tqa::make_shape(N1,N2)), _data(tqa::make_shape(size_,N1,N2)) { - mask() = order_min+size_-1; - _data() = 0; - } - tail_impl(data_type const &d, long order_min, mask_type const &om): omin(order_min), mask(om), _data(d) {} - // tail_impl(tail_impl const & x): omin(x.omin), mask(x.mask), _data(x._data){} - tail_impl(tail_impl const & x): omin(x.omin), mask(x.mask), _data(x._data){} - tail_impl(tail_impl const &) = default; - tail_impl(tail_impl &&) = default; - - friend class tail_impl; - public: + // All constructors + tail_impl(): mask(), _data(), omin(0) {} // all arrays of zero size (empty) + tail_impl(size_t N1, size_t N2, long omin_, long size_): + omin(omin_), mask(tqa::make_shape(N1,N2)), _data(tqa::make_shape(size_,N1,N2)) { + mask() = omin+size_-1; + _data() = 0; + } + tail_impl(data_type const &d, mask_type const &m, long omin_): mask(m), _data(d), omin(omin_) {} + tail_impl(tail_impl const & x): mask(x.mask), _data(x._data), omin(x.omin) {} + tail_impl(tail_impl const &) = default; + tail_impl(tail_impl &&) = default; - mv_type operator() (int n) { - if (n>this->order_max()) TRIQS_RUNTIME_ERROR<<" n > Max Order. n= "< Max Order. n= "< Max Order. n= "< Max Order. n= "< prefer serialization or hdf5 ! + void save(std::string file, bool accumulate=false) const {} - friend void h5_read (h5::group fg, std::string subgroup_name, tail_impl & t){ - auto gr = fg.open_group(subgroup_name); - // Check the attribute or throw - //auto tag_file = gr.read_triqs_hdf5_data_scheme(); - //auto tag_expected= get_triqs_hdf5_data_scheme(t); - //if (tag_file != tag_expected) - // TRIQS_RUNTIME_ERROR<< "h5_read : mismatch of the tag TRIQS_HDF5_data_scheme tag in the h5 group : found "< - void serialize(Archive & ar, const unsigned int version) { - ar & boost::serialization::make_nvp("omin",omin); - ar & boost::serialization::make_nvp("mask",mask); - ar & boost::serialization::make_nvp("data",_data); + friend std::string get_triqs_hdf5_data_scheme(tail_impl const & g) { return "TailGf"; } + + /// + friend void h5_write (h5::group fg, std::string subgroup_name, tail_impl const & t) { + auto gr = fg.create_group(subgroup_name); + // tagging the hdf5 file + //gr.write_triqs_hdf5_data_scheme(t); + h5_write(gr,"omin",t.omin); + h5_write(gr,"mask",t.mask); + h5_write(gr,"data",t._data); + } + + friend void h5_read (h5::group fg, std::string subgroup_name, tail_impl & t){ + auto gr = fg.open_group(subgroup_name); + // Check the attribute or throw + //auto tag_file = gr.read_triqs_hdf5_data_scheme(); + //auto tag_expected= get_triqs_hdf5_data_scheme(t); + //if (tag_file != tag_expected) + // TRIQS_RUNTIME_ERROR<< "h5_read : mismatch of the tag TRIQS_HDF5_data_scheme tag in the h5 group : found "< + void serialize(Archive & ar, const unsigned int version) { + ar & boost::serialization::make_nvp("omin",omin); + ar & boost::serialization::make_nvp("mask",mask); + ar & boost::serialization::make_nvp("data",_data); + } + + friend std::ostream & operator << (std::ostream & out, tail_impl const & x) { + out <<"tail/tail_view: min/smallest/max = "<< x.order_min() << " " << x.smallest_nonzero() << " "<< x.order_max(); + for (long u = x.order_min(); u <= x.order_max(); ++u) out <<"\n ... Order "< const & x) { + _data() = 0.0; + mv_type(_data(-omin, tqa::range(), tqa::range())) = x; + mask() = omin+_data.shape()[0]-1; + return *this; } - friend std::ostream & operator << (std::ostream & out, tail_impl const & x) { - out <<"tail/tail_view: min/smallest/max = "<< x.order_min() << " " << x.smallest_nonzero() << " "<< x.order_max(); - for (long u = x.order_min(); u <= x.order_max(); ++u) out <<"\n ... Order "< { - typedef tail_impl B; - friend class tail; + // ----------------------------- + // the regular class + class tail : public tail_impl { + typedef tail_impl B; + friend class tail_view; + public : + tail():B() {} + typedef tqa::mini_vector shape_type; + tail(size_t N1, size_t N2, long order_min=-1, long size=10): B(N1,N2,order_min,size) {} + tail(shape_type const & sh, long order_min=-1, long size=10): B(sh[0],sh[1],order_min,size) {} + tail(tail const & g): B(g) {} + tail(tail_view const & g): B(g) {} + tail(tail &&) = default; - public : - template tail_view(tail_impl const & t): B(t){} - tail_view(B::data_type const &d, long order_min, B::mask_type const &om): B(d, order_min, om){} - tail_view(tail_view const &) = default; - tail_view(tail_view &&) = default; - void rebind( tail_view const &X) { + // operator = for values + tail & operator = (tail_view const & rhs) { + omin = rhs.omin; + mask = rhs.mask; + _data = rhs._data; + return *this; + } + tail & operator = (tail const & rhs) { + omin = rhs.omin; + mask = rhs.mask; + _data = rhs._data; + return *this; + } + + using B::operator(); + + /// The simplest tail corresponding to : omega + static tail_view omega(size_t N1, size_t N2) { + tail t(N1, N2); + t(-1) = 1; + return t; + } + + /// The simplest tail corresponding to : omega, constructed from a shape for convenience + static tail_view omega(tail::shape_type const & sh) { return omega(sh[0],sh[1]); } + + }; + + template void assign_from_expression(tail_view & t,RHS const & rhs) { t = rhs( tail::omega(t.shape()) ); } + + inline void tail_view::rebind(tail const &X) { omin = X.omin; mask.rebind(X.mask); _data.rebind(X._data); } - inline void rebind( tail const &X); - //using B::operator=; // import operator = from impl. class or the default = is synthetized - //tail_view & operator=(const tail_view & rhs) { if (this->data.is_empty()) rebind(rhs); else B::operator=(rhs); return *this; } - - // operator = for views - tail_view & operator = (const tail_view & rhs) { - if (this->_data.is_empty()) rebind(rhs); - else { - if (rhs.omin < omin) TRIQS_RUNTIME_ERROR<<"rhs has too small omin"; - if ((_data.shape()[1] != rhs._data.shape()[1]) || (_data.shape()[2] != rhs._data.shape()[2])) - TRIQS_RUNTIME_ERROR<<"rhs has different shape"; - for (size_t i=0; i const & x) { - if (omin > 0) TRIQS_RUNTIME_ERROR<<"lhs has too large omin"; - for (size_t n=0; n { - typedef tail_impl B; - friend class tail_view; - public : - tail():B() {} - 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(tail const & g): B(g){} - tail(tail_view const & g): B(g){} - tail(tail &&) = default; - - // operator = for values - tail & operator = (tail_view const & rhs) { - omin = rhs.omin; - mask = rhs.mask; - _data = rhs._data; - return *this; - } - tail & operator = (tail const & rhs) { - omin = rhs.omin; + inline tail_view & tail_view::operator = (const tail & rhs) { + if ((_data.shape()[1] != rhs._data.shape()[1]) || (_data.shape()[2] != rhs._data.shape()[2]) || (omin != rhs.omin)) + TRIQS_RUNTIME_ERROR<<"rhs has different shape"; mask = rhs.mask; _data = rhs._data; return *this; } - using B::operator(); - - /// The simplest tail corresponding to : omega - static tail_view omega(size_t N1, size_t N2, size_t size_) { - tail t(N1, N2, size_, -1); - t(-1) = 1; - return t; + /// Slice in orbital space + template tail_view slice_target(tail_impl const & t, tqa::range R1, tqa::range R2) { + return tail_view(t.data()(tqa::range(),R1,R2), t.mask_view()(R1,R2), t.order_min()); } - /// The simplest tail corresponding to : omega, constructed from a shape for convenience - static tail_view omega(tail::shape_type const & sh, size_t size_) { return omega(sh[0],sh[1],size_);} + inline tail inverse(tail_view const & t) { + long omin1 = - t.smallest_nonzero(); + long omax1 = std::min(t.order_max() + 2*omin1, t.order_min()+long(t.size())-1); + size_t si = omax1-omin1+1; - }; + tail res(t); + res.data() = 0.0; + res.mask_view() = omax1; - template void assign_from_expression(tail_view & t,RHS const & rhs) { t = rhs( tail::omega(t.shape(),t.size())); } - - inline void tail_view::rebind( tail const &X) { - omin = X.omin; - mask.rebind(X.mask); - _data.rebind(X._data); - } - inline tail_view & tail_view::operator = (const tail & rhs) { - if (this->_data.is_empty()) rebind(rhs); - else { - if (rhs.omin < omin) TRIQS_RUNTIME_ERROR<<"rhs has too small omin"; - if ((_data.shape()[1] != rhs._data.shape()[1]) || (_data.shape()[2] != rhs._data.shape()[2])) - TRIQS_RUNTIME_ERROR<<"rhs has different shape"; - for (size_t i=0; i tail_view slice_target(tail_impl const & t, tqa::range R1, tqa::range R2) { - return tail_view(t.data()(tqa::range(),R1,R2),t.order_min(),t.mask_view()(R1,R2)); - } + inline tail mult_impl(tail_view const & l, tail_view const& r) { + if (l.shape()[1] != r.shape()[0] || l.order_min() != r.order_min()) TRIQS_RUNTIME_ERROR<< "tail multiplication : shape mismatch"; + long omin1 = l.smallest_nonzero() + r.smallest_nonzero(); + long omax1 = std::min(r.order_max()+l.smallest_nonzero(), l.order_max()+r.smallest_nonzero()); + size_t si = omax1-omin1+1; - inline tail inverse(tail_view const & t) { + tail res(l); + res.data() = 0.0; + res.mask_view() = omax1; - // find in t - long omin1 = - t.smallest_nonzero(); - long omax1 = t.order_max() + 2*omin1; - size_t si = omax1-omin1+1; - tail t_inv(t.shape(), si, omin1); - - t_inv(omin1) = inverse(t(-omin1)); - - for (size_t n=1; n= a.n_min and n-p <=b.n_max and n-p >= b.n_min - // hence p <= min ( a.n_max, n-b.n_min ) and p >= max ( a.n_min, n- b.n_max) - const long pmin = std::max(l.smallest_nonzero(), n - r.order_max() ); - const long pmax = std::min(l.order_max(), n - r.smallest_nonzero() ); - for (long p = pmin; p <= pmax; ++p) { res(n) += l(p) * r(n-p);} - } - 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) { - tail res(b); for (long n=res.order_min(); n<=res.order_max(); ++n) res(n)=a*res(n); return res; + for (long n=res.order_min(); n<=res.order_max(); ++n) { + // sum_{p}^n a_p b_{n-p}. p <= a.n_max, p >= a.n_min and n-p <=b.n_max and n-p >= b.n_min + // hence p <= min ( a.n_max, n-b.n_min ) and p >= max ( a.n_min, n- b.n_max) + const long pmin = std::max(l.smallest_nonzero(), n - r.order_max() ); + const long pmax = std::min(l.order_max(), n - r.smallest_nonzero() ); + for (long p = pmin; p <= pmax; ++p) { res(n) += l(p) * r(n-p);} + } + 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_, LocalTail>) + operator* (T1 const & a, T2 const & b) { return mult_impl(a,b); } - 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; } + 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) { return a *inverse(b); } + 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; + } - inline tail operator / (tail_view const & r, dcomplex a) { tail res(r); res.data() /=a; return res;} - inline tail operator / (dcomplex a, tail_view const & r) { return a * inverse(r); } + 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; } - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator + (T1 const & l, T2 const& r) { - using arrays::range; - if (l.shape() != r.shape()) TRIQS_RUNTIME_ERROR<< "tail addition : shape mismatch"; - long omin1 = std::min(l.smallest_nonzero(),r.smallest_nonzero()); - long omax1 = std::min(l.order_max(),r.order_max()); - size_t si = omax1-omin1+1; - tail res(l.shape(), si, omin1); - for (long i = res.order_min(); i<=res.order_max(); ++i) res(i) = l(i) + r(i); - return res; - } + template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) + operator/ (T1 const & a, T2 const & b) { return a * inverse(b); } - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator - (T1 const & l, T2 const& r) { - using arrays::range; - if (l.shape() != r.shape()) TRIQS_RUNTIME_ERROR<< "tail substraction : shape mismatch"; - long omin1 = std::min(l.smallest_nonzero(),r.smallest_nonzero()); - long omax1 = std::min(l.order_max(),r.order_max()); - size_t si = omax1-omin1+1; - tail res(l.shape(), si, omin1); - for (long i = res.order_min(); i<=res.order_max(); ++i) res(i) = l(i) - r(i); - return res; - } + inline tail operator / (tail_view const & r, dcomplex a) { tail res(r); res.data() /=a; return res;} + inline tail operator / (dcomplex a, tail_view const & r) { return a * inverse(r); } - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator + (T1 const & a, T2 const & t) { - tail res(t); - res(0) += a; - return res; - } + template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) + operator + (T1 const & l, T2 const& r) { + if (l.shape() != r.shape() || l.order_min() != r.order_min()) TRIQS_RUNTIME_ERROR<< "tail addition : shape mismatch"; + tail res(l.shape()); + res.mask_view() = std::min(l.order_max(), r.order_max()); + for (long i = res.order_min(); i<=res.order_max(); ++i) res(i) = l(i) + r(i); + return res; + } + template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) + operator - (T1 const & l, T2 const& r) { + if (l.shape() != r.shape() || l.order_min() != r.order_min()) TRIQS_RUNTIME_ERROR<< "tail addition : shape mismatch"; + tail res(l.shape()); + res.mask_view() = std::min(l.order_max(), r.order_max()); + for (long i = res.order_min(); i<=res.order_max(); ++i) res(i) = l(i) - r(i); + return res; + } - template TYPE_ENABLE_IF(tail,mpl::and_, is_scalar_or_element>) - operator + (T1 const & t, T2 const & a) { return a+t;} + template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) + operator + (T1 const & a, T2 const & t) { + tail res(t); + res(0) += a; + return res; + } - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator - (T1 const & a, T2 const & t) { return (-a) + t;} + template TYPE_ENABLE_IF(tail,mpl::and_, is_scalar_or_element>) + operator + (T1 const & t, T2 const & a) { return a+t;} - template TYPE_ENABLE_IF(tail,mpl::and_, is_scalar_or_element>) - operator - (T1 const & t, T2 const & a) { return (-a) + t;} + template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) + operator - (T1 const & a, T2 const & t) { return (-a) + t;} + + template TYPE_ENABLE_IF(tail,mpl::and_, is_scalar_or_element>) + operator - (T1 const & t, T2 const & a) { return (-a) + t;} }}} #endif