subroutine lanczos_parpack_c(MpiComm,MatVec,eval,evec,Nblock,Nitermax,v0,tol,iverbose,vrandom) !Arguments integer :: MpiComm !Interface to Matrix-Vector routine: interface subroutine MatVec(nchunk,vin,vout) integer :: nchunk complex(8),dimension(nchunk) :: vin,vout end subroutine MatVec end interface !Arguments real(8) :: eval(:) ![Neigen] complex(8) :: evec(:,:) ![Nloc,Neigen] integer,optional :: Nblock integer,optional :: Nitermax ! character(len=2),optional :: which complex(8),optional :: v0(size(evec,1)) real(8),optional :: tol logical,optional :: iverbose logical,optional :: vrandom !Dimensions: integer :: Ns integer :: Neigen integer :: maxn,maxnev,maxncv,ldv,nloc integer :: n,nconv,ncv,nev !Arrays: complex(8),allocatable :: evec_tmp(:) ![Nloc] see above complex(8),allocatable :: ax(:),d(:) complex(8),allocatable :: v(:,:) complex(8),allocatable :: workl(:),workd(:),workev(:) complex(8),allocatable :: resid(:),vec(:) 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 :: verb,vran integer :: i real(8) :: sigma,norm,norm_tmp real(8) :: tol_ character :: bmat character(len=2) :: which_ real(8),external :: dznrm2,dlapy2 real(8),allocatable :: reV(:),imV(:) integer,allocatable :: Eorder(:) !MPI logical :: mpi_master 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 ! if(MpiComm == MPI_COMM_NULL)return ! !MPI setup: mpi_master= Get_master_MPI(MpiComm) ! Ns = size(evec,1) Neigen = size(eval) call assert_shape(Evec,[Ns,Neigen],"P_arpack_d","Evec") ! maxn = Ns maxnev = Neigen ! maxncv = 10*Neigen ; if(present(Nblock))maxncv = Nblock maxitr = 512 ; if(present(Nitermax))maxitr = Nitermax which_ = 'SR' !; if(present(which))which_=which tol_ = 0d0 ; if(present(tol))tol_=tol verb = .false. ; if(present(iverbose))verb=iverbose vran = .true. ; if(present(vrandom))vran=vrandom if(verb)then ndigit=-4 logfil = 6 mcaupd=1;mnaupd=1 mcaup2=1;mnaup2=1 mceupd=4;mneupd=4 endif ! ! if(maxncv>Ns)then maxncv=Ns print*,"PARPACK WARNING Ncv > Ns: reset block size to ",Ns !BUG FIX FOR THE BLOCK RESIZE STUCK BEHAVIOR, from Xuanyu Long call MPI_ALLREDUCE(MPI_IN_PLACE,maxncv,1,MPI_INTEGER,MPI_MIN,MpiComm,ierr) endif ! ! n = maxn nev = maxnev ncv = maxncv bmat = 'I' ! !========================================================================= ! Setup distribution of data to nodes: ldv = size(evec,1) !ldv is the SMALL dimension ! if(ldv>0.AND.ldv < ncv)stop "LANCZOS_PARPACK_C error: ldv < maxNblock. Unstable! Find a way to increase ldv... (less cpu?)" if ( ldv > maxn ) then stop ' ERROR with _SDRV1: NLOC is greater than MAXNLOC ' else if ( nev > maxnev ) then stop ' ERROR with _SDRV1: NEV is greater than MAXNEV ' else if ( ncv > maxncv ) then stop ' ERROR with _SDRV1: NCV is greater than MAXNCV ' end if ! allocate(ax(ldv)) allocate(d(ncv)) allocate(resid(ldv)) allocate(v(ldv,ncv)) allocate(workd(3*ldv)) 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 = zero rd = zero lworkl = ncv*(3*ncv+5) + 10 info = 1 ido = 0 select = .true. ! ! if(present(v0))then resid = v0 else if(vran)then call mt_random(resid) else resid = one !start with unitary vector 1/sqrt(Ndim) endif norm_tmp = dot_product(resid,resid) norm = 0d0 call AllReduce_MPI(MpiComm,norm_tmp,norm) resid=resid/sqrt(norm) endif ! ishfts = 1 mode1 = 1 iparam(1) = ishfts iparam(3) = maxitr iparam(7) = mode1 ! ! ! MAIN LOOP (Reverse communication loop) do call pznaupd(MpiComm,ido,bmat,ldv,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(ldv,workd(ipntr(1)),workd(ipntr(2))) end do ! !POST PROCESSING: if(info/=0)then if(mpi_master)then write(*,'(a,i6)')'Warning/Error in PZNAUPD, info = ', info include "error_msg_arpack.h90" endif else ! call pzneupd (MpiComm,.true.,'All',select,d,v,ldv,sigma,workev,bmat,& ldv,which_,nev,tol_,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,rwork,ierr) ! do j=1,neigen eval(j)=dreal(d(j)) enddo allocate(Eorder(Neigen)) call sort_array(Eval,Eorder) ! evec=zero do j=1,neigen do i=1,ldv evec(i,j)=v(i,Eorder(j)) enddo 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 PZNEUPD, IERR = ',ierr write(*,'(a)')'Check the documentation of PZNEUPD.' else nconv = iparam(5) end if ! if(mpi_master)then include "info_msg_arpack.h90" end if ! ! ! if(mpi_master.and.nconv==0.and.verb)stop "None of the required values was 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_parpack_c