Hi All,
Thank you for reply.
> No attachment.
>Try attachment as .txt.
I'm sorry.
In the previous mail, I attached a tar+gziped file sample.tgz.
It seems, binary files are not accepted.
Now, I try files with .txt.
I renamed Makefile to Makefile.txt and sample.f to sample.txt.
And I attach these two text files. I hope I can send this time.
Best,
Takashi
On 08/06/20 06:53, NightStrike wrote:
No attachment
On Wed, Aug 5, 2020, 12:01 PM Takashi Inoue <[email protected]>
wrote:
Hi All,
This is my first post to here. Nice to meet you.
I may found a bug in mingw-w64.
In particular, it is in OpenMP schedule(dynamic, chunk-size).
Let me ask your opinion.
One of my fortran program with "!$omp do schedule(dynamic,1)"
produce slightly different result at time to time.
This does not happen with GCC on FreeBSD and NAG on Windows.
So, I think nothing is wrong in my code. (But, I'm not 100% sure.)
Mingw-w64 which I installed is the latest package for msys2.
For those who are interested in this problem,
I attach a sample code which reproduce the problem.
(This may not be appropriate for mailing list. I'm sorry).
If you run the program several times, you will get slightly
different fort.99 at time to time.
Perhaps this problem has root in POSIX threading in mingw-w64.
Any comments and suggestions will be appreciated.
Best,
Takashi
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
.SUFFIXES: .f .o
FORTRAN = gfortran
FFLAGS = -fopenmp -O3 -Wall
PROG0 = sample.exe
OBJS0 = sample.o
$(PROG0) : $(OBJS0)
$(FORTRAN) $(FFLAGS) $(OBJS0) $(LIBS) -o $(PROG0)
.f.o:
$(FORTRAN) $(FFLAGS) -c $<
clean:
rm -f *.o *~ *.exe
c--------1---------2---------3---------4---------5---------6--------7--
c Example of possible bug in "omp do schedule(dynamic,1)" of mingw-w64
c--------1---------2---------3---------4---------5---------6--------7--
program Exsample
implicit none
integer NofSTO(1), Border(1)
integer LM(2,5), NoPgf(5)
integer Nolmn(5), lmn(3,3,5)
double precision Zeta(3,1), Posi(3,1)
double precision lmncf(3,5), ai(6,5), di(6,5)
double precision Norm(5,5)
double precision bohr
integer i,j
c--------1---------2---------3---------4---------5---------6--------7--
bohr = 0.52917721042381560d0
Zeta(1,1) = 7.66d0/bohr
Zeta(2,1) = 2.25d0/bohr
Zeta(3,1) = 0d0
NofSTO(1) = 5
Border(1) = 5
Posi(1,1) = 0d0
Posi(2,1) = 0d0
Posi(3,1) = 0d0
c--------1---------2---------3---------4---------5---------6--------7--
call makeLMlmn_STO3G(LM,NoPgf,Nolmn,lmn,lmncf)
call makeKLMG_STO3G(Zeta(:,1),ai,di)
call calcnorm(1,Posi,5,Border,
& LM,NoPgf,Nolmn,lmn,lmncf,ai,di,Norm)
do j=1,5
do i=1,5
write(99,*) i,j, Norm(i,j)
end do
end do
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c Return matrix
c--------1---------2---------3---------4---------5---------6--------7--
subroutine calcnorm(Natom,Posi,NofBase,Border,
& LM,NoPgf,Nolmn,lmn,lmncf,ai,di,Norm)
implicit none
integer,intent(in):: Natom,NofBase,Border(Natom)
integer,intent(in):: LM(2,NofBase),NoPgf(NofBase)
integer,intent(in):: Nolmn(NofBase),lmn(3,3,NofBase)
double precision, intent(in):: Posi(3,Natom),lmncf(3,NofBase)
double precision, intent(in):: ai(6,NofBase),di(6,NofBase)
double precision, intent(out):: Norm(NofBase,NofBase)
double precision, parameter :: Pi = 3.1415926535897932d0
double precision XA(3),XB(3),XO(3)
double precision alP,alQ,ciP,ciQ,NrP,NrQ
double precision OvI,temp
integer lmnP(3),lmnQ(3),Lp,Lq,Mp,Mq
integer pp,qq,a,b,p,q
integer i,j,ii,jj
c--------1---------2---------3---------4---------5---------6--------7--
call zero1dim(3,XO)
!$omp parallel num_threads(4)
!$omp+private(pp,qq,i,j,ii,jj)
!$omp+private(Lp,Mp,Lq,Mq,a,b,p,q)
!$omp+private(XA,XB)
!$omp+private(alP,NrP,ciP,lmnP)
!$omp+private(alQ,NrQ,ciQ,lmnQ)
!$omp+private(temp,OvI)
!$omp do schedule(dynamic,1)
do 20 qq=1, NofBase
Lq = LM(1,qq)
Mq = LM(2,qq)
call pp2ap(Natom,Border,qq,b,q)
do i=1,3
XB(i)=Posi(i,b)
end do
do 10 pp=qq, NofBase
Lp = LM(1,pp)
Mp = LM(2,pp)
call pp2ap(Natom,Border,pp,a,p)
do i=1,3
XA(i)=Posi(i,a)
end do
c----------------------------------------------------------------
temp = 0d0
do i=1,NoPgf(pp)
alP = ai(i,pp)
NrP = 1d0/4d0**Lp *Sqrt((Pi**3)/(8*alP**(3+2*Lp)))
ciP = di(i,pp) /Sqrt(NrP)
do ii=1,Nolmn(pp)
call copylmn(lmn(:,ii,pp),lmnP)
do j=1,NoPgf(qq)
alQ = ai(j,qq)
NrQ = 1d0/4d0**Lq *Sqrt((Pi**3)/(8*alQ**(3+2*Lq)))
ciQ = di(j,qq) /Sqrt(NrQ)
do jj=1,Nolmn(qq)
call copylmn(lmn(:,jj,qq),lmnQ)
if (a.eq.b) then
call OverlapInt(XO,alP,lmnP,XO,alQ,lmnQ,OvI)
else
call OverlapInt(XA,alP,lmnP,XB,alQ,lmnQ,OvI)
end if
temp = temp + ciP*ciQ*lmncf(ii,pp)*lmncf(jj,qq)*OvI
end do
end do
end do
end do
if (Lp.eq.Lq .and. Mp.eq.Mq) then
Norm(pp,qq) = temp
Norm(qq,pp) = temp
else
Norm(pp,qq) = 0d0
Norm(qq,pp) = 0d0
end if
c----------------------------------------------------------------
10 continue
20 continue
!$omp end do
!$omp end parallel
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c Overlap Integral
c--------1---------2---------3---------4---------5---------6--------7--
subroutine OverlapInt(XA,alA,lmnA,XB,alB,lmnB,OvI)
implicit none
double precision, intent(in):: XA(3),alA
double precision, intent(in):: XB(3),alB
integer,intent(in):: lmnA(3),lmnB(3)
double precision,intent(out):: OvI
double precision, parameter :: huge = 5d2
double precision, parameter :: Pi = 3.1415926535897932d0
double precision gamma,XP(3),PA(3),PB(3),AB(3),ABsq,Shld
double precision tempx,tempy,tempz
double precision fj
integer dbFactrial
integer i
if (lmnA(1).lt.0 .or. lmnA(2).lt.0 .or. lmnA(3).lt.0 .or.
& lmnB(1).lt.0 .or. lmnB(2).lt.0 .or. lmnB(3).lt.0 )then
OvI = 0d0
return
end if
gamma = alA + alB
do i=1,3
if (XA(i).eq.XB(i)) then
XP(i) = XA(i)
PA(i) = 0d0
PB(i) = 0d0
AB(i) = 0d0
else
XP(i) = (alA*XA(i)+alB*XB(i))/gamma
PA(i) = XP(i) - XA(i)
PB(i) = XP(i) - XB(i)
AB(i) = XA(i) - XB(i)
end if
end do
ABsq = AB(1)**2 + AB(2)**2 + AB(3)**2
Shld = alA*alB*ABsq/gamma
if (Shld.gt.huge) then
OvI = 0d0
return
end if
tempx = 0d0
do i=0,(lmnA(1)+lmnB(1))/2
tempx = tempx + fj(2*i,lmnA(1),lmnB(1),PA(1),PB(1))
& *dbFactrial(2*i-1)/(2*gamma)**i
end do
tempy = 0d0
do i=0,(lmnA(2)+lmnB(2))/2
tempy = tempy + fj(2*i,lmnA(2),lmnB(2),PA(2),PB(2))
& *dbFactrial(2*i-1)/(2*gamma)**i
end do
tempz = 0d0
do i=0,(lmnA(3)+lmnB(3))/2
tempz = tempz + fj(2*i,lmnA(3),lmnB(3),PA(3),PB(3))
& *dbFactrial(2*i-1)/(2*gamma)**i
end do
OvI = (Pi/gamma)**(3d0/2d0) *tempx *tempy *tempz *Exp(-Shld)
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c f_j(l,m,a,b)
c--------1---------2---------3---------4---------5---------6--------7--
function fj(j,l,m,a,b)
implicit none
integer,intent(in):: j,l,m
double precision, intent(in):: a,b
double precision fj
integer Combini
double precision powerfunc
integer i
if (l.lt.0 .or. m.lt.0) then
write(*,*) "l, m should be non-negative"
stop
end if
if (j.lt.0 .or. j.gt.l+m) then
write(*,*) "j should be in [0:l+m]"
stop
end if
fj = 0d0
do i=Max(0,j-m),Min(j,l)
fj = fj + Combini(l,i) * powerfunc(a,l-i)
& * Combini(m,j-i)* powerfunc(b,m+i-j)
end do
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c power function x**n
c--------1---------2---------3---------4---------5---------6--------7--
function powerfunc(x,n)
implicit none
integer,intent(in):: n
double precision, intent(in):: x
double precision powerfunc
if (n.eq.0) then
powerfunc = 1d0
else
powerfunc = x**n
end if
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c The combination factor
c--------1---------2---------3---------4---------5---------6--------7--
function Combini(n,m)
implicit none
integer,intent(in):: n,m
integer Combini
integer Factrial
Combini = Factrial(n)/(Factrial(m)*Factrial(n-m))
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c The factrial n!
c--------1---------2---------3---------4---------5---------6--------7--
function Factrial(n)
implicit none
integer,intent(in):: n
integer Factrial
if(n.eq.0.or.n.eq.1)then
Factrial = 1
else if (n.eq.2) then
Factrial = 2
else if (n.eq.3) then
Factrial = 6
else if (n.eq.4) then
Factrial = 24
else if (n.eq.5) then
Factrial = 120
else if (n.eq.6) then
Factrial = 720
else if (n.eq.7) then
Factrial = 5040
else if (n.eq.8) then
Factrial = 40320
else if (n.eq.9) then
Factrial = 362880
else if (n.eq.10) then
Factrial = 3628800
else if (n.eq.11) then
Factrial = 39916800
else if (n.eq.12) then
Factrial = 479001600
else if (n.gt.12) then
write(*,*)"Too large input in Factrial. n=", n
stop
else
write(*,*)"Negative input in Factrial. n=", n
stop
end if
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c The double-factrial n!!
c--------1---------2---------3---------4---------5---------6--------7--
function dbFactrial(n)
implicit none
integer,intent(in):: n
integer dbFactrial
integer i
if(n.eq.-1 .or. n.eq.0 .or. n.eq.1)then
dbFactrial = 1
else if (n.eq.2) then
dbFactrial = 2
else if (n.eq.3) then
dbFactrial = 3
else if (n.eq.4) then
dbFactrial = 8
else if (n.eq.5) then
dbFactrial = 15
else if (n.eq.6) then
dbFactrial = 48
else if (n.eq.7) then
dbFactrial = 105
else if (n.eq.8) then
dbFactrial = 384
else if (n.eq.9) then
dbFactrial = 945
else if (n.eq.10) then
dbFactrial = 3840
else if (n.eq.11) then
dbFactrial = 10395
else if (n.eq.12) then
dbFactrial = 46080
else if (n.eq.13) then
dbFactrial = 135135
else if (n.eq.14) then
dbFactrial = 645120
else if (n.gt.14) then
dbFactrial = n
do i = 1, n/2
if(n-2*i.ne.0)then
dbFactrial = dbFactrial * (n-2*i)
end if
end do
else
write(*,*)"Negative input in dbFactrial. n=", n
stop
end if
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c Convert pp to a,p
c--------1---------2---------3---------4---------5---------6--------7--
subroutine pp2ap(Natom,Border,pp,a,p)
implicit none
integer,intent(in):: Natom,Border(Natom),pp
integer,intent(out):: a,p
integer i
a=1
do i=2,Natom
if (pp.gt.Border(i)) then
a=a+1
end if
end do
p = pp - Border(a)
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c Copy lmn A to B
c--------1---------2---------3---------4---------5---------6--------7--
subroutine copylmn(A,B)
implicit none
integer,intent(in):: A(3)
integer,intent(out)::B(3)
integer i
do i=1,3
B(i) = A(i)
end do
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c Zero-clear dimension
c--------1---------2---------3---------4---------5---------6--------7--
subroutine zero1dim(n,A)
implicit none
integer,intent(in):: n
double precision, intent(out):: A(n)
integer i
do i=1,n
A(i) = 0d0
end do
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c Make NoPgf, {LM} and {lmn} for Solid-Spherical-Harmonics
c--------1---------2---------3---------4---------5---------6--------7--
subroutine makeLMlmn_STO3G(LM,NoPgf,Nolmn,lmn,lmncf)
implicit none
integer,intent(out):: LM(2,5),NoPgf(5),Nolmn(5),lmn(3,3,5)
double precision,intent(out):: lmncf(3,5)
integer i
c-----------
do i=1,5
NoPgf(i)= 3
Nolmn(i)= 1
lmncf(1,i)= 1d0
end do
c-----------
LM(1,1)= 0
LM(2,1)= 0
lmn(1,1,1)= 0
lmn(2,1,1)= 0
lmn(3,1,1)= 0
c-----------
LM(1,2)= 0
LM(2,2)= 0
lmn(1,1,2)= 0
lmn(2,1,2)= 0
lmn(3,1,2)= 0
c-----------
LM(1,3)= 1
LM(2,3)=-1
lmn(1,1,3)= 0
lmn(2,1,3)= 1
lmn(3,1,3)= 0
c-----------
LM(1,4)= 1
LM(2,4)= 0
lmn(1,1,4)= 0
lmn(2,1,4)= 0
lmn(3,1,4)= 1
c-----------
LM(1,5)= 1
LM(2,5)= 1
lmn(1,1,5)= 1
lmn(2,1,5)= 0
lmn(3,1,5)= 0
c-----------
do i=1,5
lmn(1,2,i)= 0
lmn(2,2,i)= 0
lmn(3,2,i)= 0
lmn(1,3,i)= 0
lmn(2,3,i)= 0
lmn(3,3,i)= 0
lmncf(2,i)= 0d0
lmncf(3,i)= 0d0
end do
c-----------
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c Make the STO3G basis set
c--------1---------2---------3---------4---------5---------6--------7--
subroutine makeKLMG_STO3G(Zeta,ai,di)
implicit none
double precision,intent(in):: Zeta(3)
double precision,intent(out):: ai(6,5),di(6,5)
integer i,j
c--1S----
ai(1,1) = 2.22766d0
ai(2,1) = 4.05771d-1
ai(3,1) = 1.09818d-1
di(1,1) = 1.54329d-1
di(2,1) = 5.35328d-1
di(3,1) = 4.44635d-1
c--2S----
ai(1,2) = 9.94203d-1
ai(2,2) = 2.31031d-1
ai(3,2) = 7.51386d-2
di(1,2) =-9.99672d-2
di(2,2) = 3.99513d-1
di(3,2) = 7.00115d-1
c--2Py---
ai(1,3) = ai(1,2)
ai(2,3) = ai(2,2)
ai(3,3) = ai(3,2)
di(1,3) = 1.55916d-1
di(2,3) = 6.07684d-1
di(3,3) = 3.91957d-1
c--2Pz---
ai(1,4) = ai(1,3)
ai(2,4) = ai(2,3)
ai(3,4) = ai(3,3)
di(1,4) = di(1,3)
di(2,4) = di(2,3)
di(3,4) = di(3,3)
c--2Px---
ai(1,5) = ai(1,3)
ai(2,5) = ai(2,3)
ai(3,5) = ai(3,3)
di(1,5) = di(1,3)
di(2,5) = di(2,3)
di(3,5) = di(3,3)
do i=1,3
ai(i,1) = ai(i,1)*Zeta(1)**2
ai(i,2) = ai(i,2)*Zeta(2)**2
ai(i,3) = ai(i,3)*Zeta(2)**2
ai(i,4) = ai(i,4)*Zeta(2)**2
ai(i,5) = ai(i,5)*Zeta(2)**2
end do
do j=1,5
do i=4,6
ai(i,j) = 0d0
di(i,j) = 0d0
end do
end do
return
end
c--------1---------2---------3---------4---------5---------6--------7--
c--------1---------2---------3---------4---------5---------6--------7--
c END
c--------1---------2---------3---------4---------5---------6--------7--
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
