C  MAX_POOL — 2D max pooling, double-precision, multi-channel,
C  configurable window and stride. Non-overlapping or overlapping
C  windows supported via the S parameter.
C
C  Reference: standard CNN pooling layer; see e.g.
C  Goodfellow, Bengio, Courville, "Deep Learning" (2016), ch. 9.
C  Equivalent to cuDNN's CUDNN_POOLING_MAX with square window
C  and equal horizontal/vertical stride.
C
C  Hand-written reference for the Dark Factory's Phase 3 inference
C  kernel ladder, 2026-05-24. Public domain.
C
C  Mathematical definition (max pooling):
C    For every output position (IC, OH, OW):
C      Y(IC, OH, OW) = max over (KR, KS) of
C                        X(IC, (OH - 1) * S + KR,
C                              (OW - 1) * S + KS)
C  with KR in 1..KH, KS in 1..KW.
C
C  Shapes (column-major Fortran):
C    X is C x H  x W                       (input feature map)
C    Y is C x HO x WO                      (output feature map)
C  where HO = (H - KH) / S + 1,
C        WO = (W - KW) / S + 1.

      SUBROUTINE MAX_POOL
     +   (H, W, C, KH, KW, S, X, Y)
C  Inputs:
C    H, W      — input feature-map height and width (positive)
C    C         — channel count (positive)
C    KH, KW    — pooling window height and width (positive,
C                KH <= H, KW <= W)
C    S         — stride (positive)
C    X         — C x H x W input feature map
C  Output:
C    Y         — C x ((H - KH) / S + 1) x ((W - KW) / S + 1)
C                output feature map; each cell is the maximum
C                over its receptive-field window.
      INTEGER H, W, C, KH, KW, S
      DOUBLE PRECISION X(C, H, W)
      DOUBLE PRECISION Y(C, (H - KH) / S + 1, (W - KW) / S + 1)
      INTEGER IC, OH, OW, KR, KS, HO, WO
      DOUBLE PRECISION CUR, BEST
      HO = (H - KH) / S + 1
      WO = (W - KW) / S + 1
      DO 30 IC = 1, C
         DO 25 OH = 1, HO
            DO 20 OW = 1, WO
               BEST = X(IC, (OH - 1) * S + 1, (OW - 1) * S + 1)
               DO 15 KR = 1, KH
                  DO 10 KS = 1, KW
                     CUR = X(IC, (OH - 1) * S + KR,
     +                          (OW - 1) * S + KS)
                     IF (CUR .GT. BEST) BEST = CUR
   10             CONTINUE
   15          CONTINUE
               Y(IC, OH, OW) = BEST
   20       CONTINUE
   25    CONTINUE
   30 CONTINUE
      RETURN
      END