ExpokitOps

Fortran module ExpokitOps:Uses the Expokit package to calculate the exponential of a non-hermitian matrix.

Authors

• 13. 12. Wall

• 4. 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

ExpokitOps_f90.ZGPADM()

fortran-subroutine - ?? () Routine from Expokit for exponentiating non-Hermitian effective Hamiltonians.

Arguments

idegINTEGER, in

The degre of the diagonal Pade to be used. A value of 6 is generally satisfactory.

mINTEGER, in

Order of H.

tREAL, in

time-scale (can be < 0).

HCOMPLEX(ldh, m), in

Argument matrix to be exponentiated.

ldhINTEGER, in

Leading order (first dimension/number of rows) of the matrix H.

wspCOMPLEX(lwsp), out

Workspace with lwsp .ge. 4*m*m+ideg+1.

lwspINTEGER, in

Size of the workspace wsp.

ipivINTEGER(m), out

Workspace.

iexphINTEGER, 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).

nsINTEGER, out

Number of scaling-squaring used.

iflagINTEGER, 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

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

ExpokitOps_f90.exp_kit_complex()

fortran-subroutine - August 2017 (dj, updated) Calculate the matrix exponential of a non-hermitian matrix via EXPOKIT.

Arguments

UCOMPLEX(*, *), out

The exponential of H times the scalar t.

tREAL, in

Scaling for the matrix H, e.g. size of the time step.

HCOMPLEX(*, *), in

H is the matrix to be exponentiated, weighted with a scalar t.

dimINTEGER, in

Dimension of the square matrix.

Source Code

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