📄 init.f90
字号:
! (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)
!---------------------------------------------------------------
! parity of (x,y,z)
! even odd
IF (usites == even) THEN ! ---- ---
temp = CSHIFT (u, 1 - even_xyz, 1) ! shift = ( 0, 1)
ELSE
temp = CSHIFT (u, even_xyz, 1) ! 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
END
! }}}
! {{{ SUBROUTINE D_u_dir (v, gauge, u, sense) -- [TIMED] -- no comms
SUBROUTINE D_u_dir (v, gauge, u, sense)
!--------------------------------------
INCLUDE 'procs.inc'
INCLUDE 'consts.inc'
REAL, INTENT (OUT), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: v
REAL, INTENT (IN), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: u
!HPF$ DISTRIBUTE * (*, BLOCK, BLOCK, BLOCK, *, *) ONTO * procs :: v, u
REAL, INTENT (IN), DIMENSION (nt/2, nx, nx, nx, ngauge) :: gauge
!HPF$ DISTRIBUTE * (*, BLOCK, BLOCK, BLOCK, *) ONTO * procs :: gauge
INTEGER, INTENT (IN) :: sense
! -----------------------------------------------------------------------
! Calculates : |
! |
! v (i) = D (i,j) * u (j), |
! |
! for all sites 'i' at which {v} is defined (i.e. all even or odd |
! sites), for a particular link direction and sense. |
! |
! 'sense' = sense of the link (ij) |
! (positive = 'possense', negative = 'negsense' |
! |
! 'D' = [ g4 + i.g3 g1 + i.g2 ] |
! [-g1 + i.g2 g4 - i.g3 ] |
! |
! This subroutine is called with the arrays {u}, {v} and |
! [gauge] all aligned, such that the corresponding elements, |
! v (i), gauge (i,j) and u (j), are stored at the same array |
! positions. Thus the inner loop over sites can be VECTORISED ! |
! |
! -----------------------------------------------------------------------
! {{{ IF (sense == possense) THEN
IF (sense == possense) THEN
! ----------------------------
!----------------------------------------------------------------------
! Positive sense link |
! ------------------- |
! v (i) = D (i,j) * u (j), where site j = (i+mu). |
! |
! In this case, D(i,j) = gauge(i) (only that segment of |
! the total gauge field array corresponding to the given link |
! direction and sense is passed in). Also, the {u} vector |
! that is input has been shifted so that the sites are aligned: |
! u (i) = 'real' u (j). Thus, we calculate: |
! |
! v (i) = gauge (i) * u (i) for all sites 'i' |
!----------------------------------------------------------------------
v (:,:,:,:, 1, re) = gauge (:,:,:,:, 4) * u (:,:,:,:, 1, re) &
- gauge (:,:,:,:, 3) * u (:,:,:,:, 1, im) &
+ gauge (:,:,:,:, 1) * u (:,:,:,:, 2, re) &
- gauge (:,:,:,:, 2) * u (:,:,:,:, 2, im)
v (:,:,:,:, 1, im) = gauge (:,:,:,:, 4) * u (:,:,:,:, 1, im) &
+ gauge (:,:,:,:, 3) * u (:,:,:,:, 1, re) &
+ gauge (:,:,:,:, 1) * u (:,:,:,:, 2, im) &
+ gauge (:,:,:,:, 2) * u (:,:,:,:, 2, re)
v (:,:,:,:, 2, re) = - gauge (:,:,:,:, 1) * u (:,:,:,:, 1, re) &
- gauge (:,:,:,:, 2) * u (:,:,:,:, 1, im) &
+ gauge (:,:,:,:, 4) * u (:,:,:,:, 2, re) &
+ gauge (:,:,:,:, 3) * u (:,:,:,:, 2, im)
v (:,:,:,:, 2, im) = - gauge (:,:,:,:, 1) * u (:,:,:,:, 1, im) &
+ gauge (:,:,:,:, 2) * u (:,:,:,:, 1, re) &
+ gauge (:,:,:,:, 4) * u (:,:,:,:, 2, im) &
- gauge (:,:,:,:, 3) * u (:,:,:,:, 2, re)
! }}}
! {{{ ELSE
ELSE
! ----
!----------------------------------------------------------------------
! Negative sense link |
! ------------------- |
! v (i) = D (i,j) * u (j), where site j = (i-mu). |
! |
! In this case, D(i,j) = - herm.conj [ gauge(j,i) ] |
! = - herm.conj [ gauge (j) ] |
! |
! since only that segment of the total gauge field array |
! corresponding to the given link direction and sense is passed |
! in. Also, the {v} vector that is input has been shifted so |
! that the sites are aligned: v (j) = 'real' v (i). |
! Thus, we calculate: |
! |
! v (j) = - herm.conj [ gauge (j) ] * u (j) for all sites 'j' |
!----------------------------------------------------------------------
v (:,:,:,:, 1, re) = - gauge (:,:,:,:, 4) * u (:,:,:,:, 1, re) &
- gauge (:,:,:,:, 3) * u (:,:,:,:, 1, im) &
+ gauge (:,:,:,:, 1) * u (:,:,:,:, 2, re) &
- gauge (:,:,:,:, 2) * u (:,:,:,:, 2, im)
v (:,:,:,:, 1, im) = - gauge (:,:,:,:, 4) * u (:,:,:,:, 1, im) &
+ gauge (:,:,:,:, 3) * u (:,:,:,:, 1, re) &
+ gauge (:,:,:,:, 1) * u (:,:,:,:, 2, im) &
+ gauge (:,:,:,:, 2) * u (:,:,:,:, 2, re)
v (:,:,:,:, 2, re) = - gauge (:,:,:,:, 1) * u (:,:,:,:, 1, re) &
- gauge (:,:,:,:, 2) * u (:,:,:,:, 1, im) &
- gauge (:,:,:,:, 4) * u (:,:,:,:, 2, re) &
+ gauge (:,:,:,:, 3) * u (:,:,:,:, 2, im)
v (:,:,:,:, 2, im) = - gauge (:,:,:,:, 1) * u (:,:,:,:, 1, im) &
+ gauge (:,:,:,:, 2) * u (:,:,:,:, 1, re) &
- gauge (:,:,:,:, 4) * u (:,:,:,:, 2, im) &
- gauge (:,:,:,:, 3) * u (:,:,:,:, 2, re)
ENDIF
! }}}
END
! }}}
! {{{ SUBROUTINE solve (x, b) -- [TIMED] -- no comms in timed part
SUBROUTINE solve (x, b)
!----------------------
INCLUDE 'procs.inc'
INCLUDE 'consts.inc'
REAL, INTENT (INOUT), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: x
REAL, INTENT (IN), DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: b
!HPF$ DISTRIBUTE * (*, BLOCK, BLOCK, BLOCK, *, *) ONTO * procs :: x, b
! -------------------------------------------------------------------
! | Uses the CG algorithm to solve : |
! | [F] {x} = {b} |
! | for {x}, given {b}. |
! | |
! | {} denotes a fermion vector (defined on even sites) |
! | [F] = (- [D][D] + msq), where [D] is the Dirac matrix. |
! | |
! | At start, {x} has the initial trial solution. |
! -------------------------------------------------------------------
INTEGER niter
REAL rsq, rsqtol
REAL, DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: p, r
!HPF$ DISTRIBUTE (*, BLOCK, BLOCK, BLOCK, *, *) ONTO procs :: p, r
! ---------------------------------------------------
! | Initialise {p}, {r}, rsq, rsqtol and niter |
! | |
! | {p} = {r} = {b} - {Fx} |
! | rsq = {r}.{r} |
! ---------------------------------------------------
CALL F_u (r, x)
r = b - r
p = r
rsqtol = SUM (b * b) * (1.0E-02 ** (sigfigs))
rsq = SUM (r * r) ! hermitian scalar product, (r+, r)
niter = 0
! ------------------------------------------
! | Do CG iterations until {x} converged |
! ------------------------------------------
1 CONTINUE
niter = niter + 1
CALL CGiter (x, p, r, rsq)
IF ((rsq > rsqtol) .AND. (niter < maxiter)) go to 1
END
! }}}
! {{{ SUBROUTINE CGiter (x, p, r, rsq) [TIMED] -- comms via SUM
SUBROUTINE CGiter (x, p, r, rsq)
!-------------------------------
INCLUDE 'procs.inc'
INCLUDE 'consts.inc'
REAL, INTENT (INOUT), &
DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: x, p, r
!HPF$ DISTRIBUTE * (*, BLOCK, BLOCK, BLOCK, *, *) ONTO * procs :: x, p, r
REAL, INTENT (INOUT) :: rsq
REAL, DIMENSION (nt/2, nx, nx, nx, ncolor, ncmplx) :: Fp
!HPF$ DISTRIBUTE (*, BLOCK, BLOCK, BLOCK, *, *) ONTO procs :: Fp
REAL alpha, newrsq
! ------------------------------------------------------------
! | One iteration of the Conjugate Gradient (CG) algorithm |
! | - no preconditioning. |
! ------------------------------------------------------------
! -----------------------------------
! | alpha = (r, r) / (p, Fp) |
! -----------------------------------
CALL F_u (Fp, p)
alpha = rsq / SUM (p * Fp) ! hermitian scalar product, (p+, Fp)
! --------------------------------------
! | {new.x} = {x} + alpha * {p} |
! --------------------------------------
x = x + alpha * p
! --------------------------------------
! | {new.r} = {r} - alpha * {Fp} |
! --------------------------------------
r = r - alpha * Fp
! ----------------------------------------------------------------
! | {new.p} = ((new.r, new.r) / (r, r)) * {p} + {new.r} |
! ----------------------------------------------------------------
newrsq = SUM (r * r) ! hermitian scalar product, (r+, r)
p = r + (newrsq / rsq) * p
rsq = newrsq
END
! }}}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -