📄 init.f90
字号:
!--------------------
REAL, INTENT (OUT) :: x (4)
! --------------------------------------------------------------------
! | Returns a random 4-vector, 'x', with uniform distribution over |
! | the unit 3-sphere. |
! --------------------------------------------------------------------
INTEGER seed
REAL stored_normal, RAND
LOGICAL stored
COMMON / random / seed, stored_normal, stored
INTEGER i
REAL l, lsq
!------ Generate 'x' uniformly distributed inside a 4-cube (with sides [-1,1])
!------ Accept it if it's inside the unit 3-sphere, otherwise generate
!------ a new 'x' :
lsq = 2.0
DO WHILE ((lsq > 1.0) .OR. (lsq < 1.0e-6))
DO i = 1, 4
x(i) = (2.0 * RAND (seed)) - 1.0
ENDDO
lsq = DOT_PRODUCT (x, x)
ENDDO
!------ 'x' is now uniformly distributed INSIDE the unit 3-sphere (S3).
!------ Normalise it to length 1, so it is uniformly distributed ON S3 :
l = SQRT (lsq)
x = x / l
END
! }}}
! {{{ REAL FUNCTION normal ()
REAL FUNCTION normal ()
!----------------------
! ----------------------------------------------------------------------
! | Returns a random REAL value, normally distributed with standard |
! | deviation = SQRT (2). |
! | The algorithm actually generates 2 independent normals at once, |
! | so the second is saved and returned at the next call. |
! ----------------------------------------------------------------------
INTEGER seed
REAL stored_normal, RAND
LOGICAL stored
COMMON / random / seed, stored_normal, stored
REAL twopi
PARAMETER (twopi = 6.28318531)
REAL l, theta, x, y
IF (stored) THEN
normal = stored_normal
stored = .FALSE.
ELSE
!------ generate another 2 normals, and store one and return the other
x = 0.0
DO WHILE (x == 0.0)
x = RAND (seed)
y = RAND (seed)
ENDDO
l = SQRT (- LOG (x))
theta = twopi * y
normal = l * SIN (theta)
stored_normal = l * COS (theta)
stored = .TRUE.
ENDIF
END
! }}}
! {{{ SUBROUTINE randvec (u) -- SEQ
SUBROUTINE randvec (u)
!---------------------
INCLUDE 'procs.inc'
INCLUDE 'consts.inc'
REAL, INTENT (OUT), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: u
!HPF$ DISTRIBUTE * (*, BLOCK, BLOCK, BLOCK, *, *) ONTO * procs :: u
! ----------------------------------------------------------------------
! | Randomise the vector {u} : |
! | {u} = {z1 + i.z2} |
! | where z1 & z2 are normally distributed with std. dev. = SQRT(2) |
! ----------------------------------------------------------------------
INTEGER t, x, y, z, color, cmplx
REAL normal
DO cmplx = 1, ncmplx
DO color = 1, ncolor
DO z = 1, nx
DO y = 1, nx
DO x = 1, nx
DO t = 1, nt/2
u (t, x, y, z, color, cmplx) = normal ()
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
END
! }}}
! {{{ SUBROUTINE F_u (Fv, v) -- [TIMED] -- no comms
SUBROUTINE F_u (Fv, v)
!----------------------
INCLUDE 'procs.inc'
INCLUDE 'consts.inc'
REAL, INTENT (OUT), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: Fv
REAL, INTENT (IN), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: v
!HPF$ DISTRIBUTE * (*, BLOCK, BLOCK, BLOCK, *, *) ONTO * procs :: v, Fv
! ----------------------------------------------------------
! | Performs the [matrix] * {fermion.vector} product : |
! | |
! | {Fv} = [F] {v} |
! | = ( -[D][D] + m.sq ) {v} |
! | |
! | where {v} and {Fv} are defined on EVEN sites. |
! ----------------------------------------------------------
REAL, DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: Dv
!HPF$ DISTRIBUTE (*, BLOCK, BLOCK, BLOCK, *, *) ONTO procs :: Dv
CALL D_u (Dv, v, plus, even)
CALL D_u (Fv, Dv, minus, odd)
Fv = Fv + msq * v
END
! }}}
! {{{ SUBROUTINE D_u (Du, u, sign, usites) -- [TIMED] -- comms via CSHIFT (1)
SUBROUTINE D_u (Du, u, sign, usites)
!------------------------------------
INCLUDE 'procs.inc'
INCLUDE 'consts.inc'
REAL, INTENT (OUT), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: Du
REAL, INTENT (IN), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: u
!HPF$ DISTRIBUTE * (*, BLOCK, BLOCK, BLOCK, *, *) ONTO * procs :: Du, u
INTEGER, INTENT (IN) :: sign, usites
INCLUDE 'even.inc' ! for 'even_xyz'
INCLUDE 'gauge.inc' ! for array 'gauge'
INTEGER dir, vsites
REAL, DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: temp, v
!HPF$ DISTRIBUTE (*, BLOCK, BLOCK, BLOCK, *, *) ONTO procs :: temp, v
! {{{ < Description >
! -----------------------------------------------------------------------
! | Performs the [matrix] * {vector} product : |
! | |
! | {Du} = +/- [D] {u} |
! | |
! | [D] is the Dirac matrix |
! | {u}, {Du} are fermion vectors (defined only at even or odd sites)|
! | sign = 'plus' or 'minus' (the sign on the rhs of above eqn) |
! | usites = 'odd' or 'even' (the parity of the sites at which |
! | {u} is defined). |
! | |
! | Explanation |
! | ----------- |
! | [D] only has non-zero elements between sites at opposite ends |
! | of the same link. |
! | Consider a general site 'i'. Let '(i + mu)' represent the site|
! | at the end of the link emerging from site 'i' in the POSITIVE 'mu'|
! | direction, and '(i - mu)' represent the site at the end of the |
! | NEGATIVE 'mu' link. Then, {Du} at site 'i' is the sum of 8 terms |
! | (one for each link from 'i'): |
! | |
! | |
! | {Du}i = SUM D(i,i+mu) * u(i+mu) + SUM D(i,i-mu) * u(i-mu) |
! | +ve dirns -ve dirns |
! | mu mu |
! | |
! | |
! | Each term in the sum is evaluated (for all sites 'i') by a |
! | call to D_u_dir. D_u_dir contains a loop over sites, which may |
! | be vectorised. |
! |---------------------------------------------------------------------|
! }}}
! {{{ accumulate 'Du' for -ve directions
! {{{ < Explanation for -ve directions >
!---------------------------------------------------------------------|
! (I) Negative directions |
! ------------------- |
! |
! v (i) = D (i,j) * u (j), where site j = (i-mu). |
! |
! u(j) v(i) |
! <---- D(i,j) ----- |
! j i |
! |
! In this case : |
! D(i,j) = - herm. conj. [ gauge(j,i) ], |
! so the relevant gauge matrix is associated with the site 'j' at |
! the end of the link (where u (j) is defined). |
! To vectorise over sites, the elements of {u}, {v} and |
! [gauge] must be aligned (i.e. corresponding elements stored at |
! corresponding positions). Hence: |
! |
! (1) Call D_u_dir to evaluate : |
! |
! 'temp (j)' = - herm.conj. [gauge(j,i)] * u(j) |
! |
! for all 'usites' j (in D_u_dir). (N.B. By 'usites' we mean |
! the sites, even or odd, at which vector {u} is defined). |
! Now the result for site 'i', v (i), is held in the vector |
! {temp} at a position corresponding to site 'j'. |
! |
! (2) Shift all the 'temp (j)' in the positive direction |
! (along the link) to the corresponding sites 'i': |
! temp (j) --> v (i) (using CSHIFT). |
! |
! (3) Accumulate {v} in {Du}. |
! |
!----------------------------------------------------------------------
! }}}
!------ - t direction :
CALL D_u_dir (temp, gauge (:,:,:,:,:, tdir, usites), u, negsense)
!---------------------------------------------------------------
! Now do a +ve 't' shift of 'temp' into 'Du' (i.e. CSHIFT (shift = -1)!).
! Note that only half the 't' sites are stored. The actual 't'
! indices ('real t' below) are stored in index 'effective t' as
! shown in the following, which also shows the situation after the
! '+t' shift:
!
! effective t = | 1 | 2 | 3 | 4 | 5 |...
! real t {before +t shift = | 1 2 | 3 4 | 5 6 | 7 8 | 9 10 |...
! {after +t shift = | 1 | 2 3 | 4 5 | 6 7 | 8 9 |...
!
! Hence, only the *even* 'real t' sites have their 'effective t'
! indices shifted; the 'odd t' sites have their 'effective t' unchanged.
!
! Hence, shift temp (even t)
! -> shift temp (x,y,z) where:
! parity (x,y,z) == parity (temp) == parity (u)
!---------------------------------------------------------------
! parity of (x,y,z)
! even odd
IF (usites == even) THEN ! ---- ---
Du = CSHIFT (temp, - even_xyz, 1) ! shift = ( -1, 0)
ELSE
Du = CSHIFT (temp, even_xyz - 1, 1) ! shift = ( 0, -1)
ENDIF
!------ - x, y & z directions :
DO dir = xdir, zdir
CALL D_u_dir (temp, gauge (:,:,:,:,:, dir, usites), u, negsense)
Du = Du + CSHIFT (temp, -1, dir)
ENDDO
! }}}
! {{{ accumulate 'Du' for +ve directions
! {{{ < Explanation for +ve directions >
!----------------------------------------------------------------------
! (II) Positive directions |
! ------------------- |
! |
! v (i) = D (i,j) * u (j), where site j = (i+mu). |
! |
! v(i) u(j) |
! ----- D(i,j) ----> |
! i j |
! |
! In this case, D(i,j) = gauge(i,j), so the appropriate |
! gauge matrix is associated with site 'i' (at the start of the |
! link, the same as v (i)). |
! To vectorise over sites, the elements of {u}, {v} and |
! [gauge] must be aligned (i.e. corresponding elements stored at |
! corresponding positions). Hence: |
! |
! (1) Shift all the u (j) in the negative direction (along the link) |
! to the corresponding sites 'i': u (j) --> temp (i) |
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -