mg.h90

来自「shpf 1.9一个并行编译器」· H90 代码 · 共 957 行 · 第 1/3 页

! INITIALISATION ON FINEST LEVEL

  SELECT CASE (M)
    CASE (1);  CALL INIT1 (n1, H, U1, F1)
    CASE (2);  CALL INIT1 (n2, H, U2, F2)
    CASE (3);  CALL INIT1 (n3, H, U3, F3)
    CASE (4);  CALL INIT1 (n4, H, U4, F4)
    CASE (5);  CALL INIT1 (n5, H, U5, F5)
    CASE (6);  CALL INIT1 (n6, H, U6, F6)
    CASE (7);  CALL INIT1 (n7, H, U7, F7)
    CASE (8);  CALL INIT1 (n8, H, U8, F8)
    CASE (9);  CALL INIT1 (n9, H, U9, F9)
  END SELECT

! TRANSFER OF F AND U TO COARSER GRIDS

  DO L = 2, LIN
    SELECT CASE (M - L + 1)
      CASE (1);  CALL TRANS (n1, U1, F1, U2, F2)
      CASE (2);  CALL TRANS (n2, U2, F2, U3, F3)
      CASE (3);  CALL TRANS (n3, U3, F3, U4, F4)
      CASE (4);  CALL TRANS (n4, U4, F4, U5, F5)
      CASE (5);  CALL TRANS (n5, U5, F5, U6, F6)
      CASE (6);  CALL TRANS (n6, U6, F6, U7, F7)
      CASE (7);  CALL TRANS (n7, U7, F7, U8, F8)
      CASE (8);  CALL TRANS (n8, U8, F8, U9, F9)
    END SELECT
  ENDDO

! Zero BOUNDARY VALUES FOR MULTIGRID CORRECTIONS

  DO L = LIN + 1, M
    SELECT CASE (M - L + 1)
      CASE (1);  CALL PUTZB (n1, U1)
      CASE (2);  CALL PUTZB (n2, U2)
      CASE (3);  CALL PUTZB (n3, U3)
      CASE (4);  CALL PUTZB (n4, U4)
      CASE (5);  CALL PUTZB (n5, U5)
      CASE (6);  CALL PUTZB (n6, U6)
      CASE (7);  CALL PUTZB (n7, U7)
      CASE (8);  CALL PUTZB (n8, U8)
      CASE (9);  CALL PUTZB (n9, U9)
    END SELECT
  ENDDO

END
! }}}

! {{{  SUBROUTINE MGI(LEV, M, NY1, NY2, IGAM)  [TIMED]
SUBROUTINE MGI(LEV, M, NY1, NY2, IGAM)
  IMPLICIT NONE
  INTEGER LEV, M, NY1, NY2, IGAM
  INTEGER K, ICGAM (10), IGAML (10)
  LOGICAL IZ
  INCLUDE 'grids.inc'

!--ONE MULTIGRID ITERATION STEP

  DO K = LEV, M
    IF ((K == LEV) .OR. (K == M)) THEN
      IGAML (K) = 1
    ELSE
      IGAML (K) = IGAM
    ENDIF
    ICGAM(K) = 0
  ENDDO
  K = LEV

  DO WHILE (ICGAM (K) < IGAML (K))

!--``UP'' -- i.e. fine-to-coarse; level i = M -> 1
! using RELAX (i) then RESTR (i -> i-1)

    DO WHILE (K < M)
      IZ = K .GT. LEV .AND. ICGAM(K) .EQ. 0
      
      SELECT CASE (M - K)
        CASE (1);  CALL MGIUP (n1, F1, U2, F2, NY1, IZ)
        CASE (2);  CALL MGIUP (n2, F2, U3, F3, NY1, IZ)
        CASE (3);  CALL MGIUP (n3, F3, U4, F4, NY1, IZ)
        CASE (4);  CALL MGIUP (n4, F4, U5, F5, NY1, IZ)
        CASE (5);  CALL MGIUP (n5, F5, U6, F6, NY1, IZ)
        CASE (6);  CALL MGIUP (n6, F6, U7, F7, NY1, IZ)
        CASE (7);  CALL MGIUP (n7, F7, U8, F8, NY1, IZ)
        CASE (8);  CALL MGIUP (n8, F8, U9, F9, NY1, IZ)
      END SELECT
      
      ICGAM(K) = ICGAM(K) + 1
      K = K + 1
    ENDDO

!--EXACT SOLUTION ON COARSEST GRID (K = M)

    CALL RELAX(1, (M > LEV .AND. ICGAM(M)==0), 2, U1, F1)

    ICGAM (M) = ICGAM (M) + 1

!--``DOWN'' -- i.e. coarse-to-fine; level i = 1 -> M,
! INT2A (i -> i+1) then RELAX (i+1)

    DO WHILE ((K > LEV) .AND. (ICGAM (K) == IGAML (K)))
      ICGAM (K) = 0
      K = K - 1
      SELECT CASE (M - K)
        CASE (1);  CALL MGIDN (n1, U2, F2, U1, NY2)
        CASE (2);  CALL MGIDN (n2, U3, F3, U2, NY2)
        CASE (3);  CALL MGIDN (n3, U4, F4, U3, NY2)
        CASE (4);  CALL MGIDN (n4, U5, F5, U4, NY2)
        CASE (5);  CALL MGIDN (n5, U6, F6, U5, NY2)
        CASE (6);  CALL MGIDN (n6, U7, F7, U6, NY2)
        CASE (7);  CALL MGIDN (n7, U8, F8, U7, NY2)
        CASE (8);  CALL MGIDN (n8, U9, F9, U8, NY2)
      END SELECT
    ENDDO
  ENDDO
END
! }}}
! {{{  DOUBLE PRECISION FUNCTION GETDEF (m)
DOUBLE PRECISION FUNCTION GETDEF (m)
!---------------------------------------------------------
! Called by main program and 'mg00d ()'.
! EXTRACT MAXIMUM NORM OF THE DEFECT FROM STACK
!---------------------------------------------------------
  INTEGER  m
  DOUBLE PRECISION  defmx
  INCLUDE 'grids.inc'

  SELECT CASE (m)
    CASE (1);  GETDEF = DEFMX (n1, U1, F1)
    CASE (2);  GETDEF = DEFMX (n2, U2, F2)
    CASE (3);  GETDEF = DEFMX (n3, U3, F3)
    CASE (4);  GETDEF = DEFMX (n4, U4, F4)
    CASE (5);  GETDEF = DEFMX (n5, U5, F5)
    CASE (6);  GETDEF = DEFMX (n6, U6, F6)
    CASE (7);  GETDEF = DEFMX (n7, U7, F7)
    CASE (8);  GETDEF = DEFMX (n8, U8, F8)
    CASE (9);  GETDEF = DEFMX (n9, U9, F9)
  END SELECT
END
! }}}
! {{{  DOUBLE PRECISION FUNCTION GETDIF (m)
DOUBLE PRECISION FUNCTION GETDIF (m)
!-----------------------------------------------------------------
!  Called by the main program.
!  EXTRACT MAXIMUM NORM OF THE DIFFERENCE FROM STACK
!-----------------------------------------------------------------
! Mods
! ----
! -- Converted computed GOTO to CASE construct.
! -- Deleted dummy function argument 'sol', which isn't used.
!-----------------------------------------------------------------
  IMPLICIT NONE
  INTEGER  m
  DOUBLE PRECISION DIFMX
  INCLUDE 'grids.inc'

  SELECT CASE (m)
    CASE (1);  GETDIF = DIFMX (n1, U1)
    CASE (2);  GETDIF = DIFMX (n2, U2)
    CASE (3);  GETDIF = DIFMX (n3, U3)
    CASE (4);  GETDIF = DIFMX (n4, U4)
    CASE (5);  GETDIF = DIFMX (n5, U5)
    CASE (6);  GETDIF = DIFMX (n6, U6)
    CASE (7);  GETDIF = DIFMX (n7, U7)
    CASE (8);  GETDIF = DIFMX (n8, U8)
    CASE (9);  GETDIF = DIFMX (n9, U9)
  END SELECT

END
! }}}

