! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
! $
! $
! subroutine matvec $
! $
! $
! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!==============================================================================
!
! this subroutine calculates the multiplication of a matrix a by a vector x.
! it is suited for the use of expokit zgexpv subroutine. the calculation is the
! following
!
! a*x=y
!
! and it takes into account the sparsity of the matrix a (the same scheme is
! used in the taylor expansion)
!
!==============================================================================
subroutine matvec(x,y,n,ngrid_r,lmin,lmax,l_coupling,h_diag,h_diag1,h_diag2,&
&h_nondiag,h_laser,h_laser1,h_laser2)
implicit none
! input/output variables declaration
integer(8) :: n,ngrid_r,lmin,lmax,l_coupling
real(8) :: h_diag1,h_diag2
real(8), dimension(ngrid_r,0:lmax) :: h_diag,h_laser
real(8), dimension(0:lmax) :: h_laser1,h_laser2
real(8), dimension(ngrid_r,0:lmax,0:lmax) :: h_nondiag
complex(8), dimension(n) :: x,y
complex(8) :: wf(1:ngrid_r,0:lmax)
complex(8),dimension(1:size(wf,1),0:size(wf,2)-1) :: wf1,wf2
! subroutine variables declaration
integer :: i,j,l
y=dcmplx(0.d0,0.d0)
wf=dcmplx(0.d0,0.d0)
wf1=dcmplx(0.d0,0.d0)
wf2=dcmplx(0.d0,0.d0)
!TODO OPENMP HERE
do i=1,ngrid_r
do l=0,lmax
j=i+l*ngrid_r
wf(i,l)=x(j)
end do
end do
wf1=wf
call velmatmul(ngrid_r,wf1,wf2,h_diag,h_diag1,h_diag2,h_nondiag,&
&h_laser,h_laser1,h_laser2,lmin,lmax,l_coupling)
wf=wf2
!TODO OPENMP HERE
do i=1,ngrid_r
do l=0,lmax
j=i+l*ngrid_r
y(j)=-dcmplx(0.d0,1.d0)*wf(i,l)
end do
end do
return
end subroutine matvec
subroutine velmatmul(nr,wf1,wf2,h0matrix,h0matrix1,h0matrix2,h1matrix,&
&h2matrix,h2matrix1,h2matrix2,lmin,lmax,lmaxcoup)
implicit none
integer(8) :: nr,lmin,lmax,i,j,l,l1,l2
real(8) :: h0matrix(1:nr,0:lmax)
real(8) :: h1matrix(1:nr,0:lmax,0:lmax)
real(8) :: h0matrix1,h0matrix2
!complex(8) :: h2matrix(1:nr,0:lmax)
!complex(8) :: h2matrix1(0:lmax)
!complex(8) :: h2matrix2(0:lmax)
real(8) :: h2matrix(1:nr,0:lmax)
real(8) :: h2matrix1(0:lmax)
real(8) :: h2matrix2(0:lmax)
complex(8) :: wf1(1:nr,0:lmax)
complex(8) :: wf2(1:nr,0:lmax)
integer(8) :: n1,n2,n3,lmaxcoup
n1=nr-1
n2=nr-2
n3=nr-3
!TODO OPENMP HERE
do l=lmin,lmax
wf2(1:nr,l)=h0matrix(1:nr,l)*wf1(1:nr,l)
wf2(1:n2,l)=wf2(1:n2,l)+h0matrix2*wf1(3:nr,l)
wf2(1:n1,l)=wf2(1:n1,l)+h0matrix1*wf1(2:nr,l)
wf2(2:nr,l)=wf2(2:nr,l)+h0matrix1*wf1(1:n1,l)
wf2(3:nr,l)=wf2(3:nr,l)+h0matrix2*wf1(1:n2,l)
end do
!TODO OPENMP HERE
do l=lmin,lmax
j=l-1
if(j>=0)then
wf2(1:nr,j)=wf2(1:nr,j)+h2matrix(1:nr,j)*wf1(1:nr,l)
!wf2(1:n2,j)=wf2(1:n2,j)+h2matrix2(j)*wf1(3:nr,l)
!wf2(1:n1,j)=wf2(1:n1,j)+h2matrix1(j)*wf1(2:nr,l)
!wf2(2:nr,j)=wf2(2:nr,j)-h2matrix1(j)*wf1(1:n1,l)
!wf2(3:nr,j)=wf2(3:nr,j)-h2matrix2(j)*wf1(1:n2,l)
end if
j=l+1
if(j<=lmax)then
wf2(1:nr,j)=wf2(1:nr,j)+h2matrix(1:nr,l)*wf1(1:nr,l)
!wf2(1:nr,j)=wf2(1:nr,j)-h2matrix(1:nr,l)*wf1(1:nr,l)
!wf2(1:n2,j)=wf2(1:n2,j)+h2matrix2(l)*wf1(3:nr,l)
!wf2(1:n1,j)=wf2(1:n1,j)+h2matrix1(l)*wf1(2:nr,l)
!wf2(2:nr,j)=wf2(2:nr,j)-h2matrix1(l)*wf1(1:n1,l)
!wf2(3:nr,j)=wf2(3:nr,j)-h2matrix2(l)*wf1(1:n2,l)
end if
end do
!TODO OPENMP HERE
do l1=lmin,lmax
do l2=lmin,lmax
j=abs(l1-l2)
if(0<j.and.j<=lmaxcoup)then
do i=1,nr
wf2(i,l1)=wf2(i,l1)+h1matrix(i,l1,l2)*wf1(i,l2)
!print*,'h1,wf1',h1matrix(i,l1,l2),wf1(i,l2)
end do
!read(*,*)
end if
end do
end do
end subroutine velmatmul
|
