psym2.f90

Go to the documentation of this file.

 
 
module psym2
use pkind, only: spi,dp
use pietc, only: u0,u1,o2
implicit none
private
public:: eigensym2,invsym2,sqrtsym2,expsym2,logsym2,id2222,chol2
 
real(dp),dimension(2,2,2,2):: id
data id/u1,u0,u0,u0, u0,o2,o2,u0,  u0,o2,o2,u0, u0,u0,u0,u1/! Effective identity
 
interface eigensym2;   module procedure eigensym2,eigensym2d; end interface
interface invsym2;     module procedure invsym2,invsym2d;     end interface
interface sqrtsym2;    module procedure sqrtsym2,sqrtsym2d;   end interface
interface sqrtsym2d_e; module procedure sqrtsym2d_e;          end interface
interface sqrtsym2d_t; module procedure sqrtsym2d_t;          end interface
interface expsym2;     module procedure expsym2,expsym2d;     end interface
interface expsym2d_e;  module procedure expsym2d_e;           end interface
interface expsym2d_t;  module procedure expsym2d_t;           end interface
interface logsym2;     module procedure logsym2,logsym2d;     end interface
interface id2222;      module procedure id2222;               end interface
interface chol2;       module procedure chol2;                end interface
 
contains
 
subroutine eigensym2(em,vv,oo)!                                    [eigensym2]
implicit none
real(dp),dimension(2,2),intent(in ):: em
real(dp),dimension(2,2),intent(out):: vv,oo
real(dp):: a,b,c,d,e,f,g,h
a=em(1,1); b=em(1,2); c=em(2,2)
d=a*c-b*b! <- det(em)
e=(a+c)*o2; f=(a-c)*o2
h=sqrt(f**2+b**2)
g=sqrt(b**2+(h+abs(f))**2)
if    (g==u0)then; vv(:,1)=(/u1,u0/)
elseif(f> u0)then; vv(:,1)=(/h+f,b/)/g
else             ; vv(:,1)=(/b,h-f/)/g
endif      
vv(:,2)=(/-vv(2,1),vv(1,1)/)
oo=matmul(transpose(vv),matmul(em,vv))
oo(1,2)=u0; oo(2,1)=u0
end subroutine eigensym2
 
subroutine eigensym2d(em,vv,oo,vvd,ood,ff)!                        [eigensym2]
implicit none
real(dp),dimension(2,2),    intent(in ):: em
real(dp),dimension(2,2),    intent(out):: vv,oo
real(dp),dimension(2,2,2,2),intent(out):: vvd,ood
logical,                    intent(out):: ff
real(dp),dimension(2,2):: emd,tt,vvr
real(dp)               :: oodif,dtheta
integer(spi)           :: i,j
call eigensym2(em,vv,oo); vvr(1,:)=-vv(2,:); vvr(2,:)=vv(1,:)
oodif=oo(1,1)-oo(2,2); ff=oodif==u0; if(ff)return
ood=0
vvd=0
do j=1,2
   do i=1,2
      emd=0
      if(i==j)then
         emd(i,j)=u1
      else
         emd(i,j)=o2; emd(j,i)=o2
      endif
      tt=matmul(transpose(vv),matmul(emd,vv))
      ood(1,1,i,j)=tt(1,1)
      ood(2,2,i,j)=tt(2,2)
      dtheta=tt(1,2)/oodif
      vvd(:,:,i,j)=vvr*dtheta
   enddo
enddo
end subroutine eigensym2d
 
subroutine invsym2(em,z)!                                            [invsym2]
implicit none
real(dp),dimension(2,2),intent(in ):: em
real(dp),dimension(2,2),intent(out):: z
real(dp):: detem
z(1,1)=em(2,2); z(2,1)=-em(2,1); z(1,2)=-em(1,2); z(2,2)=em(1,1)
detem=em(1,1)*em(2,2)-em(2,1)*em(1,2)
z=z/detem
end subroutine invsym2
 
subroutine invsym2d(em,z,zd)!                                        [invsym2]
implicit none
real(dp),dimension(2,2)    ,intent(in ):: em
real(dp),dimension(2,2)    ,intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
integer(spi):: k,l
call invsym2(em,z)
call id2222(zd)
do l=1,2; do k=1,2
   zd(:,:,k,l)=-matmul(matmul(z,zd(:,:,k,l)),z)
enddo;    enddo
end subroutine invsym2d
 
subroutine sqrtsym2(em,z)!                                          [sqrtsym2]
implicit none
real(dp),dimension(2,2),intent(in ):: em
real(dp),dimension(2,2),intent(out):: z
real(dp),dimension(2,2):: vv,oo
integer(spi)           :: i
call eigensym2(em,vv,oo)
do i=1,2
if(oo(i,i)<0)stop 'In sqrtsym2; matrix em is not non-negative'
oo(i,i)=sqrt(oo(i,i)); enddo
z=matmul(vv,matmul(oo,transpose(vv)))
end subroutine sqrtsym2
 
subroutine sqrtsym2d(x,z,zd)!                                        [sqrtsym2]
implicit none
real(dp),dimension(2,2),    intent(in ):: x
real(dp),dimension(2,2),    intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
real(dp),dimension(2,2):: px
real(dp)               :: rdetx,lrdetx,htrpxs,q
rdetx=sqrt(x(1,1)*x(2,2)-x(1,2)*x(2,1)) ! <- sqrt(determinant of x)
lrdetx=sqrt(rdetx)
px=x/rdetx                  ! - preconditioned x (has unit determinant)
htrpxs= ((px(1,1)+px(2,2))/2)**2 ! - {half-trace-px}-squared
q=htrpxs-u1
if(q<.05_dp)then               ! - Taylor expansion method
   call sqrtsym2d_t(px,z,zd)
   z=z*lrdetx; zd=zd/lrdetx
else
   call sqrtsym2d_e(x,z,zd) ! - Eigen-method
endif
end subroutine sqrtsym2d
 
subroutine sqrtsym2d_e(x,z,zd)!                                  [sqrtsym2d_e]
implicit none
real(dp),dimension(2,2),    intent(in ):: x
real(dp),dimension(2,2),    intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
real(dp),dimension(2,2,2,2):: vvd,ood
real(dp),dimension(2,2)    :: vv,oo,oori,tt
integer(spi)               :: i,j
logical                    :: ff
call eigensym2(x,vv,oo,vvd,ood,ff)
z=u0; z(1,1)=sqrt(oo(1,1)); z(2,2)=sqrt(oo(2,2))
z=matmul(matmul(vv,z),transpose(vv))
oori=u0; oori(1,1)=u1/sqrt(oo(1,1)); oori(2,2)=u1/sqrt(oo(2,2))
do j=1,2
do i=1,2
   tt=matmul(vvd(:,:,i,j),transpose(vv))
   zd(:,:,i,j)=o2*matmul(matmul(matmul(vv,ood(:,:,i,j)),oori),transpose(vv))&
        +matmul(tt,z)-matmul(z,tt)
enddo
enddo
end subroutine sqrtsym2d_e
 
subroutine sqrtsym2d_t(x,z,zd)!                                  [sqrtsym2d_t]
implicit none
real(dp),dimension(2,2),    intent(in ):: x
real(dp),dimension(2,2),    intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
integer(spi),parameter     :: nit=300 ! number of iterative increments allowed
real(dp),parameter         :: crit=1.e-17
real(dp),dimension(2,2)    :: r,rp,rd,rpd
real(dp)                   :: c
integer(spi)               :: i,j,n
r=x; r(1,1)=x(1,1)-1; r(2,2)=x(2,2)-1
z=u0; z(1,1)=u1; z(2,2)=u1
rp=r
c=o2
do n=0,nit
   z=z+c*rp
   rp=matmul(rp,r)
   if(sum(abs(rp))<crit)exit
   c=-c*(n*2+1)/(2*(n+2))
