Skip to content

Commit

Permalink
Lobatto and grid making gfortran-proof
Browse files Browse the repository at this point in the history 
Lobatto and grid making gfortran-proof
  • Loading branch information
@Trovemaster
Trovemaster committed Jan 13, 2019
1 parent f7ce7fd commit 9e405eb
Show file tree
Hide file tree
Showing 8 changed files with 178 additions and 129 deletions.
206 changes: 105 additions & 101 deletions Lobatto.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,10 @@ module Lobatto

subroutine AllLobattoKE(KE,npnt2)
implicit none
real(rk) :: KE(1:npnt2,1:npnt2),a0,rmin
integer(ik):: bra,ket,Ntot,Ntotp2,npnt2
integer(ik),intent(in) :: npnt2
real(rk),intent(out) :: KE(1:npnt2,1:npnt2)
integer(ik):: bra,ket,Ntot,Ntotp2
real(rk) :: a0,rmin
real(rk),allocatable,dimension(:) :: rpt, w
real(rk), allocatable, dimension(:,:) :: derLobMat
!
Expand Down Expand Up @@ -50,10 +52,10 @@ end subroutine AllLobattoKE
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!
subroutine OlLobattoDVR(result,bra,ket,Ntot)
subroutine OlLobattoDVR(result,bra,ket)
implicit none
real(rk), intent(out):: result
integer(ik), intent(in):: bra, ket, Ntot
integer(ik), intent(in):: bra, ket
!
!
if (bra == ket) then
Expand All @@ -71,18 +73,19 @@ end subroutine OlLobattoDVR

subroutine KeLobattoDVR(result,bra,ket,Ntot,rpt,w,derLobMat)
implicit none
real(rk), intent(out):: result,rpt(0:Ntot+1),w(0:Ntot+1)
integer(ik), intent(in):: bra,ket,Ntot
real(rk), intent(out):: result,rpt(0:Ntot+1),w(0:Ntot+1)
real(rk), intent(in):: derLobMat(0:Ntot+1,0:Ntot+1)
real(rk):: pf, intresult
integer(ik):: k


pf=1d0!/sqrt(w(bra)*w(ket))
intresult=0d0;
do k=0,Ntot+1
intresult=intresult+w(k)*derLobMat(bra,k)*derLobMat(ket,k)
end do
integer(ik):: k
!
rpt = 0
w = 0
pf=1d0!/sqrt(w(bra)*w(ket))
intresult=0d0;
do k=0,Ntot+1
intresult=intresult+w(k)*derLobMat(bra,k)*derLobMat(ket,k)
end do


result=intresult*pf;
Expand Down Expand Up @@ -117,21 +120,21 @@ end subroutine rminus2LobattoDVR

subroutine derLobattoMat(result,Ntot,rpt,w)
implicit none
real(rk), intent(out) :: result(0:Ntot+1,0:Ntot+1)
integer(ik), intent(in) :: Ntot
real(rk), intent(out) :: result(0:Ntot+1,0:Ntot+1)
real(rk), intent(in) :: rpt(0:Ntot+1),w(0:Ntot+1)
real(rk) :: elementres
integer(ik) :: n,eta,k

do n=0,Ntot+1
do eta=n,Ntot+1
call derLobatto(result(n,eta),n,eta,Ntot,rpt,w)
if(n.ne.eta) result(eta,n) = (-1.0d0) * (w(eta)/w(n)) * result(n,eta)
result(n,eta) = result(n,eta) * 1.0d0/sqrt(w(n))
if(n.ne.eta) result(eta,n) = result(eta,n) * 1.0d0/sqrt(w(eta))
end do
end do
return
do n=0,Ntot+1
do eta=n,Ntot+1
call derLobatto(result(n,eta),n,eta,Ntot,rpt,w)
if(n.ne.eta) result(eta,n) = (-1.0d0) * (w(eta)/w(n)) * result(n,eta)
result(n,eta) = result(n,eta) * 1.0d0/sqrt(w(n))
if(n.ne.eta) result(eta,n) = result(eta,n) * 1.0d0/sqrt(w(eta))
end do
end do
return
end subroutine derLobattoMat


