"""
Fortran module ExpokitOps:Uses the Expokit package to calculate the exponential of a
non-hermitian matrix.
**Authors**
* M. L. Wall
* D. Jaschke
* Roger B. Sidje (rbs@maths.uq.edu.au), https://www.maths.uq.edu.au/expokit/,
* author of the original Expokit package
**Details**
The copyrights of the Expokit package apply to this module. A
copy of the Expokit copyright has been delived within OpenMPS. More info:
EXPOKIT: Software Package for Computing Matrix Exponentials.
ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
The following procedures are available.
+---------------------------+-------------+-------------+
| procedure | include.f90 | mpi.f90 |
+===========================+=============+=============+
| exp_kit | X | |
+---------------------------+-------------+-------------+
"""
[docs]def ZGPADM():
"""
fortran-subroutine - ?? ()
Routine from Expokit for exponentiating non-Hermitian effective
Hamiltonians.
**Arguments**
ideg : INTEGER, in
The degre of the diagonal Pade to be used. A value of 6 is
generally satisfactory.
m : INTEGER, in
Order of H.
t : REAL, in
time-scale (can be < 0).
H : COMPLEX(ldh, m), in
Argument matrix to be exponentiated.
ldh : INTEGER, in
Leading order (first dimension/number of rows) of the matrix H.
wsp : COMPLEX(lwsp), out
Workspace with lwsp .ge. 4*m*m+ideg+1.
lwsp : INTEGER, in
Size of the workspace wsp.
ipiv : INTEGER(m), out
Workspace.
iexph : INTEGER, out
Number such that wsp(iexph) points to exp(tH)
i.e., exp(tH) is located at wsp(iexph ... iexph+m*m-1)
NOTE: if the routine was called with wsp(iptr),
then exp(tH) will start at wsp(iptr+iexph-1).
ns : INTEGER, out
Number of scaling-squaring used.
iflag : INTEGER, out
Exit flag: iflag = 0 - no problem; iflag < 0 - problem
**Details**
Computes exp(t*H), the matrix exponential of a general complex
matrix in full, using the irreducible rational Pade approximation
to the exponential exp(z) = r(z) = (+/-)( I + 2*(q(z)/p(z)) ),
combined with scaling-and-squaring.
Roger B. Sidje (rbs@maths.uq.edu.au)
EXPOKIT: Software Package for Computing Matrix Exponentials.
ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine ZGPADM(ideg, m, t, H, ldh, wsp, lwsp, ipiv, iexph, ns, iflag)
implicit none
real(KIND=rKind) :: t
integer :: ideg, m, ldh, lwsp, iexph, ns, iflag, ipiv(m)
complex(KIND=rKind) :: H(ldh,m), wsp(lwsp)
! Local variables
! ---------------
integer :: i,j,k,icoef,mm,ih2,iodd,iused,ifree,iq,ip,iput,iget
real(KIND=rKind) :: hnorm
complex(KIND=rKind) :: cp, cq, scale, scale2, ZERO, ONE
parameter( ZERO=(0.0_rKind,0.0_rKind), ONE=(1.0_rKind,0.0_rKind) )
!intrinsic ABS, CMPLX, DBLE, INT, LOG, MAX
!--- check restrictions on input parameters ...
mm = m*m
iflag = 0
if ( ldh.lt.m ) iflag = -1
if ( lwsp.lt.4*mm+ideg+1 ) iflag = -2
if ( iflag.ne.0 ) stop 'bad sizes (in input of ZGPADM)'
!
!--- initialise pointers ...
!
icoef = 1
ih2 = icoef + (ideg+1)
ip = ih2 + mm
iq = ip + mm
ifree = iq + mm
!
!--- scaling: seek ns such that ||t*H/2^ns|| < 1/2;
! and set scale = t/2^ns ...
!
do i = 1,m
wsp(i) = ZERO
end do
do j = 1,m
do i = 1,m
wsp(i) = wsp(i) + abs( H(i,j) )
end do
end do
hnorm = 0.0_rKind
do i = 1,m
hnorm = max( hnorm,dble(wsp(i)) )
end do
hnorm = abs( t*hnorm )
if ( hnorm.eq.0.0_rKind ) stop 'Error - null H in input of ZGPADM.'
ns = max( 0,int(log(hnorm)/log(2.0_rKind))+2 )
scale = cmplx( t / dble(2.0_rKind**ns), 0.0_rKind, KIND=rKind )
scale2 = scale*scale
!
!--- compute Pade coefficients ...
!
i = ideg+1
j = 2*ideg+1
wsp(icoef) = ONE
do k = 1,ideg
wsp(icoef+k) = (wsp(icoef+k-1)*dble( i-k ))/dble( k*(j-k) )
end do
!
!--- H2 = scale2*H*H ...
!
call ZGEMM( 'n','n',m,m,m,scale2,H,ldh,H,ldh,ZERO,wsp(ih2),m )
!
!--- initialise p (numerator) and q (denominator) ...
!
cp = wsp(icoef+ideg-1)
cq = wsp(icoef+ideg)
do j = 1,m
do i = 1,m
wsp(ip + (j-1)*m + i-1) = ZERO
wsp(iq + (j-1)*m + i-1) = ZERO
end do
wsp(ip + (j-1)*(m+1)) = cp
wsp(iq + (j-1)*(m+1)) = cq
end do
!
!--- Apply Horner rule ...
!
iodd = 1
k = ideg - 1
do while(k.gt.0)
iused = iodd*iq + (1-iodd)*ip
call ZGEMM( 'n','n',m,m,m, ONE,wsp(iused),m,&
wsp(ih2),m, ZERO,wsp(ifree),m )
do j = 1,m
wsp(ifree+(j-1)*(m+1)) = wsp(ifree+(j-1)*(m+1))+wsp(icoef+k-1)
end do
ip = (1-iodd)*ifree + iodd*ip
iq = iodd*ifree + (1-iodd)*iq
ifree = iused
iodd = 1-iodd
k = k-1
end do
!
!--- Obtain (+/-)(I + 2*(p\q)) ...
!
if ( iodd.ne.0 ) then
call ZGEMM( 'n','n',m,m,m, scale,wsp(iq),m,&
H,ldh, ZERO,wsp(ifree),m )
iq = ifree
else
call ZGEMM( 'n','n',m,m,m, scale,wsp(ip),m,&
H,ldh, ZERO,wsp(ifree),m )
ip = ifree
endif
call ZAXPY( mm, -ONE,wsp(ip),1, wsp(iq),1 )
call ZGESV( m,m, wsp(iq),m, ipiv, wsp(ip),m, iflag )
if ( iflag.ne.0 ) stop 'Problem in ZGESV (within ZGPADM)'
call ZDSCAL( mm, 2.0_rKind, wsp(ip), 1 )
do j = 1,m
wsp(ip+(j-1)*(m+1)) = wsp(ip+(j-1)*(m+1)) + ONE
enddo
iput = ip
if ( ns.eq.0 .and. iodd.ne.0 ) then
call ZDSCAL( mm, -1.0_rKind, wsp(ip), 1 )
iexph = iput
return
endif
!
!-- squaring : exp(t*H) = (exp(t*H))^(2^ns) ...
!
iodd = 1
do k = 1,ns
iget = iodd*ip + (1-iodd)*iq
iput = (1-iodd)*ip + iodd*iq
call ZGEMM( 'n','n',m,m,m, ONE,wsp(iget),m, wsp(iget),m,&
ZERO,wsp(iput),m )
iodd = 1-iodd
end do
iexph = iput
end subroutine ZGPADM
"""
return
[docs]def exp_kit_complex():
"""
fortran-subroutine - August 2017 (dj, updated)
Calculate the matrix exponential of a non-hermitian matrix via EXPOKIT.
**Arguments**
U : COMPLEX(\*, \*), out
The exponential of H times the scalar t.
t : REAL, in
Scaling for the matrix H, e.g. size of the time step.
H : COMPLEX(\*, \*), in
H is the matrix to be exponentiated, weighted with a scalar t.
dim : INTEGER, in
Dimension of the square matrix.
**Source Code**
.. hidden-code-block:: fortran
:label: show / hide f90 code
subroutine exp_kit_complex(U, t, H, dim)
integer, intent(in) :: dim
complex(KIND=rKind), dimension(dim, dim), intent(out) :: U
real(KIND=rKind), intent(in) :: t
complex(KIND=rKind), dimension(dim, dim), intent(in) :: H
! Local variables
! ---------------
integer :: m, iflag, ns, iexph, lwsp
complex(KIND=rKind), ALLOCATABLE :: wsp(:)
integer, ALLOCATABLE :: ipiv(:)
real(KIND=rKind) :: nrm
integer :: i, j, expokit_degree
expokit_degree = 6
nrm = 0.0_rKind
do i = 1,size(H, 1)
do j = 1, size(H, 2)
nrm = nrm + abs(H(i, j))
end do
end do
if(nrm .eq. 0.0_rKind) then
U = 0.0_rKind
do i = 1, size(H, 1)
U(i,i) = 1.0_rKind
end do
return
end if
m = size(H, 1)
lwsp = 4 * m * m + ExpoKit_degree + 1
allocate(wsp(lwsp), ipiv(m))
call ZGPADM(ExpoKit_degree, m, t, H, m, wsp, lwsp, ipiv, iexph, ns, &
iflag)
if(iflag .lt. 0) then
print *, 'error from expokit:', iflag
stop
end if
U = reshape(wsp(iexph:iexph + m * m - 1), (/m, m/))
deallocate(wsp, ipiv)
end subroutine exp_kit_complex
"""
return