📄 cg2.f90
字号:
READ (10) gauge (:,:,:,:, igauge, idir, ipar)
ENDDO
ENDDO
ENDDO
CLOSE (10)
END
! }}}
! {{{ SUBROUTINE init_sites ()
SUBROUTINE init_sites ()
!----------------------
INCLUDE 'consts.inc'
INCLUDE 'procs.inc'
INCLUDE 'sites.inc' ! for arrays 'psi' & 'psidot'
! ---------------------------------------------------------------
! | Read the initial values of arrays 'psi' & 'psidot' from a |
! | file. |
! ---------------------------------------------------------------
OPEN (10, FILE="sites.dat", FORM="UNFORMATTED", STATUS="OLD")
DO icmplx = 1, ncmplx
DO icolor = 1, ncolor
READ (10) psi (:,:,:,:, icolor, icmplx)
ENDDO
ENDDO
DO icmplx = 1, ncmplx
DO icolor = 1, ncolor
READ (10) psidot (:,:,:,:, icolor, icmplx)
ENDDO
ENDDO
CLOSE (10)
END
! }}}
! {{{ SUBROUTINE F_u (Fv, v) -- [TIMED] -- no comms
SUBROUTINE F_u (Fv, v)
!----------------------
INCLUDE 'consts.inc'
INCLUDE 'procs.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$ ALIGN (*, :, :, :, *, *) WITH *tmpl :: Fv, v
! ----------------------------------------------------------
! | 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$ ALIGN (*, :, :, :, *, *) WITH tmpl :: 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 'consts.inc'
INCLUDE 'procs.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$ ALIGN (*, :, :, :, *, *) WITH *tmpl :: 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$ ALIGN (*, :, :, :, *, *) WITH tmpl :: 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 :
!---------------------------------------------------------------
! 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)
!---------------------------------------------------------------
CALL D_u_dir (Du, gauge (:,:,:,:,:, tdir, usites), u, negsense)
temp = CSHIFT (Du, -1, 1) ! parity of (x,y,z)
! even odd
IF (usites == even) THEN ! ---- ---
WHERE (even_xyz) Du = temp ! shift = ( -1, 0)
ELSE
WHERE (.NOT. even_xyz) Du = temp ! 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) |
! (using CSHIFT). |
! |
! (2) Perform |
! v(i) = gauge(i,j) * temp(i) |
! |
! for all 'vsites' i (in D_u_dir). (N.B. 'vsites' means the |
! sites, even or odd, at which vector {v} is defined). |
! The result v (i) for site 'i' is now held at the correct |
! position in vector {v}, so no further shift is necessary. |
! |
! (3) Accumulate {v} in {Du}. |
! |
!----------------------------------------------------------------------
! }}}
vsites = 1 - usites
!------ + t direction :
!---------------------------------------------------------------
! Now do a -ve 't' shift of 'u' into 'temp' (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 = | 2 3 | 4 5 | 6 7 | 8 9 |10 11 |...
!
! Hence, only the *odd* 'real t' sites have their 'effective t'
! indices shifted; the 'even t' sites have their 'effective t' unchanged.
!
! Hence, shift u (odd t)
! -> shift u (x,y,z) where:
! parity (x,y,z) != parity (u)
!---------------------------------------------------------------
temp = CSHIFT (u, 1, 1) ! parity of (x,y,z)
! even odd
IF (usites == even) THEN ! ---- ---
WHERE (even_xyz) temp = u ! shift = ( 0, 1)
ELSE
WHERE (.NOT. even_xyz) temp = u ! shift = ( 1, 0)
ENDIF
CALL D_u_dir (v, gauge (:,:,:,:,:, tdir, vsites), temp, possense)
Du = Du + v
!------ + x, y & z directions :
DO dir = xdir, zdir
temp = CSHIFT (u, 1, dir)
CALL D_u_dir (v, gauge (:,:,:,:,:, dir, vsites), temp, possense)
Du = Du + v
ENDDO
! }}}
!------Must still scale by +/- (1/2)
IF (sign == plus) THEN
Du = 0.5 * Du
ELSE
Du = - 0.5 * Du
ENDIF
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -