mg.h90
来自「shpf 1.9一个并行编译器」· H90 代码 · 共 957 行 · 第 1/3 页
H90
957 行
!
! with no declared relation between NPC (= n+1) and NPF (= 2*n+1).
! Also, there was no alignment between the arrays FC and FF,
! so all the assignments beolow involved 'tentative' remaps.
! I have used 'n' throughout, which is related to old params thus:
! npc = n+1
! npf = 2*n + 1
! nc = n
!
! -- It would be nice if the following were true (since then no comms):
! ALIGN uc (:, :) WITH uf (1 : 2*n+1 : 2, 1 : 2*n+1 : 2)
!
! -- Bug in last stmt of original:
! UC (1 : NPC, NPC) = UF (1 : NPF : 2, NPC)
! ^^^
! should be NPF
! -- Originally had interior & 4 edges in separate assignments
! for some reason!!:
! UC (2 : n, 2 : n) = UF (3 : 2*n-1 : 2, 3 : 2*n-1 : 2)
! UC ( 1 , 2 : n) = UF ( 1 , 3 : 2*n-1 : 2)
! UC ( n+1 , 2 : n) = UF ( 2*n+1 , 3 : 2*n-1 : 2)
! UC (1 : n+1, 1 ) = UF (1 : 2*n+1 : 2, 1 )
! UC (1 : n+1, n+1 ) = UF (1 : 2*n+1 : 2, 2*n+1 )
! I've merged them into 1 assignment!
!
! -- Q: Why is 'fc' not set at the edges??
!------------------------------------------------------------------
! }}}
!-----------------------------------------------------------------
! Called by 'init ()'. Transfers 'f' & 'u' from a fine grid of
! size (2*n+1) * (2*n+1) to a coarse grid of size (n+1) * (n+1).
!-----------------------------------------------------------------
IMPLICIT NONE
INTEGER n
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: fc, uc
DOUBLE PRECISION, DIMENSION (2*n+1, 2*n+1) :: ff, uf
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: uc, fc
!HPF$ ALIGN WITH * tmplf :: uf, ff
uc = uf (1:2*n+1:2, 1:2*n+1:2) ! could have uf (: :2, : :2),
! but 'xlf90' doesn't like it!
fc (2:n, 2:n) = 4 * ff (3:2*n-1:2, 3:2*n-1:2)
!-----------------------------------------------------------------
! N.B. For the assignment to 'uc', The exact transcription of the
! original is:
!
! uc (:, 1:n) = uf (1 : 2*n+1 : 2, 1 : 2*n-1 : 2)
! uc (:, n+1) = uf (1 : 2*n+1 : 2, n+1)
! ^^^^^
! bug ??
!
! but presumably the 2nd assignment is a bug (used NPC rather than NPF)!!!
!-----------------------------------------------------------------
END
! }}}
! {{{ SUBROUTINE PUTZB (n, u) [start of level]
SUBROUTINE PUTZB (n, u)
!-----------------------------------------------------------------
! Called from 'init () & 'fmgdn ()'. Zeroes boundary of coarser gid.
!-----------------------------------------------------------------
! Mods
! ----
! -- Replaced parameter NPC (= n+1) by 'n', and made corresponding
! changes to code and calls.
!-----------------------------------------------------------------
IMPLICIT NONE
INTEGER n
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: u
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: u
u ( 1 , 2:n) = 0.0D0
u (n+1, 2:n) = 0.0D0
u ( : , 1 ) = 0.0D0
u ( : , n+1) = 0.0D0
END
! }}}
! {{{ SUBROUTINE INT2A (n, uc, uf) [TIMED]
SUBROUTINE INT2A (n, uc, uf)
!------------------------------------------------------------------
! Called from 'mgidn ()'.
! COARSE-TO-FINE TRANSFER (BILINEAR INTERPOLATION) AND CORRECTION.
! (THIS MODIFIES ODD SITES ONLY OF uf)
!------------------------------------------------------------------
IMPLICIT NONE
INTEGER n
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: uc
DOUBLE PRECISION, DIMENSION (2*n+1, 2*n+1) :: uf
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: uc
!HPF$ ALIGN WITH * tmplf :: uf
uf (2 : 2*n : 2, 3 : 2*n-1 : 2) = &
uf (2 : 2*n : 2, 3 : 2*n-1 : 2) &
+ 0.5D0 * (uc (1 : n, 2 : n) &
+ uc (2 : n+1, 2 : n))
uf (3 : 2*n-1 : 2, 2 : 2*n : 2) = &
uf (3 : 2*n-1 : 2, 2 : 2*n : 2) &
+ 0.5D0 * (uc (2 : n, 1 : n) &
+ uc (2 : n, 2 : n+1))
END
! }}}
! {{{ SUBROUTINE INT4 (n, uc, uf, ff) [start of level]
SUBROUTINE INT4 (n, uc, uf, ff)
!------------------------------------------------------------------
! Called from 'fmgdn ()'. 4-TH ORDER FMG-INTERPOLATION
!------------------------------------------------------------------
IMPLICIT NONE
INTEGER :: n
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: uc
DOUBLE PRECISION, DIMENSION (2*n+1, 2*n+1) :: uf, ff
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: uc
!HPF$ ALIGN WITH * tmplf :: ff, uf
!---------------------------------------------------------------!
! Transfer coarse grid values to the *odd* points of fine grid !
! (interior points only) !
!---------------------------------------------------------------!
uf (3:2*n-1:2, 3:2*n-1:2) = uc (2:n, 2:n)
!----------------------------------------------------------!
! Compute the *even* points of the fine grid using the !
! rotated 5-point star. !
!----------------------------------------------------------!
uf (2:2*n:2, 2:2*n:2) = 0.25d0 * (2 * ff (2:2*n:2, 2:2*n:2) &
+ uc (1:n, 1:n) &
+ uc (1:n, 2:n+1) &
+ uc (2:n+1, 1:n) &
+ uc (2:n+1, 2:n+1))
! COMPUTATION OF THE ODD POINTS BY ONE HALF (ODD) RELAXATION STEP
! **** WHAT?? -VE ITMAX DOESN'T DO THAT !! ****
! *** BUG -- FOLLOWING CALL RETURNS IMMEDIATELY! ***
CALL RELAX (-1, .FALSE., (2*n), uf, ff)
END
!***************************************************************
!*** ORIGINALLY HAD FOLLOWING (WITH BUGS)...
!
! hn (:, 2 : n+1) = uc (:, 1 : n) + uc (:, 2 : n+1)
!
!***** FOLLOWING ASSIGNMENT TO TEMPORARY IS COMPLETELY UNECESSARY...
! ha (2:n+1, 2:n+1) = hn (1:n, 2:n+1)
!
!****** BUG IN FOLLOWING ASSIGNMENT -- DOESN'T CONFORM!!
!****** PRESUMABLY MISSED-OUT STRIDE OF 2 IN FINE GRIDS!
!
! uf (2:2*n, 2:2*n) = 0.25d0 * (2 * ff (2:2*n, 2:2*n) &
! + ha (2:n+1, 2:n+1) &
! + hn (2:n+1, 2:n+1))
!***************************************************************
! }}}
! {{{ SUBROUTINE RELAX (ITMAX, IZ, n, u, f) [TIMED]
SUBROUTINE RELAX (ITMAX, IZ, n, u, f)
!--------------------------------------------------------------
! Called from 'mgi ()', 'mgiup ()', 'mgidn ()' and 'int4 ()',
! to perform red-black relaxation on gridsize (n+1).
!--------------------------------------------------------------
IMPLICIT NONE
INTEGER ITMAX, n ! (itmax == -1) when called from 'int4()' !!!
LOGICAL IZ
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: f, u
INTEGER it
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: u, f
DO IT = 1, ITMAX
!-----------------------------!
! Relaxation of even points !
!-----------------------------!
IF (IZ .AND. (IT <= 1)) THEN ! zero starting values
u (2: n :2, 2: n :2) = 0.25 * f (2: n :2, 2: n :2)
u (3:n-1:2, 3:n-1:2) = 0.25 * f (3:n-1:2, 3:n-1:2)
ELSE ! non-zero starting values
u (2: n :2, 2: n :2) = 0.25 * (f (2: n :2, 2: n :2) &
+ u (1:n-1:2, 2: n :2) &
+ u (3:n+1:2, 2: n :2) &
+ u (2: n :2, 1:n-1:2) &
+ u (2: n :2, 3:n+1:2))
u (3:n-1:2, 3:n-1:2) = 0.25 * (f (3:n-1:2, 3:n-1:2) &
+ u (2:n-2:2, 3:n-1:2) &
+ u (4: n :2, 3:n-1:2) &
+ u (3:n-1:2, 2:n-2:2) &
+ u (3:n-1:2, 4: n :2))
ENDIF
!-----------------------------!
! Relaxation of odd points !
!-----------------------------!
u (3:n-1:2, 2: n :2) = 0.25 * (f (3:n-1:2, 2: n :2) &
+ u (2:n-2:2, 2: n :2) &
+ u (4: n :2, 2: n :2) &
+ u (3:n-1:2, 1:n-1:2) &
+ u (3:n-1:2, 3:n+1:2))
u (2: n :2, 3:n-1:2) = 0.25 * (f (2: n :2, 3:n-1:2) &
+ u (1:n-1:2, 3:n-1:2) &
+ u (3:n+1:2, 3:n-1:2) &
+ u (2: n :2, 2:n-2:2) &
+ u (2: n :2, 4: n :2))
ENDDO
END
! }}}
! {{{ SUBROUTINE RESTR (n, fc, uf, ff) [TIMED]
SUBROUTINE RESTR (n, fc, uf, ff)
!--------------------------------------------------------------------
! Called by 'mgiup ()'.
! COMPUTATION OF THE DEFECT AND FINE-TO-COARSE TRANSFER (HALF-INJECTION)
!--------------------------------------------------------------------
IMPLICIT NONE
INTEGER n
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: fc
DOUBLE PRECISION, DIMENSION (2*n+1, 2*n+1) :: ff, uf
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: fc
!HPF$ ALIGN WITH * tmplf :: ff, uf
fc (2:n, 2:n) = 2 * (ff (3 : 2*n-1 : 2, 3 : 2*n-1 : 2) &
- 4 * uf (3 : 2*n-1 : 2, 3 : 2*n-1 : 2) &
+ uf (3 : 2*n-1 : 2, 2 : 2*n-2 : 2) &
+ uf (3 : 2*n-1 : 2, 4 : 2*n : 2) &
+ uf (2 : 2*n-2 : 2, 3 : 2*n-1 : 2) &
+ uf (4 : 2*n : 2, 3 : 2*n-1 : 2))
END
! }}}
! {{{ DOUBLE PRECISION FUNCTION DEFMX (n, u, f)
DOUBLE PRECISION FUNCTION DEFMX (n, u, f)
!-----------------------------------------------------------------
! Called from 'getdef ()'.
! COMPUTES THE MAXIMUM NORM OF THE DEFECT
!-----------------------------------------------------------------
! Mods
! ----
! -- Replaced parameter NP (= n+1) by 'n', and made corresponding
! changes to code and calls.
!-----------------------------------------------------------------
IMPLICIT NONE
INTEGER n
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: f, u
DOUBLE PRECISION hsqr
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: f, u
hsqr = DBLE (FLOAT (n * n))
defmx = MAXVAL (DABS (hsqr * (f (2: n , 2: n ) &
- 4.0D0 * u (2: n , 2: n ) &
+ u (1:n-1, 2: n ) &
+ u (3:n+1, 2: n ) &
+ u (2: n , 1:n-1) &
+ u (2: n , 3:n+1)) ))
END
! }}}
! {{{ DOUBLE PRECISION FUNCTION DIFMX (n, u)
DOUBLE PRECISION FUNCTION DIFMX (n, u)
!-----------------------------------------------------------------
! Called by 'getdif (m, sol)'.
! COMPUTES THE MAXIMUM NORM OF THE DIFFERENCE BETWEEN
! SOL (=SOLUTION OF THE BVP) AND THE VALUES IN u.
! (PARALLELISED BY INLINING CALLS TO `SOL', AND CONVERTING
! TO ARRAY SYNTAX).
!-----------------------------------------------------------------
! Mods
! ----
! -- Replaced parameter NP (= n+1) by 'n', and made corresponding
! changes to code and calls.
! -- Removed dummy function parameter 'sol' as it isn't used.
!-----------------------------------------------------------------
IMPLICIT NONE
INTEGER n
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: u
INTEGER i
DOUBLE PRECISION h
DOUBLE PRECISION, DIMENSION (n+1) :: x, y
DOUBLE PRECISION, DIMENSION (n+1, n+1) :: solval
INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: u
!HPF$ ALIGN WITH tmplc :: solval
!HPF$ ALIGN WITH tmplc (:, *) :: x
!HPF$ ALIGN WITH tmplc (*, :) :: y
h = 1.0D0 / n
DO i = 1, n+1
x (i) = (i - 1) * h
y (i) = (i - 1) * h
ENDDO
solval = SIN (3.0D0 * SPREAD (x, 2, n+1) + SPREAD (y, 1, n+1))
difmx = MAXVAL (DABS (solval - u))
END
! }}}
! }}}
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?