!===============================================================================
!
! BQCD -- Berlin Quantum ChromoDynamics programme
!
! Author: Hinnerk Stueben <stueben@zib.de>
!
! Copyright (C) 1998-2006, Hinnerk Stueben, Zuse-Institut Berlin
!
!-------------------------------------------------------------------------------
!
! dsf1.F90  --  p(j,x,mu) := p(j,x,mu) - step * D_{x,mu,j} S_{f1}
!
!-------------------------------------------------------------------------------
# include "defs.h"

!-------------------------------------------------------------------------------
subroutine dsf1(p, conf, step, calc_sf, sf, para)

  use typedef_hmc
  use module_hmc_forces
  use module_mre
  use module_p_interface
  use module_switches
  use module_vol
  implicit none

  GENERATOR_FIELD, intent(inout) :: p
  type(hmc_conf),  intent(in)    :: conf
  type(hmc_para),  intent(in)    :: para
  integer,         intent(in)    :: calc_sf
  REAL,            intent(in)    :: step
  REAL,            intent(out)   :: sf

  type(type_mre),  save          :: solutions
  P_SPINCOL_FIELD, save          :: a, b
  REAL, external                 :: dotprod
  integer                        :: iterations
  external                       :: w_mult
  external                       :: w_dagger_w

  sf = ZERO

  if (switches%quenched) return

  ALLOCATE_SC_FIELD(a)
  ALLOCATE_SC_FIELD(b)

  call flip_bc(conf%u)

  call mre_get(solutions, w_mult, a, conf%phi, para, conf)
  call cg(w_dagger_w, a, conf%phi, para, conf, iterations) ! A = inv(W+ W~) Phi
  call mre_put(solutions, a, calc_sf)                      ! calc_sf <=> reset
  call w_mult(b, a, para, conf)                            ! B = W~ A

  if (calc_sf /= 0) sf = dotprod(b, b, SIZE_SC_FIELD)

  call hmc_forces_old(p)
  call dsf(p, conf, step, para, a, b)
  call hmc_forces_new(p, step, i_sf1)  

  !! call flip_bc(conf%u) <- done in dsf()

  call iteration_count_f1(iterations)
end

!===============================================================================