Expand All @@ -142,41 +145,41 @@ subroutine derLobatto(result,n,eta,Ntot,rpt,w)
! Calculates u'_n(r_eta)
!!!!!!
implicit none
real(rk), intent(out):: result
real(rk), intent(out):: result
integer(ik), intent(in):: n, Ntot,eta
real(rk), intent(in):: w(0:Ntot+1), rpt(0:Ntot+1)
real(rk):: pf, intresult,intpf, factor(0:Ntot+1),signrecorder
real(rk), intent(in):: w(0:Ntot+1), rpt(0:Ntot+1)
real(rk):: pf, intresult,intpf, factor(0:Ntot+1),signrecorder
integer(ik):: j,m,Nhalf

Nhalf=(Ntot+2)/2
if (n==eta) then
if (n==0) then
result=-1d0/(2d0*w(n))
else if (n==Ntot+1) then
result=1d0/(2d0*w(n))
else
result=0d0;
end if
else
result=1d0/((rpt(n)-rpt(eta)))
signrecorder = 1.0
do j=0,Ntot+1
if ((j .ne. n) .and. (j .ne. eta)) then
factor(j) = (rpt(eta)-rpt(j))/(rpt(n)-rpt(j))
signrecorder = signrecorder * (factor(j) / abs(factor(j)))
factor(j) = abs(factor(j))
else
factor(j) = 1.0d0
end if
end do
call QuickSort(factor,1,Ntot+2)
call Riffle(factor,Ntot+2,Nhalf)
do j = 0,Ntot+1
result = result * factor(j)
end do
result = result * signrecorder
end if
return
Nhalf=(Ntot+2)/2
if (n==eta) then
if (n==0) then
result=-1d0/(2d0*w(n))
else if (n==Ntot+1) then
result=1d0/(2d0*w(n))
else
result=0d0;
end if
else
result=1d0/((rpt(n)-rpt(eta)))
signrecorder = 1.0
do j=0,Ntot+1
if ((j .ne. n) .and. (j .ne. eta)) then
factor(j) = (rpt(eta)-rpt(j))/(rpt(n)-rpt(j))
signrecorder = signrecorder * (factor(j) / abs(factor(j)))
factor(j) = abs(factor(j))
else
factor(j) = 1.0d0
end if
end do
call QuickSort(factor,1,Ntot+2)
call Riffle(factor,Ntot+2,Nhalf)
do j = 0,Ntot+1
result = result * factor(j)
end do
result = result * signrecorder
end if
return
end subroutine derLobatto


Expand Down Expand Up @@ -232,57 +235,58 @@ end subroutine LobattoAbsWeights
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

recursive subroutine QuickSort(a, first, last)
implicit none
real(rk) a(*), x, t
integer(ik) first, last
integer(ik) i,j

x = a( (first + last) / 2)
i = first
j = last
do
do while (a(i) < x)
i = i+1
end do
do while (x < a(j))
j = j-1
end do
if (i >= j) exit
t = a(i); a(i) = a(j); a(j) = t
i = i+1
j = j-1
end do
if (first < i-1) call quicksort(a,first,i-1)
if (j+1 < last) call quicksort(a,j+1,last)
implicit none
real(rk) a(*), x, t
integer(ik) first, last
integer(ik) i,j

x = a( (first + last) / 2)
i = first
j = last
do
do while (a(i) < x)
i = i+1
end do
do while (x < a(j))
j = j-1
end do
if (i >= j) exit
t = a(i); a(i) = a(j); a(j) = t
i = i+1
j = j-1
end do
if (first < i-1) call quicksort(a,first,i-1)
if (j+1 < last) call quicksort(a,j+1,last)
end subroutine quicksort

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subroutine Riffle(a,NAbsc,Nhalf)
implicit none
real(rk) :: a(1:NAbsc),b(1:Nhalf),c(1:Nhalf)
integer(ik) NAbsc,i,Nhalf,Nhalfplusone

Nhalfplusone = Nhalf+1
if(mod(NAbsc,2).eq.0) then
do i=1,Nhalf
b(i) = a(i)
c(i) = a(NAbsc-i+1)
enddo
do i = 1,Nhalf
a(2*i) = c(i)
a(2*i-1) = b(i)
enddo
else
do i=1,Nhalfplusone
b(i) = a(i)
c(i) = a(NAbsc-i+1)
enddo
do i = 1,Nhalfplusone
if(i.ne.Nhalfplusone) a(2*i) = c(i)
a(2*i-1) = b(i)
enddo
endif
implicit none
integer(ik),intent(in) :: NAbsc,Nhalf
real(rk) :: a(1:NAbsc),b(1:Nhalf),c(1:Nhalf)
integer(ik) i,Nhalfplusone

Nhalfplusone = Nhalf+1
if(mod(NAbsc,2).eq.0) then
do i=1,Nhalf
b(i) = a(i)
c(i) = a(NAbsc-i+1)
enddo
do i = 1,Nhalf
a(2*i) = c(i)
a(2*i-1) = b(i)
enddo
else
do i=1,Nhalfplusone
b(i) = a(i)
c(i) = a(NAbsc-i+1)
enddo
do i = 1,Nhalfplusone
if(i.ne.Nhalfplusone) a(2*i) = c(i)
a(2*i-1) = b(i)
enddo
endif
end subroutine Riffle

end module Lobatto
11 changes: 7 additions & 4 deletions MAKEFILES/makefile_g95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ LIB = $(LAPACK)

###############################################################################

OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o symmetry.o dipole.o header_info.o
OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o symmetry.o dipole.o header_info.o Lobatto.o
#compilation_details.o

duo.x: $(OBJ) duo.o
Expand All @@ -27,10 +27,10 @@ duo.x: $(OBJ) duo.o
duo.o: duo.f90 $(OBJ)
$(FOR) -c duo.f90 $(FFLAGS)

grids.o: grids.f90 accuracy.o input.o
grids.o: grids.f90 accuracy.o input.o Lobatto.o
$(FOR) -c grids.f90 $(FFLAGS)

diatom.o: diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o
diatom.o: diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o Lobatto.o
$(FOR) -c diatom.f90 $(FFLAGS)

refinement.o: refinement.f90 accuracy.o input.o lapack.o diatom.o
Expand All @@ -39,7 +39,7 @@ refinement.o: refinement.f90 accuracy.o input.o lapack.o diatom.o
functions.o: functions.f90 accuracy.o input.o lapack.o
$(FOR) -c functions.f90 $(FFLAGS)

dipole.o: dipole.f90 accuracy.o input.o lapack.o diatom.o
dipole.o: dipole.f90 accuracy.o input.o lapack.o diatom.o
$(FOR) -c dipole.f90 $(FFLAGS)

accuracy.o: accuracy.f90
Expand All @@ -66,5 +66,8 @@ header_info.o: accuracy.o
atomic_and_nuclear_data.o: atomic_and_nuclear_data.f90
$(FOR) -c atomic_and_nuclear_data.f90 $(FFLAGS)

Lobatto.o: Lobatto.f90 timer.o
$(FOR) -c Lobatto.f90 $(FFLAGS)

clean:
rm -f $(OBJ) *.mod *__genmod.f90 duo.o
13 changes: 8 additions & 5 deletions MAKEFILES/makefile_gfortran
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ LIB = $(LAPACK)

###############################################################################

OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o symmetry.o dipole.o header_info.o
OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o symmetry.o dipole.o header_info.o Lobatto.o
#compilation_details.o

duo.x: $(OBJ) duo.o
Expand All @@ -29,10 +29,10 @@ duo.x: $(OBJ) duo.o
duo.o: duo.f90 $(OBJ)
$(FOR) -c duo.f90 $(FFLAGS)

