subroutine lanczos_parpack_d(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 real(8),dimension(nchunk) :: vin,vout end subroutine MatVec end interface !Arguments real(8) :: eval(:)![Neigen] real(8) :: evec(:,:)![Ns,Neigen] integer,optional :: Nblock integer,optional :: Nitermax ! character(len=2),optional :: which real(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 integer :: n,nconv,ncv,nev !Arrays: real(8),allocatable :: evec_tmp(:) ![Nloc] see above real(8),allocatable :: ax(:),d(:,:) real(8),allocatable :: resid(:),vec(:) real(8),allocatable :: workl(:),workd(:) real(8),allocatable :: v(:,:) logical,allocatable :: select(:) integer :: iparam(11) integer :: ipntr(11) !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 :: dnrm2 !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_ = 'SA' ! ; 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 ! ldv = maxn n = maxn nev = maxnev ncv = maxncv bmat = 'I' ! !========================================================================= ! Setup distribution of data to nodes: 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(resid(ldv)) allocate(workd(3*ldv)) allocate(v(ldv,ncv)) allocate(d(ncv,2)) allocate(workl(ncv*(ncv+8))) allocate(select(ncv)) ! ax =0.d0 d =0.d0 resid =0.d0 workl =0.d0 workd =0.d0 v =0.d0 lworkl = ncv*(ncv+8)+10 info = 1 ido = 0 select = .true. ! if(present(v0))then resid = v0 else allocate(vec(ldv)) if(vran)then call mt_random(vec) else vec = 1d0 !start with unitary vector 1/sqrt(Ndim) endif norm_tmp = dot_product(vec,vec) norm = 0d0 call AllReduce_MPI(MpiComm,norm_tmp,norm) resid=vec/sqrt(norm) deallocate(vec) endif ! ishfts = 1 mode1 = 1 iparam(1) = ishfts iparam(3) = maxitr iparam(7) = mode1 ! ! MAIN LOOP (Reverse communication loop) do call pdsaupd(MpiComm,ido,bmat,ldv,which_,nev,tol_,resid,& ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,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 PDSAUPD, info = ', info include "error_msg_arpack.h90" endif else ! call pdseupd (MpiComm,.true.,'All',select,d,v,ldv,sigma,bmat,& ldv,which_,nev,tol_,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,ierr) do j=1,neigen eval(j)=d(j,1) enddo ! evec=0d0 do j=1,neigen do i=1,ldv evec(i,j)=v(i,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 PDSEUPD, IERR = ',ierr write(*,'(a)')'Check the documentation of PDSEUPD.' stop else nconv = iparam(5) end if ! if(mpi_master)then include "info_msg_arpack.h90" endif ! ! if(mpi_master.and.nconv==0.and.verb)stop "None of the required values was found." endif deallocate(ax,resid,workd,v,d,workl,select) end subroutine lanczos_parpack_d