#if defined HAVE_CONFIG_H
#include "config.h"
#endif
!{\src2tex{textfont=tt}}
!!****m* ABINIT/m_BathOperator
!! NAME
!!  m_BathOperator
!! 
!! FUNCTION 
!!  Manage all stuff related to the bath for the 
!!  simgle Anderson Impurity Model
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

#include "defs.h"
MODULE m_BathOperator
USE m_MatrixHyb
USE m_Vector
USE m_VectorInt
USE m_Global
USE m_ListCdagC

IMPLICIT NONE

!!***

PRIVATE

!!****t* m_BathOperator/BathOperator
!! NAME
!!  BathOperator
!!
!! FUNCTION
!!  This structured datatype contains the necessary data
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! SOURCE

TYPE, PUBLIC :: BathOperator
  LOGICAL _PRIVATE :: set         = .FALSE.
  LOGICAL          :: MAddFlag    = .FALSE. ! Set to true if we can compute a new M (see updateDetXX)
  LOGICAL          :: MRemoveFlag = .FALSE. ! Set to true if we can compute a new M (see updateDetXX)
  LOGICAL _PRIVATE :: antiShift   = .FALSE. ! shift when M is updated with antiseg
  LOGICAL _PRIVATE :: doCheck     = .FALSE.
  INTEGER _PRIVATE :: flavors
  INTEGER          :: activeFlavor
  INTEGER _PRIVATE :: samples
  INTEGER _PRIVATE :: sizeHybrid
  INTEGER _PRIVATE :: updatePosRow
  INTEGER _PRIVATE :: updatePosCol
  INTEGER _PRIVATE :: iTech
  INTEGER _PRIVATE :: checkNumber
  DOUBLE PRECISION _PRIVATE                   :: beta
  DOUBLE PRECISION _PRIVATE                   :: dt
  DOUBLE PRECISION _PRIVATE                   :: inv_dt
  DOUBLE PRECISION _PRIVATE                   :: meanError
  DOUBLE PRECISION _PRIVATE                   :: S
  DOUBLE PRECISION _PRIVATE                   :: Stau
  DOUBLE PRECISION _PRIVATE                   :: Stilde
  TYPE(Vector)     _PRIVATE                   :: R 
  TYPE(Vector)     _PRIVATE                   :: Q 
  TYPE(Vector)     _PRIVATE                   :: Rtau
  TYPE(Vector)     _PRIVATE                   :: Qtau
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) _PRIVATE :: F ! sample,Flavors
  TYPE(MatrixHyb) , ALLOCATABLE, DIMENSION(:)            :: M  ! Flavors
  TYPE(MatrixHyb) , ALLOCATABLE, DIMENSION(:)   _PRIVATE :: M_update  ! Flavors
END TYPE BathOperator
!!***

PUBLIC  :: BathOperator_init
PUBLIC  :: BathOperator_reset
PUBLIC  :: BathOperator_activateParticle
PRIVATE :: BathOperator_hybrid
PUBLIC  :: BathOperator_getDetAdd
PUBLIC  :: BathOperator_getDetRemove
PUBLIC  :: BathOperator_getDetF
PUBLIC  :: BathOperator_setMAdd
PUBLIC  :: BathOperator_setMRemove
PUBLIC  :: BathOperator_swap
PUBLIC  :: BathOperator_initF
PUBLIC  :: BathOperator_setF
PUBLIC  :: BathOperator_printF
PUBLIC  :: BathOperator_printM
PUBLIC  :: BathOperator_destroy
PUBLIC  :: BathOperator_doCheck
PRIVATE :: BathOperator_checkM
PUBLIC  :: BathOperator_getError

CONTAINS
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_init
!! NAME
!!  BathOperator_init
!!
!! FUNCTION
!!  Initialize and allocate data
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath object
!!  flavors=numbers of flavors we have (including spin)
!!  samples=Time slices in the input file
!!  beta=inverse temperature
!!  iTech=imaginary time or frequencies
!!  It is imposes to imaginary time
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_init(op, flavors, samples, beta, iTech)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_init'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
  INTEGER           , INTENT(IN   ) :: flavors
  INTEGER           , INTENT(IN   ) :: samples
  DOUBLE PRECISION  , INTENT(IN   ) :: beta
!Local variables ------------------------------
  INTEGER           , INTENT(IN   ) :: iTech
  INTEGER                           :: it

  op%MAddFlag     = .FALSE.
  op%MRemoveFlag  = .FALSE.
  op%flavors      = flavors
  op%beta         = beta
  op%samples      = samples
  op%sizeHybrid   = samples + 1
  op%dt      = beta / DBLE(samples)
  op%inv_dt  = DBLE(samples) / beta
  op%activeFlavor= 0 
  op%updatePosRow = 0
  op%updatePosCol = 0
  op%iTech        = iTech
!#ifdef CTQMC_CHECK
  op%checkNumber  = 0
  op%meanError    = 0.d0
  op%doCheck = .FALSE.
!#endif

  FREEIF(op%F)
  MALLOC(op%F,(1:op%sizeHybrid+1,1:flavors))
  DT_FREEIF(op%M)
  DT_MALLOC(op%M,(1:flavors))
  DT_FREEIF(op%M_update)
  DT_MALLOC(op%M_update,(1:flavors))
  
  CALL Vector_init(op%R,100)
  CALL Vector_init(op%Q,100)
  CALL Vector_init(op%Rtau,100)
  CALL Vector_init(op%Qtau,100)

  DO it = 1, flavors
    CALL MatrixHyb_init(op%M(it),op%iTech,size=Global_SIZE,Wmax=samples) !FIXME Should be consistent with ListCagC
    CALL MatrixHyb_init(op%M_update(it),op%iTech,size=Global_SIZE,Wmax=samples) !FIXME Should be consistent with ListCagC
  END DO
  op%F       = 0.d0
  op%set     = .TRUE.
  
END SUBROUTINE BathOperator_init
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_reset
!! NAME
!!  BathOperator_reset
!!
!! FUNCTION
!!  Reset all internal variables
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator to reset
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_reset(op)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_reset'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER                           :: it
  op%MAddFlag     = .FALSE.
  op%MRemoveFlag  = .FALSE.
  op%activeFlavor = 0 
  op%updatePosRow = 0
  op%updatePosCol = 0
!#ifdef CTQMC_CHECK
  op%checkNumber  = 0
  op%meanError    = 0.d0
