******************************************************************************** * * * 2DHF version 1-2003 * * Copyright (C) 1996 Jacek Kobus, Leif Laaksonen, Dage Sundholm * * * * This software may be used and distributed according to the terms * * of the GNU General Public License, see README and COPYING. * * * ******************************************************************************** c ### relexch1 ### c subroutine relexch1 (iorb,cw_sor,psi,pot,excp,bpot,d,e,f3,g, & lhs,rhs,wk2) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) real*8 lhs integer*4 cw_sor c_mpi include 'mpif.h' include 'commons_mpi.inc' c_mpi include 'commons8.inc' dimension psi(*),pot(*),excp(*),cw_sor(*), & rhs(*),bpot(*),d(*),e(*),f3(*), & g(*),lhs(*),wk2(*) c integer*4 STATUS(MPI_STATUS_SIZE),STATUS1(MPI_STATUS_SIZE) dimension IREQARR(100),ISTATARR(MPI_STATUS_SIZE) parameter(iprint=0) call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr ) call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr ) if (nel.lt.2 ) return irequest=0 ngorb=i1si(iorb) c iad=0 k=0 200 continue k=k+1 if (k.le.nexchorb(iorb)) then iorb1=load(iorb,myid,k) c if (iorb1.ne.0) then if (iorb1.eq.0) goto 200 iswtch=0 c orbitals in increasing order if (iorb.lt.iorb1) then in1=iorb in2=iorb1 else in1=iorb1 in2=iorb endif ib1=i1b(in1) ng1=i1si(in1) ib2=i1b(in2) ng2=i1si(in2) ipc=in1+in2*(in2-1)/2 iwexch(ipc)=-1 ibexp=i3b(ipc) ngexp=i3si(ipc) idel=iabs(mgx(6,in1)-mgx(6,in2)) if (in1.eq.in2) idel=2*mgx(6,in1) idel1=mgx(6,in1)+mgx(6,in2) c prepare right-hand side of the poisson equation do i=1,ngorb wk2(i)=0.d0 enddo ngrid=min(ng1,ng2) call prod2 (ngrid,psi(ib1),psi(ib2),wk2) call prod2 (ngorb,wk2,g,rhs) ipex=idel-2*(idel/2) c if (mod(idel,2).eq.0) then if (ipex.eq.0) then isym= 1 else isym=-1 endif 1000 continue if (iprint.eq.1) then write (*,'("relexech1: myid,k,iorb,iorb1,ipc,ibexp",6i8)') & myid,k,iorb,iorb1,ipc,ibexp endif do itr1=1,maxsor1 do ig=1,i1ng(iorb) do i=1,i1si(iorb) wk2(i)=0.d0 enddo omega=ovfexch(ig) omega1=1.d0-omega ioffst=ioffs(ig) ioffs1=ioffst+1 ibexp1=ibexp+ioffst do i=1,4 dmu2t(i)=dmu2(i,ig) dmu1t(i)=dmu1(i,ig) enddo c prepare left-hand side of the poisson equation c include the diagonal part of the differentiation c operator in lhs if (idel.eq.0) then do i=1,ngsize(ig) lhs(ioffst+i)=f3(ioffst+i)+diag(ig) enddo else xdel=dble(idel*idel) do i=1,ngsize(ig) lhs(ioffst+i)=f3(ioffst+i)+ & xdel*e(ioffst+i)+diag(ig) enddo endif if (i1ng(iorb).eq.1) then c icase=1 ifill=1 ngrd1 =ingr1(1,ig) ngrd6a=ingr1(2,ig) ngrd6b=ingr1(3,ig) ngrd7 =ingr1(4,ig) call fill (nni,nmu(ig),isym,excp(ibexp),wk2) do itr2=1,maxsor2/lsor call sor & (wk2,lhs(ioffs1),rhs(ioffs1), & bpot(ioffs1),d(ioffs1), & cw_sor(iadext(ig)),cw_sor(iadnor(ig)), & cw_sor(iadex1(ig)),cw_sor(iadex2(ig)), & cw_sor(iadex3(ig))) enddo call refill (nni,nmu(ig),excp(ibexp),wk2) else if (ig.eq.1) then c icase=2 ifill=2 ngrd1 =ingr2(1,ig) ngrd6a=ingr2(2,ig) ngrd6b=ingr2(3,ig) ngrd7 =ingr2(4,ig) elseif (ig.ne.1.and.ig.ne.i1ng(iorb)) then c icase=2 ifill=3 ngrd1 =ingr2(1,ig) ngrd6a=ingr2(2,ig) ngrd6b=ingr2(3,ig) ngrd7 =ingr2(4,ig) elseif (ig.eq.i1ng(iorb)) then ifill=4 ngrd1 =ingr1(1,ig) ngrd6a=ingr1(2,ig) ngrd6b=ingr1(3,ig) ngrd7 =ingr1(4,ig) muoffs=iemu(ig-1)-1 endif endif enddo c if (myid.ne.0) then c irequest=ipc c irequest=irequest+1 c call MPI_ISEND(excp(ibexp),ngexp, c & MPI_DOUBLE_PRECISION, c & 0,ipc,MPI_COMM_WORLD,IREQARR(irequest),IERROR1) c write(*,'("isend",5i8,2d16.6)'),myid,irequest,iorb,iorb1, c & ibexp,excp(ibexp),excp(ibexp+100) c endif enddo if (iswtch.eq.1) goto 10 if (ilc(ipc).lt.2) goto 10 iwexch(ipc)=-2 idel=idel1 ibexp=ibexp+ngexp ipc=ipc+norb*(norb+1)/2 k=k+1 iswtch=1 goto 1000 10 continue goto 200 endif 210 continue c stop irequest=0 c send calculated values to 0 process if (myid.eq.0) then do npr=1,numprocs-1 do kk=1,nexchorb(iorb) iorb1=load(iorb,npr,kk) if (iorb1.ne.0) then ibexp=iexchadd(iorb,iorb1) ngexp=iexchlen(iorb,iorb1) iswtch=0 if (iorb.lt.iorb1) then in1=iorb in2=iorb1 else in1=iorb1 in2=iorb endif ipc=in1+in2*(in2-1)/2 irequest=irequest+1 26 continue call MPI_IRECV(excp(ibexp),ngexp, & MPI_DOUBLE_PRECISION, c & myid,ipc,MPI_COMM_WORLD,IREQARR(irequest), & npr,MPI_ANY_TAG,MPI_COMM_WORLD,IREQARR(irequest), & ISTATARR,IERROR1) call MPI_WAIT(IREQARR(irequest),ISTATARR,IERROR) if (iprint.eq.2) then write(*,'("irecv",5i8,2d16.6)'),myid,irequest, & iorb,iorb1,ibexp,excp(ibexp),excp(ibexp+100) endif if (iswtch.eq.1) goto 25 if (ilc(ipc).lt.2) goto 25 ibexp=ibexp+ngexp ipc=ipc+norb*(norb+1)/2 iswtch=1 goto 26 endif 25 continue enddo enddo else c if (iprint.ne.0) then c write(*,'("relexch1: myid,nexchorb(iorb)",4i8)'), c & myid,nexchorb(iorb) c endif do kk=1,nexchorb(iorb) iorb1=load(iorb,myid,kk) if (iprint.eq.2) then write(*,'("relexch1: myid,kk,iorb1",4i8)'), & myid,kk,iorb1 endif if (iorb1.ne.0) then ibexp=iexchadd(iorb,iorb1) ngexp=iexchlen(iorb,iorb1) iswtch=0 if (iorb.lt.iorb1) then in1=iorb in2=iorb1 else in1=iorb1 in2=iorb endif ipc=in1+in2*(in2-1)/2 16 continue irequest=irequest+1 if (iprint.eq.2) then write(*,'("isend",4i8)'),myid,irequest,ipc endif call MPI_ISEND(excp(ibexp),ngexp, & MPI_DOUBLE_PRECISION, & 0,ipc,MPI_COMM_WORLD,IREQARR(irequest),IERROR1) call MPI_WAIT(IREQARR(irequest),ISTATARR,IERROR) if (iprint.eq.2) then c write(*,'("isend",4i8)'),myid,irequest,ipc,istatus endif if (iswtch.eq.1) goto 15 if (ilc(ipc).lt.2) goto 15 ibexp=ibexp+ngexp ipc=ipc+norb*(norb+1)/2 iswtch=1 goto 16 endif 15 continue enddo endif call MPI_BARRIER( MPI_COMM_WORLD,ierr ) 1021 format(/,i5,10d12.3,/(5x,10d12.3)) return end