! {{{  SUBROUTINE FMGDN (n, uf, ff, uc)  [start of level]
SUBROUTINE FMGDN (n, uf, ff, uc)
!-----------------------------------------------------------------
! Called by 'mg00d ()' to perform one 'down' step in the full
! multigrid.
!-----------------------------------------------------------------
  IMPLICIT NONE
  INTEGER  n
  DOUBLE PRECISION, DIMENSION (n+1,   n+1)   ::  uc
  DOUBLE PRECISION, DIMENSION (2*n+1, 2*n+1) ::  ff, uf
  
  INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: uc
!HPF$ ALIGN WITH * tmplf :: uf, ff


  CALL INT4  (n, uc, uf, ff)
  CALL PUTZB (n, uc)
END
! }}}
! {{{  SUBROUTINE MGIUP (n, fc, uf, ff, NY1, IZ)  [TIMED]
SUBROUTINE MGIUP (n, fc, uf, ff, NY1, IZ)
!-----------------------------------------------------------------
!  Called by 'mgi ()' to perform one 'up' step in a multigrid iteration.
!  RELAXATIONS BEFORE COARSE-GRID CORRECTION...
!-----------------------------------------------------------------
! Mods
! ----
! -- Replaced parameters NPC (= n+1) & NPF (= 2*n+1) by 'n',
! and made corresponding changes to code and calls.
!-----------------------------------------------------------------
  IMPLICIT NONE
  INTEGER  n, NY1
  LOGICAL  IZ
  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 :: uf, ff


  CALL RELAX (NY1, IZ, (2*n), uf, ff)

!  RESIDUAL TRANSFER TO NEXT COARSER GRID...

  CALL RESTR (n, fc, uf, ff)

END
! }}}
! {{{  SUBROUTINE MGIDN (n, uf, ff, uc, NY2)      [TIMED]
SUBROUTINE MGIDN (n, uf, ff, uc, NY2)
!-----------------------------------------------------------------
!  Called by 'mgi ()' to perform one 'down' step in a multigrid iteration.
!--LINEAR INTERPOLATION TO NEXT FINER GRID
!-----------------------------------------------------------------
! Mods
! ----
! -- Replaced parameters NPF (= 2*n+1) & NPC (= n+1) by 'n',
! and made corresponding changes to code and calls.
!-----------------------------------------------------------------
  IMPLICIT NONE
  INTEGER  n, NY2
  DOUBLE PRECISION, DIMENSION (n+1,   n+1)   ::  uc
  DOUBLE PRECISION, DIMENSION (2*n+1, 2*n+1) ::  ff, uf

  INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: uc
!HPF$ ALIGN WITH * tmplf :: uf, ff


  CALL INT2A (n, uc, uf)

!--RELAXATION AFTER CGC

  CALL RELAX (NY2, .FALSE., (2*n), uf, ff)
END
! }}}

!------------------------------------------------!
! Initialisation routines, called from 'init ()' !
!------------------------------------------------!
! {{{  SUBROUTINE INIT1 (n, h, u, f)  [start of level]
SUBROUTINE INIT1 (n, h, u, f)
! {{{  *** Mods ***
!----------------------------------------------------------------
! Mods
! ----
! -- Replaced parameter NP (= n+1) by 'n'.  Made corresponding
! replacements within this function and in calls to it.
!
! -- Deleted dummy function parameters F & G, since they aren't used.
! Made corresponding change to calls.
!
! -- N.B.  The RHS's of the 1st 2 eqns below look suspicious.
! Are they correct?
!        U (  1    , 2 : n) = SIN (3.0D0 * 0.0D0 + y (2 : n))
!        U ( n+1   , 2 : n) = SIN (3.0D0 * 1.0D0 + y (2 : n))
!        U (1 : n+1,   1  ) = SIN (3.0D0 * X + 0.0D0)
!        U (1 : n+1,  n+1 ) = SIN (3.0D0 * X + 1.0D0)
!----------------------------------------------------------------
! }}}
!-------------------------------------------------------------------
! Called by 'init ()'.
!     COMPUTES INITIAL VALUES ON THE FINEST GRID
!     (PARALLELISED BY INLINING CALLS TO `F' AND `G', AND CONVERTING
!     TO ARRAY SYNTAX).
!-------------------------------------------------------------------
  IMPLICIT NONE
  INTEGER  n
  DOUBLE PRECISION  :: h
  DOUBLE PRECISION, DIMENSION (n+1, n+1) :: u, f

  INTEGER  i
  DOUBLE PRECISION, DIMENSION (n+1, n+1) :: fval
  DOUBLE PRECISION, DIMENSION (n+1)      :: x, y

  INCLUDE 'tmpl.inc'
!HPF$ ALIGN WITH * tmplc :: u, f
!HPF$ ALIGN WITH   tmplc :: fval
!HPF$ ALIGN WITH   tmplc (:, *) :: x
!HPF$ ALIGN WITH   tmplc (*, :) :: y


  DO i = 1, n+1
    x (i) = (i - 1) * h
    y (i) = (i - 1) * h
  ENDDO

  u (2:n, 2:n) = 0.0D0
  u ( 1 , 2:n) = SIN (y (2:n))
  u (n+1, 2:n) = SIN (3.0D0 + y (2:n))
  u ( : ,  1 ) = SIN (3.0D0 * x)
  u ( : , n+1) = SIN (3.0D0 * x + 1.0D0)

  fval = 10.0D0 * SIN (3.0D0 * SPREAD(x, 2, n+1) + SPREAD(y, 1, n+1))
  f (2:n,  2:n) = h*h * fval (2:n,  2:n)
END
! }}}
! {{{  SUBROUTINE TRANS (n, uc, fc, uf, ff)  [start of level]
SUBROUTINE TRANS (n, uc, fc, uf, ff)
! {{{  *** Mods ***
!------------------------------------------------------------------
! Mods
! ----
! -- Originally this had args
!
!      SUBROUTINE TRANS (NPC, NPF, UC, fc, UF, FF)