enddo
do j=1,2; do i=1,2
   rd=id(:,:,i,j)
   rpd=rd
   zd(:,:,i,j)=0
   rp=r
   c=o2
   do n=0,nit
      zd(:,:,i,j)=zd(:,:,i,j)+c*rpd
      rpd=matmul(rd,rp)+matmul(r,rpd)
      rp=matmul(r,rp)
      if(sum(abs(rp))<crit)exit
      c=-c*(n*2+1)/(2*(n+2))
   enddo
enddo; enddo
end subroutine sqrtsym2d_t
 
subroutine expsym2(em,expem)!                                        [expsym2]
implicit none
real(dp),dimension(2,2),intent(in ):: em
real(dp),dimension(2,2),intent(out):: expem
real(dp),dimension(2,2):: vv,oo
integer(spi)           :: i
call eigensym2(em,vv,oo)
do i=1,2; oo(i,i)=exp(oo(i,i)); enddo
expem=matmul(vv,matmul(oo,transpose(vv)))
end subroutine expsym2
 
subroutine expsym2d(x,z,zd)!                                         [expsym2]
implicit none
real(dp),dimension(2,2),    intent(in ):: x
real(dp),dimension(2,2),    intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
real(dp),dimension(2,2):: px
real(dp)               :: trxh,detpx
trxh=(x(1,1)+x(2,2))*o2
px=x;px(1,1)=x(1,1)-trxh;px(2,2)=x(2,2)-trxh
detpx=abs(px(1,1)*px(2,2)-px(1,2)*px(2,1))
if(detpx<.1_dp)then; call expsym2d_t(px,z,zd)
else               ; call expsym2d_e(px,z,zd)
endif
z=z*exp(trxh)
end subroutine expsym2d
 
subroutine expsym2d_e(x,z,zd)!                                    [expsym2d_e]
implicit none
real(dp),dimension(2,2),    intent(in ):: x
real(dp),dimension(2,2),    intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
real(dp),dimension(2,2,2,2):: vvd,ood
real(dp),dimension(2,2)    :: vv,oo,ooe,tt
integer(spi)               :: i,j
logical                    :: ff
call eigensym2(x,vv,oo,vvd,ood,ff)
z=u0; z(1,1)=exp(oo(1,1)); z(2,2)=exp(oo(2,2))
z=matmul(matmul(vv,z),transpose(vv))
ooe=u0; ooe(1,1)=exp(oo(1,1)); ooe(2,2)=exp(oo(2,2))
do j=1,2
do i=1,2
   tt=matmul(vvd(:,:,i,j),transpose(vv))
   zd(:,:,i,j)=matmul(matmul(matmul(vv,ood(:,:,i,j)),ooe),transpose(vv))&
        +matmul(tt,z)-matmul(z,tt)
enddo
enddo
end subroutine expsym2d_e
 
subroutine expsym2d_t(x,z,zd)!                                    [expsym2d_t]
implicit none
real(dp),dimension(2,2),    intent(in ):: x
real(dp),dimension(2,2),    intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
integer(spi),parameter     :: nit=100 ! number of iterative increments allowed
real(dp),parameter         :: crit=1.e-17_dp
real(dp),dimension(2,2)    :: xp,xd,xpd
real(dp)                   :: c
integer(spi)               :: i,j,n
z=0; z(1,1)=u1; z(2,2)=u1
xp=x
c=u1
do n=1,nit
   z=z+c*xp
   xp=matmul(xp,x)
   if(sum(abs(xp))<crit)exit
   c=c/(n+1)
