exponential of the identity

[JP]

Hi,

I am using expokit and the routine ZGEXPV. My problem can be
simplified to the following : when I call ZGEXPV to calculate the
exponential of the identity matrix, the imaginary part of the result
is not correct??
This is my test program :

*----------------------------------------------------------------------
*
* The matrix-vector multiplication routine.
*
* y:=A*x
*
* In this example, the matrix is the identity!
*
*----------------------------------------------------------------------
*
subroutine zprod( x, y )
implicit none
complex*16 x(*), y(*)

integer n, j
common /CMAT/ n

do j = 1,n
y(j) = x(j)
enddo
END
*----------------------------------------------------------------------|
program expoIdentity
implicit none
external zprod

integer n,adim
parameter( n = 500)
common /CMAT/ adim

*--- arguments variables ...
integer m, lwsp, liwsp
parameter( m = 50 )
parameter( lwsp = n*(m+2)+5*(m+2)**2+7, liwsp = m+2 )
integer iwsp(liwsp)
double precision t, tol, anorm
complex*16 v(n), w(n), wsp(lwsp)

integer i, itrace, iflag

*--- the norm of the matrix A (A is the identity!)

anorm=1.d0

*--- the operand vector v is set to 1
do i = 1,n
v(i) = (1.0d0,0.d0)
enddo

*--- set other input arguments ...
t = 1.0d0
tol = 1.0d-7

adim=n
itrace = 0

*--- compute w = exp(t*A)v with ZGEXPV ...
call ZGEXPV( n, m, t,v,w, tol, anorm,
. wsp,lwsp, iwsp,liwsp, zprod, itrace, iflag )

*---
print
9001,'----------------------------------------------------'
print 9001,'ZGEXPV has completed:'
print
9001,'----------------------------------------------------'
print 9001,'w(1:n) ='
do i = 1,1
print*,w(i)
enddo

*--- display some statistics if desired ...
print 9001,'final
report----------------------------------------'
print 9002,'||A||_inf = ',anorm
print 9003,'n =',n
print 9003,'m =',m
print 9003,'itrace =',itrace
print 9003,'iflag =',iflag
print 9003,'ibrkflag =',iwsp(6)
print 9003,'mbrkdwn =',iwsp(7)
print 9003,'nstep =',iwsp(4)
print 9003,'nreject =',iwsp(5)
print 9003,'nmult =',iwsp(1)
print 9003,'nexph =',iwsp(2)
print 9003,'nscale =',iwsp(3)

print 9002,'tol = ',tol
print 9002,'t = ',t
print 9002,'tbrkdwn = ',DBLE( wsp(7) )
print 9002,'step_min = ',DBLE( wsp(1) )
print 9002,'step_max = ',DBLE( wsp(2) )
print 9002,'max_round = ',DBLE( wsp(3) )
print 9002,'sum_round = ',DBLE( wsp(4) )
print 9002,'max_error = ',DBLE( wsp(5) )
print 9002,'sum_error = ',DBLE( wsp(6) )
print 9002,'hump = ',DBLE( wsp(9) )
print 9002,'scale-norm= ',DBLE( wsp(10) )


9001 format(A)
9002 format(A,E8.2)
9003 format(A,I9)
END

The results :

happy breakdown: mbrkdwn = 1 h = 6.206335383118183E-016
----------------------------------------------------
ZGEXPV has completed:
----------------------------------------------------
w(1:n) =
(2.71828181040259,0.199611686544197)
final report----------------------------------------
||A||_inf = 0.10E+01
n = 500
m = 50
itrace = 0
iflag = 0
ibrkflag = 1
mbrkdwn = 1
nstep = 1
nreject = 0
nmult = 1
nexph = 1
nscale = 2
tol = 0.10E-06
t = 0.10E+01
tbrkdwn = 0.00E+00
step_min = 0.10E+01
step_max = 0.10E+01
max_round = 0.00E+00
sum_round = 0.00E+00
max_error = 0.10E-06
sum_error = 0.10E-06
hump = 0.27E+01
scale-norm= 0.27E+01

I would like to have : w(1)= (2.71828181040259,0.d0)??? Has anyone an
idea why I get this result??

Best regards

JP