Commit da34051e authored by Florian Goth's avatar Florian Goth
Browse files

introduce and use a function for the frobeniusnorm

parent 650d4217
......@@ -474,7 +474,7 @@ function determinediagtype(nodes, nrnodes, mys) result(diagtype)
iszero = .true.
do i = 1, nrnodes
localzero = 1E-15*abs(nodes(i)%axy)*sqrt(2*(mys(nodes(i)%x)/abs(nodes(i)%axy) )**2 + 1)
localzero = 1E-15*abs(nodes(i)%axy)*sqrt(2.D0*(mys(nodes(i)%x)/abs(nodes(i)%axy) )**2 + 1.D0)
if (abs(mys(nodes(i)%x)) > localzero) iszero = .false.
enddo
if (iszero) then
......@@ -487,7 +487,7 @@ function determinediagtype(nodes, nrnodes, mys) result(diagtype)
istraceless = .true.
do i = 1, nrnodes
localzero = 1E-15*sqrt(mys(nodes(i)%x)**2 + mys(nodes(i)%y)**2 + dble(nodes(i)%axy*conjg(nodes(i)%axy)))
localzero = 1E-15*frobnorm(mys(nodes(i)%x), mys(nodes(i)%y), nodes(i)%axy)
if (abs(mys(nodes(i)%x) + mys(nodes(i)%y)) > localzero) istraceless = .false.
enddo
if (istraceless) then
......
......@@ -214,7 +214,7 @@ subroutine GeneralSingleColExp_init(this, nodes, nredges, mys, weight)
integer, intent(in) :: nredges
real (kind=kind(0.d0)), intent(in) :: weight
integer :: i
real (kind=kind(0.d0)) :: nf, my1, my2, localzero, md, mav, dweight
real (kind=kind(0.d0)) :: my1, my2, localzero, md, mav, dweight
allocate(this%x(2*nredges), this%y(nredges), this%c(2*nredges), this%s(nredges), this%cinv(2*nredges), this%sinv(nredges) )
allocate(this%c2(2*nredges), this%c2inv(2*nredges), this%s2(nredges), this%s2inv(nredges))
this%nrofentries = nredges
......@@ -228,8 +228,8 @@ subroutine GeneralSingleColExp_init(this, nodes, nredges, mys, weight)
!calculate Frobenius norm
my1 = mys(nodes(i)%x)
my2 = mys(nodes(i)%y)
nf = sqrt(my1*my1+my2*my2 + 2*dble(nodes(i)%axy * conjg(nodes(i)%axy))) ! dependence on weight cancels in all comps
localzero = 1E-15*nf ! definition of my local scale that defines zero
! dependence on weight cancels in all comps
localzero = 1E-15*frobnorm(my1, my2, nodes(i)%axy) ! definition of my local scale that defines zero
md = 0.5*(my1 - my2)
mav = 0.5*(my1 + my2)
if (abs(md) < localzero) then
......
......@@ -183,7 +183,7 @@ subroutine HomogeneousSingleColExp_init(this, nodes, nredges, mys, weight)
integer, intent(in) :: nredges
real (kind=kind(0.d0)), intent(in) :: weight
integer :: i
real (kind=kind(0.d0)) :: nf, my1, my2, localzero
real (kind=kind(0.d0)) :: my1, my2, localzero
allocate(this%x(2*nredges), this%y(nredges), this%c(nredges), this%s(nredges))
allocate(this%c2(nredges), this%s2(nredges))
allocate(this%c2inv(nredges), this%s2inv(nredges), this%cinv(nredges), this%sinv(nredges))
......@@ -198,8 +198,8 @@ subroutine HomogeneousSingleColExp_init(this, nodes, nredges, mys, weight)
!calculate Frobenius norm
my1 = mys(nodes(i)%x)
my2 = mys(nodes(i)%y)
nf = sqrt(my1*my1 + my2*my2 + 2*dble(nodes(i)%axy * conjg(nodes(i)%axy)))! dependence on weight drops out in all comps
localzero = 1E-15*nf ! definition of my local scale that defines zero
! dependence on weight drops out in all comparisons
localzero = 1E-15*frobnorm(my1, my2, nodes(i)%axy) ! definition of my local scale that defines zero
if (abs(my1-my2) > localzero) then
write(*,*) "[HomogeneousSingleColExp_init]: Unequal diagonals found. This should not happen here."
error stop 1
......
......@@ -272,4 +272,22 @@ pure subroutine lmultbase(c, s, x, nrofentries, mat)
endif
deallocate(xyarray, csh, snh)
end subroutine
!--------------------------------------------------------------------
!> @author
!> Florian Goth
!
!> @brief
!> A function to calculate the Frobenius norm of hermitian 2x2 matrices.
!
!> @param[in] d1 first diagonal entry
!> @param[in] d2 second diagonal entry
!> @param[in] o off-diagonal entry
!> @return The value of the frobenius norm
!--------------------------------------------------------------------
function frobnorm(d1, d2, o) result(fn)
real (kind=kind(0.d0)) :: fn, d1, d2
complex(kind=kind(0.D0)), intent(in) :: o
fn = sqrt(d1*d1+d2*d2 + 2*dble(o * conjg(o)))
end function
end module SingleColExpBase_mod
......@@ -390,7 +390,7 @@ subroutine TraceLessSingleColExp_init(this, nodes, nredges, mys, weight)
integer, intent(in) :: nredges
real (kind=kind(0.d0)), intent(in) :: weight
integer :: i
real (kind=kind(0.d0)) :: nf, my1, my2, localzero, tmp
real (kind=kind(0.d0)) :: my1, my2, localzero, tmp
! We need twice the amount of storage for the diagonal for this traceless case.
allocate(this%x(2*nredges), this%y(nredges), this%s(nredges), this%c(2*nredges))
allocate(this%c2(2*nredges), this%s2(nredges))
......@@ -405,8 +405,8 @@ subroutine TraceLessSingleColExp_init(this, nodes, nredges, mys, weight)
!calculate Frobenius norm
my1 = mys(nodes(i)%x)
my2 = mys(nodes(i)%y)
nf = sqrt(my1*my1+my2*my2 + 2*dble(nodes(i)%axy * conjg(nodes(i)%axy)))! dependence on weight cancels in all comps.
localzero = 1E-15*nf ! definition of my local scale that defines zero
! dependence on weight cancels in all comparisons.
localzero = 1E-15*frobnorm(my1, my2, nodes(i)%axy) ! definition of my local scale that defines zero
if (abs(my1+my2) > localzero) then
write(*,*) "[TraceLessSingleColExp_init]: Matrix not traceless. This should not happen here."
error stop 1
......
......@@ -245,7 +245,7 @@ subroutine ZeroDiagSingleColExp_init(this, nodes, nredges, mys, weight)
integer, intent(in) :: nredges
real (kind=kind(0.d0)), intent(in) :: weight
integer :: i
real (kind=kind(0.d0)) :: nf, my1, my2, localzero
real (kind=kind(0.d0)) :: my1, my2, localzero
allocate(this%x(2*nredges), this%y(nredges), this%c(nredges), this%s(nredges))
allocate(this%c2(nredges), this%s2(nredges))
this%nrofentries = nredges
......@@ -259,8 +259,8 @@ subroutine ZeroDiagSingleColExp_init(this, nodes, nredges, mys, weight)
!calculate Frobenius norm
my1 = mys(nodes(i)%x)
my2 = mys(nodes(i)%y)
nf = sqrt(my1*my1+my2*my2 + 2*dble(nodes(i)%axy * conjg(nodes(i)%axy)))! dependence on weight cancels in all comparisons.
localzero = 1E-15*nf ! definition of my local scale that defines zero
! dependence on weight cancels in all comparisons.
localzero = 1E-15*frobnorm(my1, my2, nodes(i)%axy) ! definition of my local scale that defines zero
if ((abs(my1) > localzero) .or. (abs(my2) > localzero)) then
write(*,*) "[ZeroDiagSingleColExp_init]: Diagonal NOT zero. This should not happen here."
error stop 1
......
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