! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! $ ! $ ! 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 |