enddo
do j=1,2; do i=1,2
   xd=id(:,:,i,j)
   xpd=xd
   zd(:,:,i,j)=0
   xp=x
   c=u1
   do n=1,nit
      zd(:,:,i,j)=zd(:,:,i,j)+c*xpd
      xpd=matmul(xd,xp)+matmul(x,xpd)
      xp=matmul(x,xp)
      if(sum(abs(xpd))<crit)exit
      c=c/(n+1)
   enddo
enddo; enddo
end subroutine expsym2d_t
 
subroutine logsym2(em,logem)!                                        [logsym2]
implicit none
real(dp),dimension(2,2),intent(in ):: em
real(dp),dimension(2,2),intent(out):: logem
real(dp),dimension(2,2):: vv,oo
integer(spi)           :: i
call eigensym2(em,vv,oo)
do i=1,2
   if(oo(i,i)<=u0)stop 'In logsym2; matrix em is not positive definite'
   oo(i,i)=log(oo(i,i))
enddo
logem=matmul(vv,matmul(oo,transpose(vv)))
end subroutine logsym2
 
subroutine logsym2d(x,z,zd)!                                         [logsym2]
use pfun, only: sinhox
implicit none
real(dp),dimension(2,2),    intent(in ):: x
real(dp),dimension(2,2),    intent(out):: z
real(dp),dimension(2,2,2,2),intent(out):: zd
real(dp),dimension(2,2):: vv,oo,d11,d12,d22,pqr,d11pqr,d12pqr,d22pqr
real(dp)               :: c,s,cc,cs,ss,c2h,p,q,r,lp,lq,L
integer(spi)           :: i
call eigensym2(x,vv,oo)
if(oo(1,1)<=u0 .or. oo(2,2)<=u0)stop 'In logsym2; matrix x is not positive definite'
c=vv(1,1); s=vv(1,2); cc=c*c; cs=c*s; ss=s*s; c2h=(cc-ss)*o2
p=u1/oo(1,1); q=u1/oo(2,2); lp=log(p); lq=log(q); oo(1,1)=-lp; oo(2,2)=-lq
z=matmul(vv,matmul(oo,transpose(vv)))
l=(lp-lq)*o2; r=sqrt(p*q)/sinhox(l)
d11(1,:)=(/ cc,cs /); d11(2,:)=(/ cs,ss/)
d12(1,:)=(/-cs,c2h/); d12(2,:)=(/c2h,cs/)
d22(1,:)=(/ ss,-cs/); d22(2,:)=(/-cs,cc/)
pqr(1,:)=(/p,r/)    ; pqr(2,:)=(/r,q/)
d11pqr=d11*pqr
d12pqr=d12*pqr
d22pqr=d22*pqr
zd(:,:,1,1)=matmul(vv,matmul(d11pqr,transpose(vv)))
zd(:,:,1,2)=matmul(vv,matmul(d12pqr,transpose(vv)))
zd(:,:,2,2)=matmul(vv,matmul(d22pqr,transpose(vv)))
zd(:,:,2,1)=zd(:,:,1,2)
end subroutine logsym2d
 
subroutine id2222(em)!                                                [id2222]
implicit none
real(dp),dimension(2,2,2,2),intent(out):: em
real(dp),dimension(2,2,2,2)            :: id
!data id/u1,u0,u0,u0, u0,o2,o2,u0,  u0,o2,o2,u0, u0,u0,u0,u1/! Effective identity
em=id
end subroutine id2222
 
subroutine chol2(s,c,ff)!                                            [chol2]
use pietc, only: u0
implicit none
real(dp),dimension(2,2),intent(in ):: s
real(dp),dimension(2,2),intent(out):: c
logical                ,intent(out):: ff
real(dp):: r
ff=s(1,1)<=u0; if(ff)return
c(1,1)=sqrt(s(1,1))
c(1,2)=u0
c(2,1)=s(2,1)/c(1,1)
r=s(2,2)-c(2,1)**2
ff=r<=u0;      if(ff)return
c(2,2)=sqrt(r)
end subroutine chol2
 
end module psym2