mg.h90

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

! {{{  mg.h90
! {{{  PROGRAM PDE2 
PROGRAM PDE2
! {{{  Description
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     MAIN PROGRAM FOR MG BENCHMARK - SIMD/VECTOR VERSION
!
!     SUPPLIED BY   MAX LEMKE
!                   SUPRENUM GMBH
!                   HOHE STRASSE 73
!                   D-5300 BONN 1
!
!     VERSION FEBRUARY 15, 1989
!
!     REFERENCES:
!               - K. STUEBEN, U. TROTTENBERG: MULTIGRID METHODS:
!                 FUNDAMENTAL ALGORITHMS, MODEL PROBLEM ANALYSIS
!                 AND APPLICATIONS, DECEMBER 84, GMD-STUDIEN NR. 96,
!                 GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG MBH,
!                 BONN
!
!               - M. LEMKE: EXPERIMENTS WITH A VECTORIZED MULTIGRID
!                 POISSON SOLVER ON THE CDC CYBER 205, CRAY X-MP AND
!                 FUJITSU VP 200;
!                 NOVEMBER 85, ARBEITSPAPIERE DER GMD, NR. 179,
!                 GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG MBH,
!                 BONN
!
!               - M. LEMKE: ERFAHRUNGEN MIT MEHRGITTERVERFAHREN FUER
!                 HELMHOLTZ-AEHNLICHE PROBLEME AUF VEKTORRECHNERN
!                 UND MEHRPROZESSOR-VEKTORRECHNERN
!                 DECEMBER 87, ARBEITSPAPIERE DER GMD, NR. 278,
!                 GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG MBH,
!                 BONN
!
!     AUTHORS: FIRST VERSION:
!                 K. WITSCH, H. FOERSTER (1984, GMD)
!              VECTORIZED (AND MULTITASKING) VERSION:
!                 M. LEMKE (1984 - 1987, GMD)
!              GENESIS BENCHMARK VERSION:
!                 M. LEMKE (1989, SUPRENUM GMBH)
!              SUBSET FORTRAN 90 REWRITE:
!                 D.B. CARPENTER (1993, UNIV. OF SOUTHAMPTON)
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!     THE TIMER (DOUBLE PRECISION IN SECONDS) IS CALLED IN MG00D.
!     THE MACHINE DEPENDENT TIMING ROUTINE HAS TO BE INCLUDED.
!
!     ATTENTION: ONLY THE MG KERNEL IS VECTORIZED. SOME OTHER ROUTINES
!                ARE PURELY SEQUENTIAL. THEY NEED NOT BE VECTORIZED
!                BECAUSE THEY ARE NOT INCLUDED IN TIMING.
!
!                THE FMG OPTION IS NOT AVAILABLE IN THIS VERSION
!
!     FOR NUMERICAL VERIFICATION:
!        THE AVERAGE CONVERGENCE RATE SHOLD BE SOMEWHERE IN BETWEEN
!        0.02 AND 0.07.
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! }}}
! {{{  decls
  IMPLICIT NONE
  DOUBLE PRECISION  GETDIF, GETDEF, DIF, DEF
  INTEGER  M, N, NP, NY1, NY2, NFMG, ncycle, IGAM
  INTEGER NOP, LEV

  DOUBLE PRECISION DEF1, DEF2, CONVR
  DOUBLE PRECISION  RMFLOP
  DOUBLE PRECISION SEC2
  COMMON /TT/ SEC2
  COMMON /DEFECT/ DEF1, DEF2

  INCLUDE "runparms.inc"
  INCLUDE "procs.inc"
! }}}

  SEC2=1.
! CALL HEADER(6)
  PRINT *, ' ************* MULTIGRID BENCHMARK ******************* '

! LOOP OVER THE FINEST GRIDS IN ORDER TO GET RESULTS FOR DIFFERENT
! PROBLEM SIZES

  ncycle = niter        ! initial value = 10  (originally 60)
  DO M = mmin, MMAX     ! 9,9.  Originally, 3,8

    NY1 = 2
    NY2 = 1
    NFMG = 0
    IGAM = 1
  
    CALL MG00D (M, NY1, NY2, NFMG, NCYCLE, IGAM)
  
    N = 2**M
    NP = N + 1
    DIF = GETDIF (M)
    DEF = GETDEF (M)
    CONVR=(DEF/DEF1)**(0.1)
    WRITE (6, 9100) NP, DIF, DEF, CONVR
  
    NOP = 1
    DO LEV = 2, M
      NOP = NOP + (2**LEV-1)*(2**LEV-1)*(2.75+9.)           &
                        + (2**(LEV-1)-1)*(2**(LEV-1)-1)*7   &
                        + (2**(LEV-1)-1)*(2**(LEV-1)-1)*6
    ENDDO
    NOP = NOP + (2**M-1)*(2**M-1)*1.75
  
    RMFLOP = (ncycle * NOP) / (10**6 * SEC2)
    WRITE(6, 900) nproc0, nproc1, (SEC2 / NCYCLE), ncycle, RMFLOP, nop
    PRINT *, ' +++++++++++++++++++++++++++++++++++++++++++++++++++++ '

!    NCYCLE = MAX(NCYCLE/4,10)
  ENDDO
  
  STOP
 9100 FORMAT (/' NUMBER OF POINTS PER DIRECTION = ', I12 /     &
               ' MAXIMUM NORM OF THE ERROR  =     ', D12.4 /   &
               ' MAXIMUM NORM OF THE DEFECT =     ', D12.4/    &
               ' AVERAGE CONV. RATE             = ', D12.4/)
  900 FORMAT(/ '(', i2, ',', i2, ') processor array' /         &
               f8.3, ' sec / cycle  (average of ', i3, ' cycles).'/ &
               f9.5, ' Mflops      (', i15, ' operations).'/)
END

! }}}
! {{{  DOUBLE PRECISION FUNCTION F(X,Y) -- could be elemental!
DOUBLE PRECISION FUNCTION F(X,Y)
  IMPLICIT NONE
  DOUBLE PRECISION X, Y

  F = 10.0D0 * SIN (3.0D0 * X + Y)
END
! }}}
! {{{  DOUBLE PRECISION FUNCTION G (X, Y) -- could be elemental!
DOUBLE PRECISION FUNCTION G (X, Y)
  IMPLICIT NONE
  DOUBLE PRECISION X, Y

  G = SIN (3.0D0 * X + Y)
END
! }}}
! {{{  SUBROUTINE MG00D (M, NY1, NY2, NFMG, NCYCLE, IGAM) [once per level]
SUBROUTINE MG00D (M, NY1, NY2, NFMG, NCYCLE, IGAM)

! {{{  Description
!
!     MULTIGRID MODULE FOR THE FAST SOLUTION OF POISSON'S EQUATION
!     WITH DIRICHLET BOUNDARY CONDITIONS ON THE UNIT SQUARE
!
!          DIFFERENTIAL EQUATION :   -DELTA U(X,Y) = F(X,Y)
!          BOUNDARY CONDITION :             U(X,Y) = G(X,Y)
!
!     CYCLE STRUCTURE:
!
!                COARSENING: H, 2H, 4H, ...
!                DIFFERENCE OPERATORS: USUAL 5-POINT STARS ON ALL GRIDS
!                RELAXATION: RED-BLACK
!                FINE-TO-COARSE: HALF INJECTION
!                COARSE-TO-FINE: BILINEAR INTERPOLATION
!
!     FULL MULTIGRID MODE:
!
!                FULL MULTIGRID INTERPOLATION USES GRID EQUATION
!                (4-TH ORDER)
!
!     INPUT:
!
!        M       NUMBER OF GRIDS (0 < M < 11). FOR GIVEN M:
!                NUMBER OF P0INTS ON THE FINEST GRID (INCLUDING BOUNDARY
!                POINTS) = 2**M + 1 IN BOTH DIRECTIONS
!        NY1     NUMBER OF RELAXATIONS BEFORE COARSE-GRID CORRECTION
!                (NY1 > 0)
!        NY2     NUMBER OF RELAXATIONS AFTER COARSE-GRID CORRECTION
!                (NY2 > 0)
!        NCYCLE  NUMBER OF MULTIGRID ITERATIONS (NCYCLE > 0)
!        NFMG    .EQ. 0 : NCYCLE MULTIGRID ITERATIONS ARE PERFORMED
!                .NE. 0 : FMG-VERSION IS PERFORMED PLUS NCYCLE-1
!                         ADDITIONAL MG ITERATIONS AFTERWARDS
!        IGAM    TYPE OF CYCLING (IGAM > 0).
!                E.G.: IGAM=1 FOR V-CYCLES, IGAM=2 FOR W-CYCLES
!        W       DOUBLE PRECISION WORK ARRAY OF DIMENSION IDIM
!        IDIM    DIMENSION OF W. APPROXIMATELY: IDIM > 2.8*4**M
!
!     EXTERNALS:
!
!        F       DOUBLE PRECISION FUNCTION F(X,Y), RIGHT HAND SIDE OF
!                THE DIFFERENTIAL EQUATION
!        G       DOUBLE PRECISION FUNCTION G(X,Y), BOUNDARY VALUES
!
!     REMARK:    GRID #1 AND #M ARE THE FINEST AND COARSEST GRID USED,
!                RESPECTIVELY
!
!     OUTPUT:
!
!        IER     ERROR INDICATOR
!                = 0  NO ERRORS
!                = 1  INSUFFICIENT MEMORY , I.E. IDIM TOO SMALL.
!                     IN THIS CASE IDIM IS USED AS OUTPUT PARAMETER
!                     TO SHOW THE MINIMAL DIMENSION
!                = 2  M, NY1, NY2, NCYCLE OR IGAM WRONG
!        JDIM    ONLY IN CASE IER=1 : MINIMAL LENGTH OF W
!        INITF   CF. DESCRIPTION OF W
!        W       W CONTAINS THE DISCRETE APPROXIMATION TO THE GIVEN
!                BOUNDARY VALUE PROBLEM ON THE FINEST GRID. THE GRID
!                VALUES ARE STORED ROWWISE FROM LEFT TO RIGHT AND
!                FROM BOTTOM TO TOP.
!                I.E. THE GRID VALUE CORRESPONDING TO THE GRID POINT
!                (XI,YJ) :
!
!                     XI = (I-1)*H,  YJ = (J-1)*H   (H=1/N, N=2**M)
!
!                IS STORED AT
!
!                        W((J-1)*(N+1)+I)     (0 < I,J < N+2)
!
!                THE CORRESPONDING VALUES OF THE RIGHT HAND SIDE ARE
!                STORED IN THE SAME MANNER AT
!
!                        W((J-1)*(N+1)+I+INITF).
!
!                THE REMAINING STORAGE CONTAINS COARSE GRID VALUES.
!
!     REMARK:    THIS IS A SPECIALIZED PROGRAM FROM THE MG00 PROGRAM
!                COLLECTION. SEE: FOERSTER, WITSCH IN
!                'ON EFFICIENT MULTIGRID SOFTWARE FOR ELLIPTIC
!                PROBLEMS ON RECTANGULAR DOMAINS', MATH. COMPUT.
!                SIMULATION XXIII, PP. 293-298, 1981
!
!
! }}}
  IMPLICIT NONE
  INTEGER M, NY1, NY2, NFMG, NCYCLE, IGAM
  INTEGER LIN, IT, LEV
  DOUBLE PRECISION  DEF1, DEF2
  DOUBLE PRECISION  SEC1, SEC2
  COMMON /TT/ SEC2
  COMMON /DEFECT/ DEF1,DEF2
  DOUBLE PRECISION GETDEF
  INCLUDE 'grids.inc'
  

! LIN = LEVEL FROM WHICH INITIAL GUESS PROPAGATED DOWN

  IF (NFMG == 0) THEN
    LIN = 1
  ELSE
    LIN = M
  ENDIF

!       SET UP ALL GRID VALUES NEEDED

  CALL INIT (M, LIN)

! (FULL) MULTIGRID PROCEDURE...
! INITIAL GUESS FROM COARSEST GRID

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

! MULTIGRID CYCLES

  DEF1 = GETDEF (M)
  
!----> TIMING
  
  CALL TIMER (SEC1)
  DO IT = 1, NCYCLE
    CALL MGI (1, M, NY1, NY2, IGAM)
  ENDDO
  CALL TIMER (SEC2)
  SEC2 = SEC2 - SEC1
  
!----> TIMING
  
END
! }}}

! {{{  SUBROUTINE INIT (M, LIN)   [start of level]
SUBROUTINE INIT (M, LIN)
! {{{  *** Mods ***
!---------------------------------------------------------------
! Mods
! ----
! -- Changed computed GOTOs, e.g.:
!        GOTO (101, 102, 103, 104, 105, 106, 107, 108) M
! to equivalent 'SELECT CASE' constructs.
!
! -- Deleted dummy function parameters F & G, since they aren't used.
! Made corresponding change to calls.
!---------------------------------------------------------------
! }}}
!---------------------------------------------------------------
! Called by 'mg00d ()'.
!     COMPUTES INITIAL VALUES ON ALL GRIDS.
!     LIN = LEVEL FROM WHICH INITIAL GUESS PROPAGATED DOWN
!---------------------------------------------------------------
  IMPLICIT NONE
  INTEGER M, LIN

  INTEGER L
  DOUBLE PRECISION H
  INCLUDE 'grids.inc'

  H = 1.0D0 / DBLE (FLOAT (2**M))