Commit 44fb3f0d authored by Florian Goth's avatar Florian Goth
Browse files

dutifully introduce and use a function for the reverse adjoint

parent 982e0935
Pipeline #13041 passed with stages
in 25 minutes and 19 seconds
......@@ -65,6 +65,8 @@ module Exponentials_mod
procedure :: lmult_T => EulerExp_lmult_T
procedure :: adjoint => EulerExp_adjoint
procedure :: adjoint_T => EulerExp_adjoint_T
procedure :: reverseadjoint => EulerExp_reverseadjoint
procedure :: reverseadjoint_T => EulerExp_reverseadjoint_T
procedure :: adjoint_over_two => EulerExp_adjoint_over_two
procedure :: adjoint_over_two_T => EulerExp_adjoint_over_two_T
procedure :: rmultinv_T => EulerExp_rmultinv_T
......@@ -97,6 +99,7 @@ module Exponentials_mod
procedure :: rmultinv => FullExp_rmultinv
procedure :: lmult_T => FullExp_lmult_T
procedure :: adjoint => FullExp_adjoint
procedure :: reverseadjoint => FullExp_reverseadjoint
procedure :: adjoint_over_two => FullExp_adjoint_over_two
end type FullExp
......@@ -244,6 +247,16 @@ subroutine FullExp_adjoint(this, mat)
enddo
end subroutine FullExp_adjoint
subroutine FullExp_reverseadjoint(this, mat)
class(FullExp) :: this
complex(kind=kind(0.D0)), intent(inout) :: mat(:,:)
integer :: i
do i = this%evals-1, 1, -2
call this%stages(i+1)%reverseadjoint_T(mat)
call this%stages(i)%reverseadjoint(mat)
enddo
end subroutine FullExp_reverseadjoint
subroutine FullExp_lmultinv(this, mat)
class(FullExp) :: this
complex(kind=kind(0.D0)), intent(inout) :: mat(:,:)
......@@ -370,6 +383,24 @@ subroutine EulerExp_adjoint_T(this, mat)
enddo
end subroutine EulerExp_adjoint_T
subroutine EulerExp_reverseadjoint(this, mat)
class(EulerExp) :: this
complex(kind=kind(0.D0)), dimension(:, :) :: mat
integer :: i
do i = this%nrofcols, 1, -1
call this%singleexps(i)%dat%reverseadjointaction(mat)
enddo
end subroutine EulerExp_reverseadjoint
subroutine EulerExp_reverseadjoint_T(this, mat)
class(EulerExp) :: this
complex(kind=kind(0.D0)), dimension(:, :) :: mat
integer :: i
do i = 1,this%nrofcols
call this%singleexps(i)%dat%reverseadjointaction(mat)
enddo
end subroutine EulerExp_reverseadjoint_T
subroutine EulerExp_adjoint_over_two(this, mat)
class(EulerExp) :: this
complex(kind=kind(0.D0)), dimension(:, :) :: mat
......
......@@ -46,6 +46,7 @@ module GeneralSingleColExp_mod
procedure :: rmultinv => GeneralSingleColExp_rmultinv
procedure :: adjoint_over_two => GeneralSingleColExp_adjoint_over_two
procedure :: adjointaction => GeneralSingleColExp_adjointaction
procedure :: reverseadjointaction => GeneralSingleColExp_reverseadjointaction
end type
contains
......@@ -92,6 +93,14 @@ subroutine GeneralSingleColExp_adjointaction(this, mat)
call this%rmultinv(mat)
end subroutine GeneralSingleColExp_adjointaction
subroutine GeneralSingleColExp_reverseadjointaction(this, mat)
class(GeneralSingleColExp), intent(in) :: this
complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
call this%lmultinv(mat)
call this%rmult(mat)
end subroutine GeneralSingleColExp_reverseadjointaction
!--------------------------------------------------------------------
!> @author
!> Florian Goth
......
......@@ -48,6 +48,7 @@ module HomogeneousSingleColExp_mod
procedure :: rmultinv => HomogeneousSingleColExp_rmultinv
procedure :: adjoint_over_two => HomogeneousSingleColExp_adjoint_over_two
procedure :: adjointaction => HomogeneousSingleColExp_adjointaction
procedure :: reverseadjointaction => HomogeneousSingleColExp_reverseadjointaction
end type
contains
......@@ -94,6 +95,14 @@ subroutine HomogeneousSingleColExp_adjointaction(this, mat)
call this%rmultinv(mat)
end subroutine HomogeneousSingleColExp_adjointaction
subroutine HomogeneousSingleColExp_reverseadjointaction(this, mat)
class(HomogeneousSingleColExp), intent(in) :: this
complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
call this%lmultinv(mat)
call this%rmult(mat)
end subroutine HomogeneousSingleColExp_reverseadjointaction
subroutine HomogeneousSingleColExp_adjoint_over_two(this, mat)
class(HomogeneousSingleColExp), intent(in) :: this
complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
......
......@@ -54,6 +54,7 @@ module SingleColExpBase_mod
procedure(rmultinvinterface), deferred :: rmultinv
procedure(lmultinvinterface), deferred :: lmultinv
procedure(adjointactioninterface), deferred :: adjointaction
procedure(reverseadjointactioninterface), deferred :: reverseadjointaction
procedure(adjointactionovertwointerface), deferred :: adjoint_over_two
procedure(initinterface), deferred :: init
procedure(deallocinterface), deferred :: dealloc
......@@ -155,6 +156,19 @@ module SingleColExpBase_mod
class(SingleColExpBase), intent(inout) :: this
end subroutine
!--------------------------------------------------------------------
!> @brief
!> Perform the similarity transform e^{T} arg e^{-T}
!
!> @param[in] this
!> @param[inout] mat the matrix that we intend to transform.
!--------------------------------------------------------------------
subroutine reverseadjointactioninterface(this, mat)
import SingleColExpBase
class(SingleColExpBase), intent(in) :: this
Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: mat
end subroutine
!--------------------------------------------------------------------
!> @brief
!> Perform the similarity transform e^{-T} arg e^{T}
......@@ -167,7 +181,7 @@ module SingleColExpBase_mod
class(SingleColExpBase), intent(in) :: this
Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: mat
end subroutine
!--------------------------------------------------------------------
!> @brief
!> Perform the similarity transform e^{-T/2} arg e^{T/2}
......
......@@ -46,6 +46,7 @@ module TraceLessSingleColExp_mod
procedure :: rmultinv => TraceLessSingleColExp_rmultinv
procedure :: adjoint_over_two => TraceLessSingleColExp_adjoint_over_two
procedure :: adjointaction => TraceLessSingleColExp_adjointaction
procedure :: reverseadjointaction => TraceLessSingleColExp_reverseadjointaction
end type
contains
......@@ -217,6 +218,30 @@ subroutine TraceLessSingleColExp_adjointaction(this, mat)
call this%rmultinv(mat)
end subroutine TraceLessSingleColExp_adjointaction
!--------------------------------------------------------------------
!> @author
!> Florian Goth
!
!> @brief
!> The routines for moving to an adjoint representation : out = this.mat.this^(-1)
!> If needed we could instead calculate an eigendecomposition and use that.
!> We could really invest in a diagonal calculation at every multiplication
!> The stability of this topic has been discussed in
!> Hargreaves, G. (2005). Topics in matrix computations: Stability and efficiency of algorithms (Doctoral dissertation, University of Manchester).
!> and "Unifying unitary and hyperbolic transformations Adam Bojanczyka, Sanzheng Qiaob;;1, Allan O. Steinhardt"
!> For the future we might want to look into fast hyperbolic rotations of Hargreaves, G. (2005).
!
!> @param[in] this The exponential that we consider
!> @param[inout] mat the matrix that we modify.
!--------------------------------------------------------------------
subroutine TraceLessSingleColExp_reverseadjointaction(this, mat)
class(TraceLessSingleColExp), intent(in) :: this
complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
call this%lmultinv(mat)
call this%rmult(mat)
end subroutine TraceLessSingleColExp_reverseadjointaction
subroutine TraceLessSingleColExp_adjoint_over_two(this, mat)
class(TraceLessSingleColExp), intent(in) :: this
complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
......
......@@ -46,6 +46,7 @@ module ZeroDiagSingleColExp_mod
procedure :: rmultinv => ZeroDiagSingleColExp_rmultinv
procedure :: adjoint_over_two => ZeroDiagSingleColExp_adjoint_over_two
procedure :: adjointaction => ZeroDiagSingleColExp_adjointaction
procedure :: reverseadjointaction => ZeroDiagSingleColExp_reverseadjointaction
end type
contains
......@@ -137,6 +138,30 @@ subroutine ZeroDiagSingleColExp_adjointaction(this, mat)
call this%rmultinv(mat)
end subroutine ZeroDiagSingleColExp_adjointaction
!--------------------------------------------------------------------
!> @author
!> Florian Goth
!
!> @brief
!> The routines for moving to an adjoint representation : out = this.mat.this^(-1)
!> If needed we could instead calculate an eigendecomposition and use that.
!> We could really invest in a diagonal calculation at every multiplication
!> The stability of this topic has been discussed in
!> Hargreaves, G. (2005). Topics in matrix computations: Stability and efficiency of algorithms (Doctoral dissertation, University of Manchester).
!> and "Unifying unitary and hyperbolic transformations Adam Bojanczyka, Sanzheng Qiaob;;1, Allan O. Steinhardt"
!> For the future we might want to look into fast hyperbolic rotations of Hargreaves, G. (2005).
!
!> @param[in] this The exponential that we consider.
!> @param[inout] mat the matrix that we modify.
!--------------------------------------------------------------------
subroutine ZeroDiagSingleColExp_reverseadjointaction(this, mat)
class(ZeroDiagSingleColExp), intent(in) :: this
complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
call this%lmultinv(mat)
call this%rmult(mat)
end subroutine ZeroDiagSingleColExp_reverseadjointaction
subroutine ZeroDiagSingleColExp_adjoint_over_two(this, mat)
class(ZeroDiagSingleColExp), intent(in) :: this
complex(kind=kind(0.D0)), dimension(:, :), intent(inout) :: mat
......
......@@ -45,6 +45,7 @@ module ContainerElementBase_mod
procedure(rmultinvinterface), deferred :: rmultinv
procedure(lmultinvinterface), deferred :: lmultinv
procedure(adjointactioninterface), deferred :: adjoint
procedure(reverseadjointinterface), deferred :: reverseadjoint
procedure(adjoint_over_twointerface), deferred :: adjoint_over_two
procedure(dump), deferred :: dump
procedure(dealloc), deferred :: dealloc
......@@ -135,6 +136,19 @@ module ContainerElementBase_mod
class(ContainerElementBase), intent(in) :: this
Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: arg
end subroutine
!--------------------------------------------------------------------
!> @brief
!> Perform the similarity transform e^{T} arg e^{-T}
!
!> @param[in] this
!> @param[inout] the matrix that we intend to transform.
!--------------------------------------------------------------------
subroutine reverseadjointinterface(this, arg)
import ContainerElementBase
class(ContainerElementBase), intent(in) :: this
Complex(kind=kind(0.d0)), intent(inout), dimension(:,:) :: arg
end subroutine
!--------------------------------------------------------------------
!> @brief
......
......@@ -249,6 +249,46 @@
end Subroutine Hop_mod_mmthl
!--------------------------------------------------------------------
Subroutine Hop_mod_adjoint (In,nf)
!e^U M e^{-U}
! InOut: In = IN * e^{ -dtau T }
Implicit none
Complex (Kind=Kind(0.d0)), intent(INOUT) :: IN(:,:)
Integer :: nf
!Local
Integer :: nc
class(ContainerElementBase), pointer :: dummy
do nc = Ncheck, 1, -1
dummy => ExpOpT_vec(nf)%at(nc)
call dummy%adjoint(In)
Enddo
end Subroutine Hop_mod_adjoint
!----------------------------------------------------------------------
Subroutine Hop_mod_reverseadjoint (In, nf)
!e^{-U} M e^U
Implicit none
Complex (Kind=Kind(0.d0)), intent(INOUT) :: IN(:,:)
Integer :: nf
!Local
Integer :: nc
class(ContainerElementBase), pointer :: dummy
do nc = Ncheck, 1, -1
dummy => ExpOpT_vec(nf)%at(nc)
call dummy%reverseadjoint(In)
Enddo
end Subroutine Hop_mod_reverseadjoint
!----------------------------------------------------------------------
Subroutine Hop_mod_mmthlc (In,nf)
......
......@@ -59,6 +59,7 @@ module OpTTypes_mod
procedure :: rmultinv => RealExpOpT_rmultinv ! right multiplication with Op_T inverse
procedure :: lmultinv => RealExpOpT_lmultinv
procedure :: adjoint => RealExpOpT_adjoint
procedure :: reverseadjoint => RealExpOpT_reverseadjoint
procedure :: adjoint_over_two => RealExpOpT_adjoint_over_two
procedure :: dump => RealExpOpT_dump ! dump matrices for debugging to screen
end type RealExpOpT
......@@ -84,6 +85,7 @@ module OpTTypes_mod
procedure :: rmultinv => CmplxExpOpT_rmultinv ! right multiplication with Op_T inverse
procedure :: lmultinv => CmplxExpOpT_lmultinv
procedure :: adjoint => CmplxExpOpT_adjoint
procedure :: reverseadjoint => CmplxExpOpT_reverseadjoint
procedure :: adjoint_over_two => CmplxExpOpT_adjoint_over_two
procedure :: dump => CmplxExpOpT_dump ! dump matrices for debugging to screen
end type CmplxExpOpT
......@@ -158,6 +160,19 @@ contains
Endif
end subroutine
subroutine RealExpOpT_reverseadjoint(this, arg)
class(RealExpOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
Integer :: n1, n2
n1 = size(arg,1)
n2 = size(arg,2)
If ( this%g*this%g > this%Zero ) then
call ZDSLSYMM('L', 'U', this%Ndim_hop, n1, n2, this%invmat, this%P, arg)
call ZDSLSYMM('R', 'U', this%Ndim_hop, n1, n2, this%mat, this%P, arg)
Endif
end subroutine
subroutine RealExpOpT_rmult(this, arg)
class(RealExpOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
......@@ -265,7 +280,20 @@ contains
call ZSLHEMM('R', 'U', this%Ndim_hop, n1, n2, this%invmat, this%P, arg)
Endif
end subroutine
subroutine CmplxExpOpT_reverseadjoint(this, arg)
class(CmplxExpOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
Integer :: n1, n2
n1 = size(arg,1)
n2 = size(arg,2)
If ( dble(this%g*conjg(this%g)) > this%Zero ) then
call ZSLHEMM('L', 'U', this%Ndim_hop, n1, n2, this%invmat, this%P, arg)
call ZSLHEMM('R', 'U', this%Ndim_hop, n1, n2, this%mat, this%P, arg)
Endif
end subroutine
subroutine CmplxExpOpT_rmult(this, arg)
class(CmplxExpOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
......
......@@ -128,8 +128,9 @@ Contains
! Wrap up, upgrade ntau1. with B^{1}(tau1)
NTAU1 = NTAU + 1
Do nf = 1,N_FL
CALL HOP_MOD_mmthr (GR(:,:,nf), nf )
CALL HOP_MOD_mmthl_m1(GR(:,:,nf), nf )
call HOP_MOD_adjoint(GR(:,:,nf), nf )
! CALL HOP_MOD_mmthr (GR(:,:,nf), nf )
! CALL HOP_MOD_mmthl_m1(GR(:,:,nf), nf )
Enddo
Do n = Nt_sequential_start,Nt_sequential_end
Do nf = 1, N_FL
......@@ -273,10 +274,11 @@ Contains
enddo
enddo
DO nf = 1,N_FL
Call Hop_mod_mmthl (GR(:,:,nf), nf)
Call Hop_mod_mmthr_m1(GR(:,:,nf), nf)
call Hop_mod_reverseadjoint(GR(:,:,nf), nf)
! Call Hop_mod_mmthl (GR(:,:,nf), nf)
! Call Hop_mod_mmthr_m1(GR(:,:,nf), nf)
enddo
end SUBROUTINE WRAPGRDO
......
......@@ -58,8 +58,8 @@ module mscbOpT_mod
procedure :: lmult => CmplxmscbOpT_lmult ! left multiplication with Op_T
procedure :: rmultinv => CmplxmscbOpT_rmultinv ! right multiplication with Op_T inverse
procedure :: lmultinv => CmplxmscbOpT_lmultinv ! left multiplication with Op_T inverse
procedure :: adjointaction => CmplxmscbOpT_adjointaction
procedure :: adjoint => CmplxmscbOpT_adjoint
procedure :: reverseadjoint => CmplxmscbOpT_reverseadjoint
procedure :: adjoint_over_two => CmplxmscbOpT_adjoint_over_two
procedure :: dump => CmplxmscbOpT_dump ! dump matrices for debugging to screen
end type CmplxmscbOpT
......@@ -84,9 +84,9 @@ module mscbOpT_mod
procedure :: lmult => CmplxEulermscbOpT_lmult
procedure :: rmultinv => CmplxEulermscbOpT_rmultinv ! right multiplication with Op_T inverse
procedure :: lmultinv => CmplxEulermscbOpT_lmultinv
procedure :: adjointaction => CmplxEulermscbOpT_adjointaction
procedure :: adjoint => CmplxEulermscbOpT_adjointaction
procedure :: adjoint_over_two => CmplxEulermscbOpT_adjointaction
procedure :: adjoint => CmplxEulermscbOpT_adjoint
procedure :: reverseadjoint => CmplxEulermscbOpT_reverseadjoint
procedure :: adjoint_over_two => CmplxEulermscbOpT_adjoint_over_two
procedure :: dump => CmplxEulermscbOpT_dump ! dump matrices for debugging to screen
end type CmplxEulermscbOpT
......@@ -168,16 +168,16 @@ contains
deallocate(diags)
end subroutine
subroutine CmplxmscbOpT_adjointaction(this, arg)
subroutine CmplxmscbOpT_reverseadjoint(this, arg)
class(CmplxmscbOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
!FIXME: P
If ( dble(this%g*conjg(this%g)) > this%Zero ) then
call this%fe%adjoint_over_two(arg)
call this%fe%reverseadjoint(arg)
Endif
end subroutine
subroutine CmplxmscbOpT_adjoint(this, arg)
class(CmplxmscbOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
......@@ -292,8 +292,30 @@ contains
call dealloc_graphdata(gd)
deallocate(diags)
end subroutine
subroutine CmplxEulermscbOpT_adjointaction(this, arg)
subroutine CmplxEulermscbOpT_reverseadjoint(this, arg)
class(CmplxEulermscbOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
!FIXME: P
If ( dble(this%g*conjg(this%g)) > this%Zero ) then
call this%ee%reverseadjoint(arg)
Endif
end subroutine
subroutine CmplxEulermscbOpT_adjoint(this, arg)
class(CmplxEulermscbOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
!FIXME: P
If ( dble(this%g*conjg(this%g)) > this%Zero ) then
call this%ee%adjoint(arg)
Endif
end subroutine
subroutine CmplxEulermscbOpT_adjoint_over_two(this, arg)
class(CmplxEulermscbOpT), intent(in) :: this
Complex(kind=kind(0.D0)), intent(inout), dimension(:,:) :: arg
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment