subroutine lanczos_arpack_c(MatVec,eval,evec,Nblock,Nitermax,bmat,v0,tol,iverbose)
!Interface to Matrix-Vector routine:
interface
subroutine MatVec(Nloc,vin,vout)
integer :: Nloc
complex(8),dimension(Nloc) :: vin
complex(8),dimension(Nloc) :: vout
end subroutine MatVec
end interface
!Arguments
real(8) :: eval(:)![Neigen]
complex(8) :: evec(:,:)![Ns,Neigen]
integer,optional :: Nblock
integer,optional :: Nitermax
! character(len=2),optional :: which
character(len=1),optional :: bmat
complex(8),optional :: v0(size(evec,1))
real(8),optional :: tol
logical,optional :: iverbose
!Dimensions:
integer :: Ns
integer :: Neigen
integer :: maxn,maxnev,maxncv,ldv
integer :: n,nconv,ncv,nev
!Arrays:
complex(8),allocatable :: ax(:),d(:)
complex(8),allocatable :: v(:,:)
complex(8),allocatable :: workl(:),workd(:),workev(:)
complex(8),allocatable :: resid(:)
real(8),allocatable :: rwork(:),rd(:,:)
logical,allocatable :: select(:)
integer :: iparam(11)
integer :: ipntr(14)
!Control Vars:
integer :: ido,ierr,info,ishfts,j,lworkl,maxitr,mode1
logical :: rvec,verb
integer :: i
real(8) :: sigma
real(8) :: tol_
character :: bmat_
character(len=2) :: which_
real(8),external :: dznrm2,dlapy2
real(8),allocatable :: reV(:),imV(:)
integer,allocatable :: Eorder(:)
integer :: logfil, ndigit, mgetv0,&
msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
common /debug/logfil, ndigit, mgetv0,&
msaupd, msaup2, msaitr, mseigt, msapps, msgets, mseupd,&
mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, mneupd,&
mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd
!
Ns = size(evec,1)
Neigen = size(eval)
call assert_shape(Evec,[Ns,Neigen],"Arpack_c","Evec")
!
maxn = Ns
maxnev = Neigen
maxncv = 10*Neigen ; if(present(Nblock))maxncv = Nblock
maxitr = 512 ; if(present(Nitermax))maxitr = Nitermax
bmat_ = 'I' ; if(present(bmat))bmat_=bmat
which_='SR' ! ; if(present(which))which_=which
tol_ = 0d0 ; if(present(tol))tol_=tol
verb =.false. ; if(present(iverbose))verb=iverbose
!
if(bmat_/='I'.AND.bmat_/='G')stop "ARPACK: selected *bmat* is wrong. can be [I, G]"
! which_/="LM" .OR. &
! which_/="SM" .OR. &
! which_/="LI" .OR. &
! which_/="SI"
! if(which_/="LR" .OR. which_/="SR")then
! stop "ARPACK: selected *which_* is wrong. can be [LM, SM, LR, SR, LI, SI]"
! endif
if(verb)then
ndigit=-4
logfil = 6
mcaupd=1;mnaupd=1
mcaup2=1;mnaup2=1
mceupd=4;mneupd=4
endif
!
ldv = Ns ; if(maxncv>Ns)maxncv=Ns
n = maxn
nev = maxnev
ncv = maxncv
!
allocate(ax(n))
allocate(d(ncv))
allocate(resid(n))
allocate(v(ldv,ncv))
allocate(workd(3*n))
allocate(workev(3*ncv))
allocate(workl(ncv*(3*ncv+5) + 10))
allocate(rwork(ncv))
allocate(rd(ncv,3))
allocate(select(ncv))
ax = zero
d = zero
v = zero
workl = zero
workd = zero
resid = zero
workev = zero
rwork = 0d0
rd = 0d0
lworkl = ncv*(3*ncv+5) + 10
info = 1
ido = 0
ishfts = 1
iparam(1) = ishfts
iparam(3) = maxitr
mode1 = 1
iparam(7) = mode1
if(present(v0))then
resid=v0
else
call mt_random(resid)
endif
resid=resid/sqrt(dot_product(resid,resid))
!
!MAIN LOOP, REVERSE COMMUNICATION
do
call znaupd(ido,bmat_,n,which_,nev,tol_,resid,ncv,v,ldv,&
iparam,ipntr,workd,workl,lworkl,rwork,info)
if(ido/=-1.AND.ido/=1)exit
! Perform matrix vector multiplication:
! y <--- OP*x ; workd(ipntr(1))=input, workd(ipntr(2))=output
call MatVec(N,workd(ipntr(1)),workd(ipntr(2)) )
end do
!
!POST PROCESSING:
if(info/=0)then
write(*,'(a,i6)')'Warning/Error in ZNAUPD, info = ', info
include "error_msg_arpack.h90"
else
rvec = .true.
call zneupd (rvec,'A',select,d,v,ldv,sigma,workev,bmat_,n,which_,&
nev,tol_,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,rwork,ierr)
! Eigenvalues are returned in the first column of the two dimensional
! array D and the corresponding eigenvectors are returned in the first
! NCONV (=IPARAM(5)) columns of the two dimensional array V if requested.
! Otherwise, an orthogonal basis for the invariant subspace corresponding
! to the eigenvalues in D is returned in V.
do j=1,neigen
eval(j)=dreal(d(j))
end do
allocate(Eorder(Neigen))
call sort_array(Eval,Eorder)
do j=1,neigen
evec(:,j)=v(:,Eorder(j))
enddo
!
! Compute the residual norm || A*x - lambda*x ||
! for the NCONV accurately computed eigenvalues and eigenvectors.
if(ierr/=0)then
write(*,'(a,i6)')'Error with ZNEUPD, IERR = ',ierr
write(*,'(a)')'Check the documentation of ZNEUPD.'
stop
else
nconv = iparam(5)
end if
!
include "info_msg_arpack.h90"
!
!if(nconv == 0) stop "ARPACK:: no converged eigenvalues have been found."
!
endif
deallocate(ax,d,resid,v,workd,workev,workl,rwork,rd,select)
!
contains
!+------------------------------------------------------------------+
!PURPOSE : Sort an array, gives the new ordering of the label.
!+------------------------------------------------------------------+
subroutine sort_array(array,order2,no_touch_array)
implicit none
real(8),dimension(:) :: array
real(8),dimension(size(array)) :: backup
integer,dimension(size(array)) :: order
integer,dimension(size(array)),optional :: order2
integer :: i
logical,optional :: no_touch_array
do i=1,size(order)
order(i)=i
enddo
call qsort_sort( array, order, 1, size(array) )
if(.not.present(no_touch_array))then
do i=1,size(order)
backup(i)=array(order(i))
enddo
array=backup
endif
if(present(order2)) order2=order
end subroutine sort_array
!---------------------------------------------!
recursive subroutine qsort_sort( array, order, left, right )
implicit none
real(8), dimension(:) :: array
integer, dimension(:) :: order
integer :: left
integer :: right
integer :: i
integer :: last
if ( left .ge. right ) return
call qsort_swap( order, left, qsort_rand(left,right) )
last = left
do i = left+1, right
if ( compare(array(order(i)), array(order(left)) ) .lt. 0 ) then
last = last + 1
call qsort_swap( order, last, i )
endif
enddo
call qsort_swap( order, left, last )
call qsort_sort( array, order, left, last-1 )
call qsort_sort( array, order, last+1, right )
end subroutine qsort_sort
!---------------------------------------------!
subroutine qsort_swap( order, first, second )
implicit none
integer, dimension(:) :: order
integer :: first, second
integer :: tmp
tmp = order(first)
order(first) = order(second)
order(second) = tmp
end subroutine qsort_swap
!---------------------------------------------!
integer function qsort_rand( lower, upper )
implicit none
integer :: lower, upper
real(4) :: r
r=drand()
qsort_rand = lower + nint(r * (upper-lower))
end function qsort_rand
!---------------------------------------------!
function drand()
implicit none
real(8) :: drand
real(4) :: r
call random_number(r)
drand=dble(r)
end function drand
!---------------------------------------------!
function compare(f,g)
implicit none
real(8) :: f,g
integer :: compare
if(f<g) then
compare=-1
else
compare=1
endif
end function compare
end subroutine lanczos_arpack_c