!#endif
  op%doCheck = .FALSE.
  CALL Vector_clear(op%R)
  CALL Vector_clear(op%Q)
  CALL Vector_clear(op%Rtau)
  CALL Vector_clear(op%Qtau)

  DO it = 1, op%flavors
    CALL MatrixHyb_clear(op%M(it)) !FIXME Should be consistent with ListCagC
  END DO
  op%F       = 0.d0

END SUBROUTINE BathOperator_reset
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_activateParticle
!! NAME
!!  BathOperator_activateParticle
!!
!! FUNCTION
!!  Just save on wicht flavor we are working
!!  It is better to use the macro defined in defs.h
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  flavor=the flavor to activate
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_activateParticle(op,flavor)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_activateParticle'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
!Local variables ------------------------------
  INTEGER           , INTENT(IN   ) :: flavor

  IF ( flavor .GT. op%flavors ) &
    CALL ERROR("BathOperator_activateParticle : out of range      ")
  IF ( op%set .EQV. .TRUE. .AND. ALLOCATED(op%M) ) THEN 
    op%activeFlavor =  flavor
    op%MAddFlag     = .FALSE.
    op%MRemoveFlag  = .FALSE.
  ELSE
    CALL ERROR("BathOperator_activateParticle : not allocated      ")
  END IF
END SUBROUTINE BathOperator_activateParticle
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_hybrid
!! NAME
!!  BathOperator_hybrid
!!
!! FUNCTION
!!  Compute the hybridization for the active flavor
!!  at time time
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  time=time  F(time)
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION BathOperator_hybrid(op,time)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_hybrid'
!End of the abilint section

  TYPE(BathOperator), INTENT(IN) :: op
  DOUBLE PRECISION  , INTENT(IN) :: time
#include "BathOperator_hybrid.h"

  IF ( op%activeFlavor .LE. 0 ) &
    CALL ERROR("BathOperator_hybrid : no active hybrid func        ")
#include "BathOperator_hybrid"
  BathOperator_hybrid = hybrid

END FUNCTION BathOperator_hybrid
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_getDetAdd
!! NAME
!!  BathOperator_getDetAdd
!!
!! FUNCTION
!!  Compute the determinant ratio when a (anti)segment
!!  is trying to be added and store some array for setMadd
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  CdagC_1=segment to be added
!!  position=ordered position of the Cdag time
!!  particle=full list of CdagC for activeFlavor
!!
!! OUTPUT
!!  BathOperator_getDetAdd=the det 
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE
DOUBLE PRECISION  FUNCTION BathOperator_getDetAdd(op,CdagC_1, position, particle)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_getDetAdd'
!End of the abilint section

  TYPE(BathOperator)      , INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(1:2), INTENT(IN   ) :: CdagC_1
  INTEGER                 , INTENT(IN   ) :: position  
  TYPE(ListCdagC), INTENT(IN   ) :: particle   
!Local variables-------------------------------
  INTEGER                                 :: it1
  INTEGER                                 :: it2
  INTEGER                                 :: it3
  INTEGER                                 :: tail
  INTEGER                                 :: new_tail
  DOUBLE PRECISION                        :: C
  DOUBLE PRECISION                        :: Cbeta
  DOUBLE PRECISION                        :: Cibeta
  DOUBLE PRECISION                        :: Cdag
  DOUBLE PRECISION                        :: Cdagbeta
  DOUBLE PRECISION                        :: beta
  DOUBLE PRECISION                        :: ratio
  DOUBLE PRECISION                        :: time
!  TYPE(CdagC)    , POINTER, DIMENSION(:)  :: list => NULL()
#include "BathOperator_hybrid.h"

  op%antiShift = .FALSE.
  beta     = op%beta
  C        =  CdagC_1(C_)
!  Cbeta    = C.MOD.beta
  MODCYCLE(C,beta,Cbeta)
  Cdag     =  CdagC_1(Cdag_)
!  cdagbeta = Cdag.MOD.beta
  MODCYCLE(Cdag,beta,CdagBeta)
!  IF ( Cdag .GE. beta ) &
!    CALL ERROR("BathOperator_getDetAdd : bad case ...              ")
  IF ( op%activeFlavor .LE. 0 ) &
    CALL ERROR("BathOperator_getDetAdd : no active hybrid function ")

  tail =  particle%tail
  new_tail = tail+1