grids.o: grids.f90 accuracy.o input.o
grids.o: grids.f90 accuracy.o input.o Lobatto.o
$(FOR) -c grids.f90 $(FFLAGS)

diatom.o: diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o
diatom.o: diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o Lobatto.o
$(FOR) -c diatom.f90 $(FFLAGS)

refinement.o: refinement.f90 accuracy.o input.o lapack.o diatom.o
Expand All @@ -41,7 +41,7 @@ refinement.o: refinement.f90 accuracy.o input.o lapack.o diatom.o
functions.o: functions.f90 accuracy.o input.o lapack.o
$(FOR) -c functions.f90 $(FFLAGS)

dipole.o: dipole.f90 accuracy.o input.o lapack.o diatom.o
dipole.o: dipole.f90 accuracy.o input.o lapack.o diatom.o
$(FOR) -c dipole.f90 $(FFLAGS)

accuracy.o: accuracy.f90
Expand All @@ -67,6 +67,9 @@ header_info.o: accuracy.o

atomic_and_nuclear_data.o: atomic_and_nuclear_data.f90
$(FOR) -c atomic_and_nuclear_data.f90 $(FFLAGS)


Lobatto.o: Lobatto.f90 timer.o
$(FOR) -c Lobatto.f90 $(FFLAGS)

clean:
rm -f $(OBJ) *.mod *__genmod.f90 duo.o
16 changes: 11 additions & 5 deletions MAKEFILES/makefile_nag
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,19 +19,19 @@ LIB = $(LAPACK)

###############################################################################

OBJ = grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o symmetry.o dipole.o header_info.o
OBJ = atomic_and_nuclear_data.o grids.o accuracy.o lapack.o timer.o input.o diatom.o refinement.o functions.o symmetry.o dipole.o header_info.o Lobatto.o
#compilation_details.o

duo.x: $(OBJ) duo.o
$(FOR) -o duo$(PLAT).x $(OBJ) $(FFLAGS) duo.o $(LIB)
$(FOR) -o j-duo$(PLAT).x $(OBJ) $(FFLAGS) duo.o $(LIB)

duo.o: duo.f90 $(OBJ)
$(FOR) -c duo.f90 $(FFLAGS)

grids.o: grids.f90 accuracy.o input.o
grids.o: grids.f90 accuracy.o input.o Lobatto.o
$(FOR) -c grids.f90 $(FFLAGS)

diatom.o: diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o
diatom.o: diatom.f90 accuracy.o input.o lapack.o functions.o symmetry.o atomic_and_nuclear_data.o Lobatto.o
$(FOR) -c diatom.f90 $(FFLAGS)

refinement.o: refinement.f90 accuracy.o input.o lapack.o diatom.o
Expand All @@ -40,7 +40,7 @@ refinement.o: refinement.f90 accuracy.o input.o lapack.o diatom.o
functions.o: functions.f90 accuracy.o input.o lapack.o
$(FOR) -c functions.f90 $(FFLAGS)

dipole.o: dipole.f90 accuracy.o input.o lapack.o diatom.o
dipole.o: dipole.f90 accuracy.o input.o lapack.o diatom.o
$(FOR) -c dipole.f90 $(FFLAGS)

accuracy.o: accuracy.f90
Expand All @@ -64,5 +64,11 @@ input.o: input.f90
header_info.o: accuracy.o
$(FOR) -c header_info.f90 $(FFLAGS)

atomic_and_nuclear_data.o: atomic_and_nuclear_data.f90
$(FOR) -c atomic_and_nuclear_data.f90 $(FFLAGS)

Lobatto.o: Lobatto.f90 timer.o
$(FOR) -c Lobatto.f90 $(FFLAGS)

clean:
rm -f $(OBJ) *.mod *__genmod.f90 duo.o
Loading

0 comments on commit 9e405eb

Please sign in to comment.