!  list => particle%list
  
  IF ( ((C .GT. Cdag) .AND. (position .EQ. -1)) &
       .OR. ((C .LT. Cdag) .AND. (tail .EQ. 0))) THEN ! Possible only if it is a segment
    op%updatePosRow = tail + 1
    op%updatePosCol = tail + 1
  ELSE
    op%updatePosRow  = ABS(position)
    op%updatePosCol  = ABS(position)
  END IF
  
  ! If antisegment, the det ratio has to be by -1 ( sign of the signature of one
  ! permutation line in the matrix
  IF ( C .LT. Cdag .AND. tail .GT. 0) THEN ! if antiseg
  !  ratio = -ratio 
    op%updatePosRow  = (op%updatePosRow + 1) !position in [1;tail]
    IF ( CdagBeta .LT. particle%list(op%updatePosCol,Cdag_) ) op%antiShift = .TRUE.
  END IF

!  CALL Vector_setSize(op%R,tail)
!  CALL Vector_setSize(op%Q,tail)
  Vector_QuickResize(op%R,new_tail)
  Vector_QuickResize(op%Q,new_tail)
  Vector_QuickResize(op%Rtau,new_tail)
  Vector_QuickResize(op%Qtau,new_tail)

  DO it1 = 1, tail
    it2 = it1 + ( 1+SIGN(1,it1-op%updatePosRow) )/2
    it3 = it1 + ( 1+SIGN(1,it1-op%updatePoscol) )/2

    op%Rtau%vec(it2)= C - particle%list(it1,Cdag_)
    !op%Rtau%vec(it1)= C - particle%list(it1,Cdag_)
    time = Cbeta - particle%list(it1,Cdag_)
#include "BathOperator_hybrid"
    op%R%vec(it1) = hybrid
!    op%R%vec(it) = BathOperator_hybrid(op, Cbeta - list(it)%Cdag)
!    Cibeta = list(it)%C.MOD.beta
    MODCYCLE(particle%list(it1,C_),beta,Cibeta)
    time = Cibeta - Cdagbeta
    op%Qtau%vec(it3)= time
    !op%Qtau%vec(it1)= time
#include "BathOperator_hybrid"
    op%Q%vec(it1) = hybrid
    !op%Q%vec(it3) = hybrid
!    Q(it) = BathOperator_hybrid(op, Cibeta - Cdagbeta)
  END DO
  ! Compute S
  op%Stau = C - Cdagbeta 
  op%Rtau%vec(op%updatePosRow) = op%Stau
  op%Qtau%vec(op%updatePosCol) = op%Rtau%vec(op%updatePosRow)

  time = Cbeta-Cdagbeta
#include "BathOperator_hybrid"
  op%S = hybrid

  !ratio = op%S - DOT_PRODUCT(MATMUL(op%R%vec(1:tail),op%M(op%activeFlavor)%mat(1:tail,1:tail)),op%Q%vec(1:tail))
  ratio = 0.d0
  DO it1 = 1, tail
    time = 0.d0
    DO it2 = 1, tail
      time = time + op%R%vec(it2) * op%M(op%activeFlavor)%mat(it2,it1)
    END DO
    ratio = ratio + op%Q%vec(it1) * time
  END DO
  ratio = op%S - ratio

  op%Stilde = 1.d0 / ratio

  ! This IF is the LAST "NON CORRECTION" in my opinion this should not appears.
!  IF ( MAX(C,Cdag) .GT. op%beta ) THEN
!    WRITE(*,*) op%Stilde
!    op%Stilde = - ABS(op%Stilde)
!  END IF
  BathOperator_getDetAdd = ratio
  op%MAddFlag   = .TRUE.
!#ifdef CTQMC_CHECK
!  op%ListCdagC = particle
!!write(*,*) op%Stilde
!!write(*,*) op%antishift
!!write(*,*)    op%updatePosRow 
!!write(*,*)    op%updatePosCol 
!#endif

END FUNCTION BathOperator_getDetAdd
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_getDetRemove
!! NAME
!!  BathOperator_getDetRemove
!!
!! FUNCTION
!!  Compute the determinant ratio when a (anti)segment
!!  is trying to be removed 
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  position=position of segment to be removed
!!
!! OUTPUT
!!  BathOperator_getDetRemove=the det 
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION BathOperator_getDetRemove(op,position)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_getDetRemove'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
!Local arguments-------------------------------
  INTEGER           , INTENT(IN   ) :: position  
  INTEGER                           :: ABSposition  
  INTEGER                           :: tail

  IF ( op%activeFlavor .LE. 0 ) &
    CALL ERROR("BathOperator_getDetRemove : no active hybrid fun  ")

  op%antiShift = .FALSE.
  tail         = op%M(op%activeFlavor)%tail
  ABSposition  = ABS(position)
  IF ( ABSposition .GT. tail ) &
    CALL ERROR("BathOperator_getDetRemove : position > M size     ")
  op%updatePosCol = ABSposition
  op%antiShift    = .FALSE.
  IF ( position .GT. 0 ) THEN
    op%updatePosRow = ABSposition
  ELSE
    op%updatePosRow = ABSposition+1
    IF ( ABSposition .EQ. tail ) THEN 
      op%antiShift = .TRUE.
      op%updatePosRow = 1 !ABSposition - 1
!      op%updatePosRow = ABSposition    
!      IF ( op%updatePosCol .EQ. 0) op%updatePosCol = tail
    END IF
  ENDIF
  op%Stilde                 = op%M(op%activeflavor)%mat(op%updatePosRow,op%updatePosCol) 
  op%MRemoveFlag            = .TRUE.
  BathOperator_getDetRemove = op%Stilde
!#ifdef CTQMC_CHECK
!  op%ListCdagC = particle
!!write(*,*) op%updatePosRow, op%updatePosCol, position
!!CALL ListCdagC_print(particle)
!#endif

END FUNCTION BathOperator_getDetRemove
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_getDetF
!! NAME
!!  BathOperator_getDetF
!!
!! FUNCTION
!!  Compute the determinant of the F matrix
!!  using the hybridization of flavor and the 
!!  segments of particle
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  flavor=hybridization function to take
!!  particles=segments to use
!!
!! OUTPUT
!!  BathOperator_getDetF=the det 
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION BathOperator_getDetF(op,flavor,particle)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_getDetF'
!End of the abilint section

  TYPE(BathOperator)       , INTENT(INOUT)      :: op
  INTEGER                  , INTENT(IN   )  :: flavor
  TYPE(ListCdagC), OPTIONAL, INTENT(IN   )  :: particle
!Local arguments-------------------------------
  INTEGER :: iCdag
  INTEGER :: iC
  INTEGER :: tail
  DOUBLE PRECISION :: time
  DOUBLE PRECISION :: tC
  DOUBLE PRECISION :: tCdag
  DOUBLE PRECISION :: beta
  DOUBLE PRECISION :: mbeta_two
  DOUBLE PRECISION :: signe
  DOUBLE PRECISION :: inv_dt
#include "BathOperator_hybrid.h"

  BathOperator_getDetF = 1.d0 ! pour eviter des divisions par 0
  IF ( PRESENT( particle ) ) THEN
    tail = particle%tail
    activeF = flavor
    beta = op%beta
    mbeta_two = -beta*0.5d0
    inv_dt =  op%inv_dt
    CALL MatrixHyb_setSize(op%M_update(flavor),tail)
    DO iCdag = 1, tail
      tCdag  = particle%list(iCdag,Cdag_)
      DO iC  = 1, tail
        !tC   = particle%list(C_,iC).MOD.beta
        MODCYCLE(particle%list(iC,C_),beta,tC)
        time = tC - tCdag
#include "BathOperator_hybrid"
        op%M_update(flavor)%mat(iC,iCdag) = hybrid 
      END DO
    END DO
    ! mat_tau needs to be transpose of ordered time mat (way of measuring
    ! G(tau))
    DO iC  = 1, tail
      tC   = particle%list(iC,C_)
      DO iCdag = 1, tail
        tCdag  = particle%list(iCdag,Cdag_)
        time = tC - tCdag
        signe = SIGN(1.d0,time)
        time = time + (signe-1.d0)*mbeta_two
        op%M_update(flavor)%mat_tau(iCdag,iC) = INT( ( time * inv_dt ) + 1.5d0 )
      END DO
    END DO
    CALL MatrixHyb_inverse(op%M_update(flavor),BathOperator_getDetF) ! calcul le det de la matrice et l'inverse
  ELSE
    CALL MatrixHyb_getDet(op%M(flavor),BathOperator_getDetF) ! det M = 1/detF !
    BathOperator_getDetF = 1.d0 / BathOperator_getDetF
  ENDIF
END FUNCTION BathOperator_getDetF
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_setMAdd
!! NAME
!!  BathOperator_setMAdd
!!
!! FUNCTION
!!  Update de M matrix inserting a row and a column
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  particle=segments of active flavor
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_setMAdd(op,particle) 

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_setMAdd'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
  TYPE(ListCdagC)   , INTENT(IN   ) :: particle
!Local variables ------------------------------
  INTEGER                           :: tail
  INTEGER                           :: new_tail
  INTEGER                           :: col
  INTEGER                           :: col_move
  INTEGER                           :: row_move
  INTEGER                           :: row
  INTEGER                           :: positionRow
  INTEGER                           :: positionCol
  INTEGER                           :: aF
  DOUBLE PRECISION                  :: Stilde
  DOUBLE PRECISION                  :: time
  DOUBLE PRECISION                  :: mbeta_two
  DOUBLE PRECISION                  :: inv_dt
  TYPE(Vector) :: vec_tmp
  TYPE(VectorInt) :: vecI_tmp
  INTEGER :: m
  INTEGER :: count
  INTEGER :: i
  INTEGER :: j
  INTEGER :: p

  IF ( op%MAddFlag .EQV. .FALSE. ) &
    CALL ERROR("BathOperator_setMAdd : MAddFlag turn off           ")
  af = op%activeFlavor
  IF ( aF .LE. 0 ) &
    CALL ERROR("BathOperator_setMAdd : no active hybrid function   ")
  tail     =  op%M(aF)%tail
  new_tail =  tail + 1
!CALL matrix_print(M)

  positionRow =  op%updatePosRow
  positionCol =  op%updatePosCol
  Stilde      =  op%Stilde
!  write(6,*) "before", positionRow, positionCol
  !CALL MatrixHyb_print(op%M(aF),opt_print=1)
  CALL MatrixHyb_setSize(op%M(aF),new_tail)

  ! Compute Qtilde with Q
  !op%Q%vec(1:tail) = (-1.d0) * MATMUL(op%M(aF)%mat(1:tail,1:tail),op%Q%vec(1:tail)) * Stilde
  op%Q%vec(1:tail) = MATMUL(op%M(aF)%mat(1:tail,1:tail),op%Q%vec(1:tail))
  !op%Q%vec(PositionRow:new_tail) = EOSHIFT(op%Q%vec(PositionRow:new_tail), SHIFT=-1, BOUNDARY=-1.d0, DIM=1)
!  op%Qtau%vec(PositionCol:new_tail) = EOSHIFT(op%Qtau%vec(PositionCol:new_tail), SHIFT=-1, BOUNDARY=1.d0, DIM=1)
!  op%Qtau%vec(PositionCol) = op%Stau

  !Compute Rtilde with R and without multiplying by Stilde
  !op%R%vec(1:tail) = (-1.d0) * MATMUL(op%R%vec(1:tail),op%M(aF)%mat(1:tail,1:tail))
  op%R%vec(1:tail) = MATMUL(op%R%vec(1:tail),op%M(aF)%mat(1:tail,1:tail))
  !op%R%vec(PositionCol:new_tail) = EOSHIFT(op%R%vec(PositionCol:new_tail), SHIFT=-1, BOUNDARY=-1.d0, DIM=1)
!  op%Rtau%vec(PositionRow:new_tail) = EOSHIFT(op%Rtau%vec(PositionRow:new_tail), SHIFT=-1, BOUNDARY=1.d0, DIM=1)
!  op%Rtau%vec(PositionRow) = op%Stau

  !Compute the new M matrix
  !op%M(aF)%mat(PositionRow:new_tail,1:new_tail) = &
  !                   EOSHIFT(op%M(aF)%mat(PositionRow:new_tail,1:new_tail),SHIFT=-1, BOUNDARY=0.d0, DIM=1)
  !op%M(aF)%mat(1:new_tail,PositionCol:new_tail) = &
  !                   EOSHIFT(op%M(aF)%mat(1:new_tail,PositionCol:new_tail),SHIFT=-1, BOUNDARY=0.d0, DIM=2)
! ! op%M(aF)%mat(1:new_tail,1:new_tail) =  op%M(aF)%mat(1:new_tail,1:new_tail) + &
! ! Stilde * MATMUL(RESHAPE(op%Q%vec(1:new_tail),(/ new_tail,1 /)),RESHAPE(op%R%vec(1:new_tail),(/ 1,new_tail /)))

  !op%M(aF)%mat_tau(PositionRow:new_tail,1:new_tail) = &
  !                   EOSHIFT(op%M(aF)%mat_tau(PositionRow:new_tail,1:new_tail),SHIFT=-1, BOUNDARY=0, DIM=1)
  !op%M(aF)%mat_tau(1:new_tail,PositionCol:new_tail) = &
  !                   EOSHIFT(op%M(aF)%mat_tau(1:new_tail,PositionCol:new_tail),SHIFT=-1, BOUNDARY=0, DIM=2)

  mbeta_two = -op%beta*0.5d0
  inv_dt = op%inv_dt
  !Shift mat_tau
  !update old m
  DO col=tail,1,-1
    col_move = col +  ( 1+SIGN(1,col-PositionCol) )/2
    DO row=tail,1,-1
      row_move = row +  ( 1+SIGN(1,row-PositionRow) )/2
      op%M(aF)%mat_tau(row_move,col_move) = op%M(aF)%mat_tau(row,col)
      op%M(aF)%mat(row_move,col_move) = op%M(aF)%mat(row,col) + op%Q%vec(row)*op%R%vec(col) * Stilde
    END DO
  END DO
  ! Add new stuff for new row
  DO row = 1, tail
    row_move = row +  ( 1+SIGN(1,row-PositionRow) )/2
    op%M(aF)%mat(row_move,PositionCol) = -op%Q%vec(row)*Stilde
    time = op%Rtau%vec(row)
    time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
    op%M(aF)%mat_tau(row,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
  END DO
  ! Add last time missing in the loops
  time = op%Rtau%vec(new_tail)
  time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
  op%M(aF)%mat_tau(new_tail,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
  ! Add new stuff for new col
  DO col = 1, tail 
    col_move = col +  ( 1+SIGN(1,col-PositionCol) )/2
    op%M(aF)%mat(PositionRow,col_move) = -op%R%vec(col)*Stilde
    time = op%Qtau%vec(col)
    time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
    op%M(aF)%mat_tau(PositionRow,col) = INT ( (time*inv_dt) +1.5d0 )
  END DO
  ! Add last time missing in the loops
  time = op%Qtau%vec(new_tail)
  time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
  op%M(aF)%mat_tau(PositionRow,new_tail) = INT ( (time*inv_dt) +1.5d0 )

  op%M(aF)%mat(PositionRow,PositionCol) = Stilde

  !CALL MatrixHyb_print(op%M(aF),opt_print=1)

!  DO col = 1, new_tail
!    time = op%Rtau%vec(col)
!    time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
!    op%M(aF)%mat_tau(col,PositionCol) = INT ( (time*inv_dt) +1.5d0 )
!    time = op%Qtau%vec(col)
!    time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
!    op%M(aF)%mat_tau(PositionRow,Col) = INT ( (time*inv_dt) +1.5d0 )
!    time = op%R%vec(col)*Stilde
!    DO row = 1, new_tail
!      op%M(aF)%mat(row,col) = op%M(aF)%mat(row,col) + op%Q%vec(row)*time
!    END DO
!  END DO

  !col_move = new_tail
  !col      = tail
  !DO col_move = new_tail, 1, -1
  !  IF ( col_move .EQ. positionCol ) THEN
  !    ! on calcule rajoute Q tilde
  !    !row_move = new_tail
  !    row      = tail 
  !    DO row_move = new_tail, 1, -1
  !      ! calcul itau
  !      IF ( row_move .EQ. positionRow ) THEN
  !        op%M(aF)%mat(row_move,col_move) = Stilde
  !        !time = op%Stau
  !      ELSE
  !        op%M(aF)%mat(row_move,col_move) = -op%Q%vec(row)*Stilde
  !        !time = op%Rtau%vec(row_move)
  !        row      = row      - 1 
  !      END IF
  !      !time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
  !      !op%M(aF)%mat_tau(row_move,col_move) = INT ( (time*inv_dt) +1.5d0 )
  !    END DO
  !    ! realignement des indices
  !  ELSE
  !    ! on calcule Ptilde
  !    !row_move = new_tail
  !    row      = tail 
  !    DO row_move = new_tail, 1, -1
  !      IF ( row_move .EQ. positionRow ) THEN
  !        op%M(aF)%mat(row_move,col_move) = -op%R%vec(col) * Stilde
  !        ! calcul itau
  !        !time = op%Qtau%vec(col_move)
  !        !time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
  !        !op%M(aF)%mat_tau(row_move,col_move) = INT ( (time*inv_dt) +1.5d0 )
  !      ELSE
  !        op%M(aF)%mat(row_move,col_move) = op%M(aF)%mat(row,col) + op%Q%vec(row)*op%R%vec(col)*Stilde
  !        ! copy itau
  !        !op%M(aF)%mat_tau(row_move,col_move) = op%M(aF)%mat_tau(row,col)
  !        row      = row      - 1 
  !      END IF
  !    END DO
  !    col      = col      - 1
  !  END IF
  !END DO
!  write(6,*) "after"
!  CALL MatrixHyb_print(op%M(aF),opt_print=1)
!CALL matrix_inverse(M)
!CALL MatrixHyb_print(M)
!CALL matrix_inverse(M)

  IF ( op%antiShift .EQV. .TRUE. ) THEN ! antisegment
    CALL Vector_init(vec_tmp,new_tail)
    CALL VectorInt_init(vecI_tmp,new_tail)
  ! Shift if necessary according to op%antishift
  ! shift DIM=2 (col)
    p = new_tail - 1
    m = 1
    count = 0
    DO WHILE ( count .NE. new_tail )
      vec_tmp%vec(1:new_tail) = op%M(aF)%mat(1:new_tail,m)
      vecI_tmp%vec(1:new_tail) = op%M(aF)%mat_tau(1:new_tail,m)
      i = m
      !j = m+p
      MODCYCLE(m+p, new_tail, j)
      DO WHILE (j .NE. m)
        op%M(aF)%mat(1:new_tail,i) = op%M(aF)%mat(1:new_tail,j)
        op%M(aF)%mat_tau(1:new_tail,i) = op%M(aF)%mat_tau(1:new_tail,j)
        i = j
        MODCYCLE(j+p, new_tail, j)
        count = count+1
      END DO
      op%M(aF)%mat(1:new_tail,i) = vec_tmp%vec(1:new_tail)
      op%M(aF)%mat_tau(1:new_tail,i) = vecI_tmp%vec(1:new_tail)
      count = count+1
      m = m+1
    END DO
    ! shift DIM=1 (row)
    p = new_tail - 1
    m = 1
    count = 0
    DO WHILE ( count .NE. new_tail)
      vec_tmp%vec(1:new_tail) = op%M(aF)%mat(m,1:new_tail)
      vecI_tmp%vec(1:new_tail) = op%M(aF)%mat_tau(m,1:new_tail)
      i = m
      !j = m+p
      MODCYCLE(m+p, new_tail, j)
      DO WHILE ( j .NE. m )
        op%M(aF)%mat(i,1:new_tail) = op%M(aF)%mat(j,1:new_tail)
        op%M(aF)%mat_tau(i,1:new_tail) = op%M(aF)%mat_tau(j,1:new_tail)
        i = j
        MODCYCLE(j+p, new_tail, j)
        count = count+1
      END DO
      op%M(aF)%mat(i,1:new_tail) = vec_tmp%vec(1:new_tail)
      op%M(aF)%mat_tau(i,1:new_tail) = vecI_tmp%vec(1:new_tail)
      count = count+1
      m = m+1
    END DO
    CALL Vector_destroy(vec_tmp)
    CALL VectorInt_destroy(vecI_tmp)
    !op%M(aF)%mat(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=-1, DIM=1) ! Shift to the bottom
    !op%M(aF)%mat(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=-1, DIM=2) ! Shift to the right
    !op%M(aF)%mat_tau(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=-1, DIM=1) ! Shift to the bottom
    !op%M(aF)%mat_tau(1:new_tail,1:new_tail) = CSHIFT(op%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=-1, DIM=2) ! Shift to the right
!CALL matrix_print(M)
  END IF

  IF ( op%doCheck .EQV. .TRUE.) THEN
!#ifdef CTQMC_CHECK
  CALL BathOperator_checkM(op,particle)
!#endif
  END IF

  op%MAddFlag = .FALSE.

END SUBROUTINE BathOperator_setMAdd
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_setMRemove
!! NAME
!!  BathOperator_setMRemove
!!
!! FUNCTION
!!  delete one row and one column of the M matrix
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  particle=segments of the active flavor
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_setMRemove(op,particle) 

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_setMRemove'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT)  :: op
  TYPE(ListCdagC)   , INTENT(IN   )  :: particle
!Local variables ------------------------------
  INTEGER                            :: tail
  INTEGER                            :: new_tail
  INTEGER                            :: col
  INTEGER                            :: col_move
  INTEGER                            :: row_move
  INTEGER                            :: row
  INTEGER                            :: positionCol
  INTEGER                            :: positionRow
  INTEGER                            :: aF
  INTEGER                              :: m
  INTEGER                              :: count
  INTEGER                              :: i
  INTEGER                              :: j
  INTEGER                              :: p
  DOUBLE PRECISION                   :: invStilde
  DOUBLE PRECISION                   :: invStilde2
  TYPE(VectorInt) :: vecI_tmp
  TYPE(Vector)    :: vec_tmp

  IF ( op%MRemoveFlag .EQV. .FALSE. ) &
    CALL ERROR("BathOperator_setMRemove : MRemoveFlag turn off     ")
  af = op%activeFlavor
  IF ( aF .LE. 0 ) &
    CALL ERROR("BathOperator_setMRemove : no active hybrid func    ")
  tail        =  op%M(aF)%tail
  new_tail    =  tail - 1
  positionCol =  op%updatePosCol
  positionRow =  op%updatePosRow
  invStilde   = 1.d0 / op%Stilde

!  write(6,*) "before", positionRow, positionCol
!  CALL MatrixHyb_print(op%M(aF),opt_print=1)

!  IF ( new_tail .EQ. 0 ) THEN
!!    IF ( op%antiShift .EQV. .TRUE.  ) THEN
!!      op%M(aF)%mat(1,1) = 1.d0/BathOperator_Hybrid(op, op%beta)
!!      op%MRemoveFlag = .FALSE.
!!      RETURN
!!    END IF
!    CALL MatrixHyb_clear(op%M(aF))
!    op%MRemoveFlag = .FALSE.
!    RETURN
!  END IF

!  CALL Vector_setSize(op%Q,new_tail)
!  CALL Vector_setSize(op%R,new_tail)
  Vector_QuickResize(op%Q,new_tail)
  Vector_QuickResize(op%R,new_tail)

!  We use R and Q as op%R%vec and op%Q%vec
!  op%R%vec => op%R
!  op%Q%vec => op%Q

  row      = 1
  !row_move = 1
  col      = 1
  !col_move = 1
  DO row_move = 1, new_tail
    IF ( row .EQ. positionRow ) row = row + 1
    IF ( col .EQ. positionCol ) col = col + 1
    !col = row_move + (1+SIGN(1,row_move-positionCol))/2
    !row = row_move + (1+SIGN(1,row_move-positionRow))/2
    op%R%vec(row_move) = op%M(aF)%mat(positionRow,col)
    op%Q%vec(row_move) = op%M(aF)%mat(row,positionCol)
    row      = row + 1 
    col      = col + 1
  END DO
!!    op%R%vec(1:positionCol-1) = op%M(aF)%mat(positionRow,1:positionCol-1)
!!    op%R%vec(positionCol:new_tail) = op%M(aF)%mat(positionRow,positionCol+1:tail)
!!    op%Q%vec(1:positionRow-1) = op%M(aF)%mat(1:positionRow-1,positionCol)
!!    op%Q%vec(positionRow:new_tail) = op%M(aF)%mat(positionRow+1:tail,positionCol)
!write(*,*) positionRow, positionCol
!CALL MatrixHyb_print(M)
!CALL Vector_print(op%R)
!CALL Vector_print(op%Q)
!CALL ListCdagC_print(op%ListCdagC)

  col      = 1
  DO col_move = 1, new_tail 
    IF ( col_move .EQ. positionCol ) col = col + 1
    !col = col_move + (1+SIGN(1,col_move-positionCol))/2
    row      = 1
    invStilde2 = invStilde * op%R%vec(col_move)
    DO row_move = 1, new_tail
      IF ( row_move .EQ. positionRow ) row = row + 1
      !row = row_move + (1+SIGN(1,row_move-positionRow))/2
      op%M(aF)%mat(row_move,col_move) = op%M(aF)%mat(row,col) &
                                      - op%Q%vec(row_move)*invStilde2
      op%M(aF)%mat_tau(row_move,col_move) = op%M(aF)%mat_tau(row,col)
      row      = row      + 1
    END DO
    col      = col      + 1 
  END DO
  CALL MatrixHyb_setSize(op%M(aF),new_tail)

  IF ( op%antiShift .EQV. .TRUE. ) THEN ! antisegment
    ! Shift if necessary according to op%antishift
    ! shift DIM=2 (col)
    CALL Vector_init(vec_tmp,new_tail)
    CALL VectorInt_init(vecI_tmp,new_tail)
    p = 1
    m = 1
    count = 0
    DO WHILE ( count .NE. new_tail )
      vec_tmp%vec(1:new_tail) = op%M(aF)%mat(1:new_tail,m)
      vecI_tmp%vec(1:new_tail) = op%M(aF)%mat_tau(1:new_tail,m)
      i = m
      !j = m+p
      MODCYCLE(m+p, new_tail, j)
      DO WHILE (j .NE. m)
        op%M(aF)%mat(1:new_tail,i) = op%M(aF)%mat(1:new_tail,j)
        op%M(aF)%mat_tau(1:new_tail,i) = op%M(aF)%mat_tau(1:new_tail,j)
        i = j
        MODCYCLE(j+p, new_tail, j)
        count = count+1
      END DO
      op%M(aF)%mat(1:new_tail,i) = vec_tmp%vec(1:new_tail)
      op%M(aF)%mat_tau(1:new_tail,i) = vecI_tmp%vec(1:new_tail)
      count = count+1
      m = m+1
    END DO
    CALL Vector_destroy(vec_tmp)
    CALL VectorInt_destroy(vecI_tmp)
    !op%M(aF)%mat(1:new_tail,1:new_tail) = &
    !           CSHIFT(op%M(aF)%mat(1:new_tail,1:new_tail), SHIFT=1, DIM=2) ! Shift to the top
    !op%M(aF)%mat_tau(1:new_tail,1:new_tail) = &
    !           CSHIFT(op%M(aF)%mat_tau(1:new_tail,1:new_tail), SHIFT=1, DIM=2) ! Shift to the top
  END IF
!  write(6,*) "after "
!  CALL MatrixHyb_print(op%M(aF),opt_print=1)

  IF ( op%doCheck .EQV. .TRUE. ) THEN
!#ifdef CTQMC_CHECK
  CALL BathOperator_checkM(op,particle)
!#endif
  END IF

  op%MRemoveFlag = .FALSE.

END SUBROUTINE BathOperator_setMRemove
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_swap
!! NAME
!!  BathOperator_swap
!!
!! FUNCTION
!!  Recompute 2 M matrix swaping the segments
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  iflavor1=flavor to swap with the next one
!!  iflavor2=favor to swap with the previous one
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_swap(op, flavor1, flavor2)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_swap'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
  INTEGER           , INTENT(IN   ) :: flavor1
  INTEGER           , INTENT(IN   ) :: flavor2

  !CALL MatrixHyb_print(op%M(flavor1),234)
  op%M(flavor1) = op%M_update(flavor1)
  !CALL MatrixHyb_print(op%M(flavor1),234)
  !CALL MatrixHyb_print(op%M(flavor2),234)
  op%M(flavor2) = op%M_update(flavor2)
  !CALL MatrixHyb_print(op%M(flavor2),234)

END SUBROUTINE BathOperator_swap
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_initF
!! NAME
!!  BathOperator_initF
!!
!! FUNCTION
!!  Copy input hybridization functions from a file
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  ifstream=file stream to read F
!!
!! OUTPUT
!!  argout(sizeout)=description
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_initF(op,ifstream)

!Arguments ----------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_initF'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
  INTEGER           , INTENT(IN   ) :: ifstream
!Local variables ----------------
  INTEGER                           :: flavor
  INTEGER                           :: sample

  IF ( op%set .EQV. .FALSE. ) &
    CALL ERROR("BathOperator_initF : BathOperator not set         ")

  DO flavor=1,op%flavors
    DO sample = 1, op%sizeHybrid
      READ(ifstream,*) op%F(sample,flavor)
    END DO
  END DO
END SUBROUTINE BathOperator_initF
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_setF
!! NAME
!!  BathOperator_setF
!!
!! FUNCTION
!!  Copy F from input array
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  F=array of the hybridization function
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_setF(op,F)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_setF'
!End of the abilint section

  TYPE(BathOperator)               , INTENT(INOUT) :: op
  DOUBLE PRECISION, DIMENSION(:,:) , INTENT(IN   ) :: F
!Arguments ------------------------------------
  INTEGER                                          :: flavor
  INTEGER                                          :: sample
  INTEGER                                          :: length

  IF ( op%set .EQV. .FALSE. ) &
    CALL ERROR("BathOperator_setF : BathOperator not set          ")

 length  = SIZE(F)
  IF ( length .NE. (op%flavors * op%sizeHybrid) ) &
    CALL ERROR("BathOperator_setF : wrong input F                 ")

  DO flavor=1,op%flavors
    DO sample = 1, op%sizeHybrid
    op%F(sample,flavor) = F(sample,flavor)
    END DO
  END DO
END SUBROUTINE BathOperator_setF
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_printF
!! NAME
!!  BathOperator_printF
!!
!! FUNCTION
!!  print F function
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  ostream=file stream to write in
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_printF(op,ostream)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_printF'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
  INTEGER,OPTIONAL  , INTENT(IN   ) :: ostream
!Local variables ------------------------------
  CHARACTER(LEN=4)                  :: aflavor
  CHARACTER(LEN=50)                  :: string
  INTEGER                           :: flavor
  INTEGER                           :: sample
  INTEGER                           :: ostream_val

  IF ( PRESENT(ostream) ) THEN 
    ostream_val = ostream
  ELSE  
    ostream_val = 65
    OPEN(UNIT=ostream_val, FILE="F.dat")
  END IF

  WRITE(aflavor,'(I4)') op%flavors+1
  string = '(1x,'//TRIM(ADJUSTL(aflavor))//'E22.14)'
  DO sample = 1, op%sizeHybrid
    WRITE(ostream_val,string) (sample-1)*op%dt, (op%F(sample,flavor), flavor=1,op%flavors)
  END DO
  !CALL FLUSH(ostream_val)

  IF ( .NOT. PRESENT(ostream) ) &
    CLOSE(ostream_val)

END SUBROUTINE BathOperator_printF
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_printM
!! NAME
!!  BathOperator_printM
!!
!! FUNCTION
!!  print M =F^{-1} matrix
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  ostream=file stream to write in
!!
!! OUTPUT
!!  argout(sizeout)=description
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_printM(op,ostream)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_printM'
!End of the abilint section

  TYPE(BathOperator), INTENT(IN) :: op
  INTEGER, OPTIONAL , INTENT(IN) :: ostream
!Local variables ------------------------------
  INTEGER                        :: ostream_val

  IF ( op%activeFlavor .LE. 0 ) &
    CALL ERROR("BathOperator_printM : no active hybrid function    ")
  ostream_val = 6
  IF ( PRESENT(ostream) ) ostream_val = ostream
  CALL MatrixHyb_print(op%M(op%activeFlavor),ostream_val)
END SUBROUTINE BathOperator_printM
!!***

!!****f* ABINIT/m_BathOperator/ BathOperator_destroy
!! NAME
!!   BathOperator_destroy
!!
!! FUNCTION
!!  Deallocate and reset every thing
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE  BathOperator_destroy(op)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_destroy'
!End of the abilint section

  TYPE(BathOperator), INTENT(INOUT) :: op
  INTEGER  :: it

  DO it = 1, op%flavors
    CALL MatrixHyb_destroy(op%M(it))
    CALL MatrixHyb_destroy(op%M_update(it))
  END DO

  CALL Vector_destroy(op%R)
  CALL Vector_destroy(op%Q)
  CALL Vector_destroy(op%Rtau)
  CALL Vector_destroy(op%Qtau)
  FREEIF(op%F)
  DT_FREEIF(op%M)
  DT_FREEIF(op%M_update)

  op%MAddFlag     = .FALSE.
  op%MRemoveFlag  = .FALSE.
  op%flavors      = 0 
  op%beta         = 0.d0
  op%dt      = 0.d0
  op%inv_dt  = 0.d0
  op%samples      = 0
  op%sizeHybrid   = 0
  op%activeFlavor = 0 
  op%updatePosRow = 0
  op%updatePosCol = 0

END SUBROUTINE BathOperator_destroy
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_doCheck
!! NAME
!!  BathOperator_doCheck
!!
!! FUNCTION
!!  Just store if we perfom check for updates of M
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  opt_check=second bit should be one
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_doCheck(op,opt_check)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_doCheck'
!End of the abilint section

  TYPE(BathOperator) , INTENT(INOUT) :: op
  INTEGER            , INTENT(IN   ) :: opt_check
  
  IF ( opt_check .GE. 2 ) &
    op%doCheck = .TRUE.
END SUBROUTINE BathOperator_doCheck
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_checkM
!! NAME
!!  BathOperator_checkM
!!
!! FUNCTION
!!  compute from scratch the M matrix and compar it
!!  with the already computed M matrix
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!  particle=list of all segments of the active flavor
!!
!! OUTPUT
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

SUBROUTINE BathOperator_checkM(op,particle)

!Arguments ------------------------------------

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_checkM'
!End of the abilint section

  TYPE(BathOperator) , INTENT(INOUT) :: op
  TYPE(ListCdagC)    , INTENT(IN   ) :: particle
!Local variables ------------------------------
!  TYPE(MatrixHyb)                    :: checkMatrix
  LOGICAL :: checkTau
  INTEGER :: tail
  INTEGER :: iC
  INTEGER :: iCdag
  INTEGER :: aF
  CHARACTER(LEN=4) :: a
  DOUBLE PRECISION :: time
  DOUBLE PRECISION :: beta
  DOUBLE PRECISION :: mbeta_two
  DOUBLE PRECISION :: erreur
  DOUBLE PRECISION :: tc
  DOUBLE PRECISION :: tCdag
  DOUBLE PRECISION :: sumMmat
  DOUBLE PRECISION :: sumCheck
#include "BathOperator_hybrid.h"

  aF = op%activeFlavor
  !Construction de la matrix
  tail = particle%tail
!  CALL MatrixHyb_init(checkMatrix,op%iTech,size=tail,Wmax=op%samples)
!  CALL MatrixHyb_setSize(checkMatrix,tail)
  CALL MatrixHyb_setSize(op%M_update(aF),tail)
  beta   =  op%beta
  mbeta_two = -beta*0.5d0
  op%checkNumber = op%checkNumber + 1
  IF ( tail .NE. op%M(aF)%tail ) THEN
    CALL WARN("BathOperator_checkM : tails are different          ")
    RETURN
  END IF

!CALL ListCdagC_print(particle)
  DO iCdag = 1, tail
    tCdag  = particle%list(iCdag,Cdag_)
    DO iC  = 1, tail
      !tC   = particle%list(C_,iC).MOD.beta
      MODCYCLE(particle%list(iC,C_),beta,tC)
      time = tC - tCdag
#include "BathOperator_hybrid"
      op%M_update(aF)%mat(iC,iCdag) = hybrid

      time = time + ( SIGN(1.d0,time) - 1.d0 )*mbeta_two
      op%M_update(aF)%mat_tau(iCdag,iC) = INT ( (time*op%inv_dt) +1.5d0 ) 
    END DO
  END DO

!    CALL MatrixHyb_Print(checkMatrix)
  !Inversion de la matrix
  CALL MatrixHyb_inverse(op%M_update(aF))
!    CALL MatrixHyb_Print(checkMatrix)

  !Comparaison
  sumMmat =0.d0
  sumCheck=0.d0
  erreur = 0.d0
  checkTau = .FALSE.
  DO iCdag = 1, tail
    Do iC =1, tail
      op%M_update(aF)%mat(iC,iCdag) = ABS((op%M_update(aF)%mat(iC, iCdag) - op%M(aF)%mat(iC,iCdag))/op%M(aF)%mat(iC,iCdag))
      IF ( op%M_update(aF)%mat(iC,iCdag) .GT. erreur ) erreur = op%M_update(aF)%mat(ic,iCdag)
      IF ( op%M_update(aF)%mat_tau(iC,iCdag) .NE. op%M(aF)%mat_tau(iC,iCdag) ) checkTau = .TRUE.
    END DO
  END DO

  IF ( checkTau .EQV. .TRUE. ) THEN
    CALL WARN("BathOperator_checkM : mat_tau differs should be")
    CALL MatrixHyb_print(op%M_update(aF),opt_print=1)
    CALL WARN("BathOperator_checkM : whereas it is")
    CALL MatrixHyb_print(op%M(aF),opt_print=1)
  END IF
  op%meanError = op%meanError + erreur
  IF ( erreur .GT. 1.d0 ) THEN 
    WRITE(a,'(I4)') INT(erreur*100.d0)
!    CALL MatrixHyb_Print(op%M(aF)
    CALL WARN("BathOperator_checkM : "//a//"%                        ") 
  END IF
!  CALL MatrixHyb_destroy(checkMatrix)
END SUBROUTINE BathOperator_checkM
!!***

!!****f* ABINIT/m_BathOperator/BathOperator_getError
!! NAME
!!  BathOperator_getError
!!
!! FUNCTION
!!  compute a percentage error / checkM
!!
!! COPYRIGHT
!!  Copyright (C) 2013-2014 ABINIT group (J. Bieder)
!!  This file is distributed under the terms of the
!!  GNU General Public License, see ~abinit/COPYING
!!  or http://www.gnu.org/copyleft/gpl.txt .
!!
!! INPUTS
!!  op=bath operator
!!
!! OUTPUT
!!  BathOperator_getError=Error in percent
!!
!! SIDE EFFECTS
!!
!! NOTES
!!
!! PARENTS
!!  Will be filled automatically by the parent script
!!
!! CHILDREN
!!  Will be filled automatically by the parent script
!!
!! SOURCE

DOUBLE PRECISION FUNCTION BathOperator_getError(op)


!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'BathOperator_getError'
!End of the abilint section

  TYPE(BathOperator), INTENT(IN) :: op

  IF ( op%doCheck .EQV. .TRUE. ) THEN
    BathOperator_getError = op%meanError / DBLE(op%checkNumber)
  ELSE
    BathOperator_getError = 0.d0
  END IF
END FUNCTION BathOperator_getError
!!***
!#endif

END MODULE m_BathOperator
!!***
