From d9b4070e59d52b7af085fa456916684f51e6157b Mon Sep 17 00:00:00 2001
From: Riccardo Orsi <riccardoorsi@icloud.com>
Date: Sun, 16 Mar 2025 21:17:32 +0100
Subject: [PATCH] Integrated new features into environment

---
 CMakeLists.txt                             |   8 ++
 src/nf.f90                                 |   4 +
 src/nf/nf_avgpool1d_layer.f90              |  66 +++++++++
 src/nf/nf_avgpool1d_layer_submodule.f90    |  87 ++++++++++++
 src/nf/nf_avgpool2d_layer.f90              |  66 +++++++++
 src/nf/nf_avgpool2d_layer_submodule.f90    |  94 +++++++++++++
 src/nf/nf_avgpool3d_layer.f90              |  67 +++++++++
 src/nf/nf_avgpool3d_layer_submodule.f90    | 102 ++++++++++++++
 src/nf/nf_layer_constructors.f90           |  92 +++++++++++++
 src/nf/nf_layer_constructors_submodule.f90 | 119 ++++++++++++++++
 src/nf/nf_layer_submodule.f90              | 149 +++++++++++++++++++++
 src/nf/nf_maxpool3d_layer.f90              |  79 +++++++++++
 src/nf/nf_maxpool3d_layer_submodule.f90    | 114 ++++++++++++++++
 13 files changed, 1047 insertions(+)
 create mode 100644 src/nf/nf_avgpool1d_layer.f90
 create mode 100644 src/nf/nf_avgpool1d_layer_submodule.f90
 create mode 100644 src/nf/nf_avgpool2d_layer.f90
 create mode 100644 src/nf/nf_avgpool2d_layer_submodule.f90
 create mode 100644 src/nf/nf_avgpool3d_layer.f90
 create mode 100644 src/nf/nf_avgpool3d_layer_submodule.f90
 create mode 100644 src/nf/nf_maxpool3d_layer.f90
 create mode 100644 src/nf/nf_maxpool3d_layer_submodule.f90

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 667b5d73..c892a494 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,6 +17,12 @@ include(cmake/compilers.cmake)
 add_library(neural-fortran
   src/nf.f90
   src/nf/nf_activation.f90
+  src/nf/nf_avgpool1d_layer.f90
+  src/nf/nf_avgpool1d_layer_submodule.f90
+  src/nf/nf_avgpool2d_layer.f90
+  src/nf/nf_avgpool2d_layer_submodule.f90
+  src/nf/nf_avgpool3d_layer.f90
+  src/nf/nf_avgpool3d_layer_submodule.f90
   src/nf/nf_base_layer.f90
   src/nf/nf_conv1d_layer.f90
   src/nf/nf_conv1d_layer_submodule.f90
@@ -55,6 +61,8 @@ add_library(neural-fortran
   src/nf/nf_maxpool1d_layer_submodule.f90
   src/nf/nf_maxpool2d_layer.f90
   src/nf/nf_maxpool2d_layer_submodule.f90
+  src/nf/nf_maxpool3d_layer.f90
+  src/nf/nf_maxpool3d_layer_submodule.f90
   src/nf/nf_metrics.f90
   src/nf/nf_multihead_attention.f90
   src/nf/nf_multihead_attention_submodule.f90
diff --git a/src/nf.f90 b/src/nf.f90
index f644826d..f04bd2f0 100644
--- a/src/nf.f90
+++ b/src/nf.f90
@@ -3,6 +3,9 @@ module nf
   use nf_datasets_mnist, only: label_digits, load_mnist
   use nf_layer, only: layer
   use nf_layer_constructors, only: &
+    avgpool1d, &
+    avgpool2d, &
+    avgpool3d, &
     conv1d, &
     conv2d, &
     dense, &
@@ -15,6 +18,7 @@ module nf
     locally_connected1d, &
     maxpool1d, &
     maxpool2d, &
+    maxpool3d, &
     reshape, &
     self_attention
   use nf_loss, only: mse, quadratic
diff --git a/src/nf/nf_avgpool1d_layer.f90 b/src/nf/nf_avgpool1d_layer.f90
new file mode 100644
index 00000000..631f4895
--- /dev/null
+++ b/src/nf/nf_avgpool1d_layer.f90
@@ -0,0 +1,66 @@
+module nf_avgpool1d_layer
+    !! This module provides the 1-d average pooling layer.
+  
+    use nf_base_layer, only: base_layer
+    implicit none
+  
+    private
+    public :: avgpool1d_layer
+  
+    type, extends(base_layer) :: avgpool1d_layer
+        integer :: channels
+        integer :: width      ! Length of the input along the pooling dimension
+        integer :: pool_size
+        integer :: stride
+  
+        ! Gradient for the input (same shape as the input).
+        real, allocatable :: gradient(:,:)
+        ! Output after pooling (dimensions: (channels, new_width)).
+        real, allocatable :: output(:,:)
+    contains
+        procedure :: init
+        procedure :: forward
+        procedure :: backward
+    end type avgpool1d_layer
+  
+    interface avgpool1d_layer
+        pure module function avgpool1d_layer_cons(pool_size, stride) result(res)
+            !! `avgpool1d` constructor function.
+            integer, intent(in) :: pool_size
+                !! Width of the pooling window.
+            integer, intent(in) :: stride
+                !! Stride of the pooling window.
+            type(avgpool1d_layer) :: res
+        end function avgpool1d_layer_cons
+    end interface avgpool1d_layer
+  
+    interface
+        module subroutine init(self, input_shape)
+            !! Initialize the `avgpool1d` layer instance with an input shape.
+            class(avgpool1d_layer), intent(in out) :: self
+                !! `avgpool1d_layer` instance.
+            integer, intent(in) :: input_shape(:)
+                !! Array shape of the input layer, expected as (channels, width).
+        end subroutine init
+  
+        pure module subroutine forward(self, input)
+            !! Run a forward pass of the `avgpool1d` layer.
+            class(avgpool1d_layer), intent(in out) :: self
+                !! `avgpool1d_layer` instance.
+            real, intent(in) :: input(:,:)
+                !! Input data (output of the previous layer), with shape (channels, width).
+        end subroutine forward
+  
+        pure module subroutine backward(self, input, gradient)
+            !! Run a backward pass of the `avgpool1d` layer.
+            class(avgpool1d_layer), intent(in out) :: self
+                !! `avgpool1d_layer` instance.
+            real, intent(in) :: input(:,:)
+                !! Input data (output of the previous layer).
+            real, intent(in) :: gradient(:,:)
+                !! Gradient from the downstream layer, with shape (channels, pooled width).
+        end subroutine backward
+    end interface
+  
+  end module nf_avgpool1d_layer
+  
\ No newline at end of file
diff --git a/src/nf/nf_avgpool1d_layer_submodule.f90 b/src/nf/nf_avgpool1d_layer_submodule.f90
new file mode 100644
index 00000000..8e022fbf
--- /dev/null
+++ b/src/nf/nf_avgpool1d_layer_submodule.f90
@@ -0,0 +1,87 @@
+submodule(nf_avgpool1d_layer) nf_avgpool1d_layer_submodule
+  implicit none
+
+contains
+
+  pure module function avgpool1d_layer_cons(pool_size, stride) result(res)
+    implicit none
+    integer, intent(in) :: pool_size
+    integer, intent(in) :: stride
+    type(avgpool1d_layer) :: res
+
+    res % pool_size = pool_size
+    res % stride    = stride
+  end function avgpool1d_layer_cons
+
+
+  module subroutine init(self, input_shape)
+    implicit none
+    class(avgpool1d_layer), intent(in out) :: self
+    integer, intent(in) :: input_shape(:)
+    ! input_shape is expected to be (channels, width)
+
+    self % channels = input_shape(1)
+    ! The new width is the integer division of the input width by the stride.
+    self % width    = input_shape(2) / self % stride
+
+    ! Allocate the gradient array corresponding to the input dimensions.
+    allocate(self % gradient(input_shape(1), input_shape(2)))
+    self % gradient = 0
+
+    ! Allocate the output array (after pooling).
+    allocate(self % output(self % channels, self % width))
+    self % output = 0
+  end subroutine init
+
+
+  pure module subroutine forward(self, input)
+    implicit none
+    class(avgpool1d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:)
+    integer :: input_width
+    integer :: i, n
+    integer :: ii, iend
+    integer :: iextent
+
+    input_width = size(input, dim=2)
+    ! Ensure we only process complete pooling regions.
+    iextent = input_width - mod(input_width, self % stride)
+
+    ! Loop over the input with a step size equal to the stride and over all channels.
+    do concurrent (i = 1:iextent: self % stride, n = 1:self % channels)
+      ! Compute the index in the pooled (output) array.
+      ii = (i - 1) / self % stride + 1
+      ! Determine the ending index of the current pooling region.
+      iend = min(i + self % pool_size - 1, input_width)
+
+      ! Compute the average over the pooling region.
+      self % output(n, ii) = sum(input(n, i:iend)) / (iend - i + 1)
+    end do
+  end subroutine forward
+
+
+  pure module subroutine backward(self, input, gradient)
+    implicit none
+    class(avgpool1d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:)
+    real, intent(in) :: gradient(:,:)
+    integer :: channels, pooled_width
+    integer :: i, n, j, istart, iend
+    real :: scale_factor
+
+    channels    = size(gradient, dim=1)
+    pooled_width = size(gradient, dim=2)
+
+    ! The gradient for average pooling is distributed evenly over the pooling window.
+    do concurrent (n = 1:channels, i = 1:pooled_width)
+      istart = (i - 1) * self % stride + 1
+      iend = min(istart + self % pool_size - 1, size(input, dim=2))
+      scale_factor = 1.0 / (iend - istart + 1)
+
+      do j = istart, iend
+        self % gradient(n, j) = gradient(n, i) * scale_factor
+      end do
+    end do
+  end subroutine backward
+
+end submodule nf_avgpool1d_layer_submodule
\ No newline at end of file
diff --git a/src/nf/nf_avgpool2d_layer.f90 b/src/nf/nf_avgpool2d_layer.f90
new file mode 100644
index 00000000..fcf8b970
--- /dev/null
+++ b/src/nf/nf_avgpool2d_layer.f90
@@ -0,0 +1,66 @@
+module nf_avgpool2d_layer
+    !! This module provides the 2-d average pooling layer.
+  
+    use nf_base_layer, only: base_layer
+    implicit none
+  
+    private
+    public :: avgpool2d_layer
+  
+    type, extends(base_layer) :: avgpool2d_layer
+        integer :: channels
+        integer :: height      ! Height of the input
+        integer :: width       ! Width of the input
+        integer :: pool_size  ! Pooling window size (height, width)
+        integer :: stride    ! Stride (height, width)
+  
+        ! Gradient for the input (same shape as the input: channels, height, width).
+        real, allocatable :: gradient(:,:,:)
+        ! Output after pooling (dimensions: (channels, new_height, new_width)).
+        real, allocatable :: output(:,:,:)
+    contains
+        procedure :: init
+        procedure :: forward
+        procedure :: backward
+    end type avgpool2d_layer
+  
+    interface avgpool2d_layer
+        pure module function avgpool2d_layer_cons(pool_size, stride) result(res)
+            !! `avgpool2d` constructor function.
+            integer, intent(in) :: pool_size
+                !! Pooling window size (height, width).
+            integer, intent(in) :: stride
+                !! Stride (height, width).
+            type(avgpool2d_layer) :: res
+        end function avgpool2d_layer_cons
+    end interface avgpool2d_layer
+  
+    interface
+        module subroutine init(self, input_shape)
+            !! Initialize the `avgpool2d` layer instance with an input shape.
+            class(avgpool2d_layer), intent(in out) :: self
+                !! `avgpool2d_layer` instance.
+            integer, intent(in) :: input_shape(:)
+                !! Array shape of the input layer, expected as (channels, height, width).
+        end subroutine init
+  
+        pure module subroutine forward(self, input)
+            !! Run a forward pass of the `avgpool2d` layer.
+            class(avgpool2d_layer), intent(in out) :: self
+                !! `avgpool2d_layer` instance.
+            real, intent(in) :: input(:,:,:)
+                !! Input data (output of the previous layer), with shape (channels, height, width).
+        end subroutine forward
+  
+        pure module subroutine backward(self, input, gradient)
+            !! Run a backward pass of the `avgpool2d` layer.
+            class(avgpool2d_layer), intent(in out) :: self
+                !! `avgpool2d_layer` instance.
+            real, intent(in) :: input(:,:,:)
+                !! Input data (output of the previous layer).
+            real, intent(in) :: gradient(:,:,:)
+                !! Gradient from the downstream layer, with shape (channels, pooled_height, pooled_width).
+        end subroutine backward
+    end interface
+  
+end module nf_avgpool2d_layer
diff --git a/src/nf/nf_avgpool2d_layer_submodule.f90 b/src/nf/nf_avgpool2d_layer_submodule.f90
new file mode 100644
index 00000000..3029b826
--- /dev/null
+++ b/src/nf/nf_avgpool2d_layer_submodule.f90
@@ -0,0 +1,94 @@
+submodule(nf_avgpool2d_layer) nf_avgpool2d_layer_submodule
+  implicit none
+
+contains
+
+  pure module function avgpool2d_layer_cons(pool_size, stride) result(res)
+    implicit none
+    integer, intent(in) :: pool_size
+    integer, intent(in) :: stride
+    type(avgpool2d_layer) :: res
+
+    res % pool_size = pool_size
+    res % stride    = stride
+  end function avgpool2d_layer_cons
+
+
+  module subroutine init(self, input_shape)
+    implicit none
+    class(avgpool2d_layer), intent(in out) :: self
+    integer, intent(in) :: input_shape(:)
+    ! input_shape is expected to be (channels, width, height)
+
+    self % channels = input_shape(1)
+    self % width    = input_shape(2) / self % stride
+    self % height   = input_shape(3) / self % stride
+
+    ! Allocate the gradient array corresponding to the input dimensions.
+    allocate(self % gradient(input_shape(1), input_shape(2), input_shape(3)))
+    self % gradient = 0
+
+    ! Allocate the output array (after pooling).
+    allocate(self % output(self % channels, self % width, self % height))
+    self % output = 0
+  end subroutine init
+
+
+  pure module subroutine forward(self, input)
+    implicit none
+    class(avgpool2d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:,:)
+    integer :: input_width, input_height
+    integer :: i, j, n
+    integer :: ii, jj, iend, jend
+    integer :: iextent, jextent
+
+    input_width  = size(input, dim=2)
+    input_height = size(input, dim=3)
+    
+    ! Ensure we only process complete pooling regions.
+    iextent = input_width - mod(input_width, self % stride)
+    jextent = input_height - mod(input_height, self % stride)
+
+    ! Loop over the input with a step size equal to the stride and over all channels.
+    do concurrent (i = 1:iextent:self % stride, j = 1:jextent:self % stride, n = 1:self % channels)
+      ii = (i - 1) / self % stride + 1
+      jj = (j - 1) / self % stride + 1
+      
+      iend = min(i + self % pool_size - 1, input_width)
+      jend = min(j + self % pool_size - 1, input_height)
+      
+      ! Compute the average over the pooling region.
+      self % output(n, ii, jj) = sum(input(n, i:iend, j:jend)) / ((iend - i + 1) * (jend - j + 1))
+    end do
+  end subroutine forward
+
+
+  pure module subroutine backward(self, input, gradient)
+    implicit none
+    class(avgpool2d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:,:)
+    real, intent(in) :: gradient(:,:,:)
+    integer :: channels, pooled_width, pooled_height
+    integer :: i, j, n, x, y, istart, iend, jstart, jend
+    real :: scale_factor
+
+    channels      = size(gradient, dim=1)
+    pooled_width  = size(gradient, dim=2)
+    pooled_height = size(gradient, dim=3)
+
+    ! The gradient for average pooling is distributed evenly over the pooling window.
+    do concurrent (n = 1:channels, i = 1:pooled_width, j = 1:pooled_height)
+      istart = (i - 1) * self % stride + 1
+      iend   = min(istart + self % pool_size - 1, size(input, dim=2))
+      jstart = (j - 1) * self % stride + 1
+      jend   = min(jstart + self % pool_size - 1, size(input, dim=3))
+      scale_factor = 1.0 / ((iend - istart + 1) * (jend - jstart + 1))
+
+      do concurrent (x = istart:iend, y = jstart:jend)
+        self % gradient(n, x, y) = gradient(n, i, j) * scale_factor
+      end do
+    end do
+  end subroutine backward
+
+end submodule nf_avgpool2d_layer_submodule
diff --git a/src/nf/nf_avgpool3d_layer.f90 b/src/nf/nf_avgpool3d_layer.f90
new file mode 100644
index 00000000..580547e1
--- /dev/null
+++ b/src/nf/nf_avgpool3d_layer.f90
@@ -0,0 +1,67 @@
+module nf_avgpool3d_layer
+    !! This module provides the 3-d average pooling layer.
+  
+    use nf_base_layer, only: base_layer
+    implicit none
+  
+    private
+    public :: avgpool3d_layer
+  
+    type, extends(base_layer) :: avgpool3d_layer
+        integer :: channels
+        integer :: depth       ! Depth of the input
+        integer :: height      ! Height of the input
+        integer :: width       ! Width of the input
+        integer :: pool_size(3)  ! Pooling window size (depth, height, width)
+        integer :: stride(3)     ! Stride (depth, height, width)
+  
+        ! Gradient for the input (same shape as the input).
+        real, allocatable :: gradient(:,:,:,:)
+        ! Output after pooling (dimensions: (channels, new_depth, new_height, new_width)).
+        real, allocatable :: output(:,:,:,:)
+    contains
+        procedure :: init
+        procedure :: forward
+        procedure :: backward
+    end type avgpool3d_layer
+  
+    interface avgpool3d_layer
+        pure module function avgpool3d_layer_cons(pool_size, stride) result(res)
+            !! `avgpool3d` constructor function.
+            integer, intent(in) :: pool_size(3)
+                !! Depth, height, and width of the pooling window.
+            integer, intent(in) :: stride(3)
+                !! Stride of the pooling window (depth, height, width).
+            type(avgpool3d_layer) :: res
+        end function avgpool3d_layer_cons
+    end interface avgpool3d_layer
+  
+    interface
+        module subroutine init(self, input_shape)
+            !! Initialize the `avgpool3d` layer instance with an input shape.
+            class(avgpool3d_layer), intent(in out) :: self
+                !! `avgpool3d_layer` instance.
+            integer, intent(in) :: input_shape(:)
+                !! Array shape of the input layer, expected as (channels, depth, height, width).
+        end subroutine init
+  
+        pure module subroutine forward(self, input)
+            !! Run a forward pass of the `avgpool3d` layer.
+            class(avgpool3d_layer), intent(in out) :: self
+                !! `avgpool3d_layer` instance.
+            real, intent(in) :: input(:,:,:,:)
+                !! Input data (output of the previous layer), with shape (channels, depth, height, width).
+        end subroutine forward
+  
+        pure module subroutine backward(self, input, gradient)
+            !! Run a backward pass of the `avgpool3d` layer.
+            class(avgpool3d_layer), intent(in out) :: self
+                !! `avgpool3d_layer` instance.
+            real, intent(in) :: input(:,:,:,:)
+                !! Input data (output of the previous layer).
+            real, intent(in) :: gradient(:,:,:,:)
+                !! Gradient from the downstream layer, with shape (channels, pooled_depth, pooled_height, pooled_width).
+        end subroutine backward
+    end interface
+  
+end module nf_avgpool3d_layer
diff --git a/src/nf/nf_avgpool3d_layer_submodule.f90 b/src/nf/nf_avgpool3d_layer_submodule.f90
new file mode 100644
index 00000000..f2bd09dd
--- /dev/null
+++ b/src/nf/nf_avgpool3d_layer_submodule.f90
@@ -0,0 +1,102 @@
+submodule(nf_avgpool3d_layer) nf_avgpool3d_layer_submodule
+  implicit none
+
+contains
+
+  pure module function avgpool3d_layer_cons(pool_size, stride) result(res)
+    implicit none
+    integer, intent(in) :: pool_size(3)
+    integer, intent(in) :: stride(3)
+    type(avgpool3d_layer) :: res
+
+    res % pool_size = pool_size
+    res % stride    = stride
+  end function avgpool3d_layer_cons
+
+
+  module subroutine init(self, input_shape)
+    implicit none
+    class(avgpool3d_layer), intent(in out) :: self
+    integer, intent(in) :: input_shape(:)
+    ! input_shape is expected to be (channels, depth, height, width)
+
+    self % channels = input_shape(1)
+    self % depth    = input_shape(2) / self % stride(1)
+    self % height   = input_shape(3) / self % stride(2)
+    self % width    = input_shape(4) / self % stride(3)
+
+    ! Allocate the gradient array corresponding to the input dimensions.
+    allocate(self % gradient(input_shape(1), input_shape(2), input_shape(3), input_shape(4)))
+    self % gradient = 0
+
+    ! Allocate the output array (after pooling).
+    allocate(self % output(self % channels, self % depth, self % height, self % width))
+    self % output = 0
+  end subroutine init
+
+
+  pure module subroutine forward(self, input)
+    implicit none
+    class(avgpool3d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:,:,:)
+    integer :: input_depth, input_height, input_width
+    integer :: i, j, k, n
+    integer :: ii, jj, kk, iend, jend, kend
+    integer :: kdepth, jheight, kwidth
+
+    input_depth   = size(input, dim=2)
+    input_height  = size(input, dim=3)
+    input_width   = size(input, dim=4)
+    
+    ! Ensure we only process complete pooling regions.
+    kdepth = input_depth - mod(input_depth, self % stride(1))
+    jheight = input_height - mod(input_height, self % stride(2))
+    kwidth = input_width - mod(input_width, self % stride(3))
+
+    ! Loop over the input with a step size equal to the stride and over all channels.
+    do concurrent (i = 1:kdepth:self % stride(1), j = 1:jheight:self % stride(2), k = 1:kwidth:self % stride(3), n = 1:self % channels)
+      ii = (i - 1) / self % stride(1) + 1
+      jj = (j - 1) / self % stride(2) + 1
+      kk = (k - 1) / self % stride(3) + 1
+      
+      iend = min(i + self % pool_size(1) - 1, input_depth)
+      jend = min(j + self % pool_size(2) - 1, input_height)
+      kend = min(k + self % pool_size(3) - 1, input_width)
+      
+      ! Compute the average over the pooling region.
+      self % output(n, ii, jj, kk) = sum(input(n, i:iend, j:jend, k:kend)) / ((iend - i + 1) * (jend - j + 1) * (kend - k + 1))
+    end do
+  end subroutine forward
+
+
+  pure module subroutine backward(self, input, gradient)
+    implicit none
+    class(avgpool3d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:,:,:)
+    real, intent(in) :: gradient(:,:,:,:)
+    integer :: channels, pooled_depth, pooled_height, pooled_width
+    integer :: i, j, k, n, x, y, z, istart, iend, jstart, jend, kstart, kend
+    real :: scale_factor
+
+    channels      = size(gradient, dim=1)
+    pooled_depth  = size(gradient, dim=2)
+    pooled_height = size(gradient, dim=3)
+    pooled_width  = size(gradient, dim=4)
+
+    ! The gradient for average pooling is distributed evenly over the pooling window.
+    do concurrent (n = 1:channels, i = 1:pooled_depth, j = 1:pooled_height, k = 1:pooled_width)
+      istart = (i - 1) * self % stride(1) + 1
+      iend   = min(istart + self % pool_size(1) - 1, size(input, dim=2))
+      jstart = (j - 1) * self % stride(2) + 1
+      jend   = min(jstart + self % pool_size(2) - 1, size(input, dim=3))
+      kstart = (k - 1) * self % stride(3) + 1
+      kend   = min(kstart + self % pool_size(3) - 1, size(input, dim=4))
+      scale_factor = 1.0 / ((iend - istart + 1) * (jend - jstart + 1) * (kend - kstart + 1))
+
+      do concurrent (x = istart:iend, y = jstart:jend, z = kstart:kend)
+        self % gradient(n, x, y, z) = gradient(n, i, j, k) * scale_factor
+      end do
+    end do
+  end subroutine backward
+
+end submodule nf_avgpool3d_layer_submodule
diff --git a/src/nf/nf_layer_constructors.f90 b/src/nf/nf_layer_constructors.f90
index d3f06ca3..7c4c42b0 100644
--- a/src/nf/nf_layer_constructors.f90
+++ b/src/nf/nf_layer_constructors.f90
@@ -9,6 +9,9 @@ module nf_layer_constructors
 
   private
   public :: &
+    avgpool1d, &
+    avgpool2d, &
+    avgpool3d, &
     conv1d, &
     conv2d, &
     dense, &
@@ -19,6 +22,7 @@ module nf_layer_constructors
     locally_connected1d, &
     maxpool1d, &
     maxpool2d, &
+    maxpool3d, &
     reshape, &
     self_attention, &
     embedding, &
@@ -179,6 +183,72 @@ module function flatten() result(res)
         !! Resulting layer instance
     end function flatten
 
+    module function avgpool1d(pool_size, stride) result(res)
+      !! 1-d avgpooling layer constructor.
+      !!
+      !! This layer is for downscaling other layers, typically `conv1d`.
+      !!
+      !! Example:
+      !!
+      !! ```
+      !! use nf, only :: avgpool1d, layer
+      !! type(layer) :: avgpool1d_layer
+      !! avgpool1d_layer = avgpool1d(pool_size=2)
+      !! avgpool1d_layer = avgpool1d(pool_size=2, stride=3)
+      !! ```
+      integer, intent(in) :: pool_size
+        !! Width of the pooling window, commonly 2
+      integer, intent(in), optional :: stride
+        !! Stride of the pooling window, commonly equal to `pool_size`;
+        !! Defaults to `pool_size` if omitted.
+      type(layer) :: res
+        !! Resulting layer instance
+    end function avgpool1d
+
+    module function avgpool2d(pool_size, stride) result(res)
+      !! 2-d avgpooling layer constructor.
+      !!
+      !! This layer is for downscaling other layers, typically `conv2d`.
+      !!
+      !! Example:
+      !!
+      !! ```
+      !! use nf, only :: avgpool2d, layer
+      !! type(layer) :: avgpool2d_layer
+      !! avgpool2d_layer = avgpool2d(pool_size=2)
+      !! avgpool2d_layer = avgpool2d(pool_size=2, stride=3)
+      !! ```
+      integer, intent(in) :: pool_size
+        !! Width of the pooling window, commonly 2
+      integer, intent(in), optional :: stride
+        !! Stride of the pooling window, commonly equal to `pool_size`;
+        !! Defaults to `pool_size` if omitted.
+      type(layer) :: res
+        !! Resulting layer instance
+    end function avgpool2d
+
+    module function avgpool3d(pool_size, stride) result(res)
+      !! 3-d avgpooling layer constructor.
+      !!
+      !! This layer is for downscaling other layers, typically `conv3d`.
+      !!
+      !! Example:
+      !!
+      !! ```
+      !! use nf, only :: avgpool3d, layer
+      !! type(layer) :: avgpool3d_layer
+      !! avgpool3d_layer = avgpool3d(pool_size=2)
+      !! avgpool3d_layer = avgpool3d(pool_size=2, stride=3)
+      !! ```
+      integer, intent(in) :: pool_size
+        !! Width of the pooling window, commonly 2
+      integer, intent(in), optional :: stride
+        !! Stride of the pooling window, commonly equal to `pool_size`;
+        !! Defaults to `pool_size` if omitted.
+      type(layer) :: res
+        !! Resulting layer instance
+    end function avgpool3d
+
     module function conv1d(filters, kernel_size, activation) result(res)
       !! 1-d convolutional layer constructor.
       !!
@@ -304,6 +374,28 @@ module function maxpool2d(pool_size, stride) result(res)
         !! Resulting layer instance
     end function maxpool2d
 
+    module function maxpool3d(pool_size, stride) result(res)
+      !! 3-d maxpooling layer constructor.
+      !!
+      !! This layer is for downscaling other layers, typically `conv3d`.
+      !!
+      !! Example:
+      !!
+      !! ```
+      !! use nf, only :: maxpool3d, layer
+      !! type(layer) :: maxpool3d_layer
+      !! maxpool3d_layer = maxpool3d(pool_size=2)
+      !! maxpool3d_layer = maxpool3d(pool_size=2, stride=3)
+      !! ```
+      integer, intent(in) :: pool_size
+        !! Width of the pooling window, commonly 2
+      integer, intent(in), optional :: stride
+        !! Stride of the pooling window, commonly equal to `pool_size`;
+        !! Defaults to `pool_size` if omitted.
+      type(layer) :: res
+        !! Resulting layer instance
+    end function maxpool3d
+
     module function linear2d(out_features) result(res)
       !! Rank-2 (sequence_length, out_features) linear layer constructor.
       !! sequence_length is determined at layer initialization, based on the
diff --git a/src/nf/nf_layer_constructors_submodule.f90 b/src/nf/nf_layer_constructors_submodule.f90
index 1665d38a..d6e05cc2 100644
--- a/src/nf/nf_layer_constructors_submodule.f90
+++ b/src/nf/nf_layer_constructors_submodule.f90
@@ -1,6 +1,9 @@
 submodule(nf_layer_constructors) nf_layer_constructors_submodule
 
   use nf_layer, only: layer
+  use nf_avgpool1d_layer, only: avgpool1d_layer
+  use nf_avgpool2d_layer, only: avgpool2d_layer
+  use nf_avgpool3d_layer, only: avgpool3d_layer
   use nf_conv1d_layer, only: conv1d_layer
   use nf_conv2d_layer, only: conv2d_layer
   use nf_dense_layer, only: dense_layer
@@ -12,6 +15,7 @@
   use nf_locally_connected1d_layer, only: locally_connected1d_layer
   use nf_maxpool1d_layer, only: maxpool1d_layer
   use nf_maxpool2d_layer, only: maxpool2d_layer
+  use nf_maxpool3d_layer, only: maxpool3d_layer
   use nf_reshape2d_layer, only: reshape2d_layer
   use nf_reshape3d_layer, only: reshape3d_layer
   use nf_linear2d_layer, only: linear2d_layer
@@ -140,6 +144,92 @@ module function flatten() result(res)
   end function flatten
 
 
+  module function avgpool1d(pool_size, stride) result(res)
+    integer, intent(in) :: pool_size
+    integer, intent(in), optional :: stride
+    integer :: stride_
+    type(layer) :: res
+
+    if (pool_size < 2) &
+      error stop 'pool_size must be >= 2 in a avgpool1d layer'
+
+    ! Stride defaults to pool_size if not provided
+    if (present(stride)) then
+      stride_ = stride
+    else
+      stride_ = pool_size
+    end if
+
+    if (stride_ < 1) &
+      error stop 'stride must be >= 1 in a avgpool1d layer'
+
+    res % name = 'avgpool1d'
+
+    allocate( &
+      res % p, &
+      source=avgpool1d_layer(pool_size, stride_) &
+    )
+
+  end function avgpool1d
+
+  module function avgpool2d(pool_size, stride) result(res)
+    integer, intent(in) :: pool_size
+    integer, intent(in), optional :: stride
+    integer :: stride_
+    type(layer) :: res
+
+    if (pool_size < 2) &
+      error stop 'pool_size must be >= 2 in a avgpool2d layer'
+
+    ! Stride defaults to pool_size if not provided
+    if (present(stride)) then
+      stride_ = stride
+    else
+      stride_ = pool_size
+    end if
+
+    if (stride_ < 1) &
+      error stop 'stride must be >= 1 in a avgpool2d layer'
+
+    res % name = 'avgpool2d'
+
+    allocate( &
+      res % p, &
+      source=avgpool2d_layer(pool_size, stride_) &
+    )
+
+  end function avgpool2d
+
+
+  module function avgpool3d(pool_size, stride) result(res)
+    integer, intent(in) :: pool_size
+    integer, intent(in), optional :: stride
+    integer :: stride_
+    type(layer) :: res
+
+    if (pool_size < 2) &
+      error stop 'pool_size must be >= 2 in a avgpool3d layer'
+
+    ! Stride defaults to pool_size if not provided
+    if (present(stride)) then
+      stride_ = stride
+    else
+      stride_ = pool_size
+    end if
+
+    if (stride_ < 1) &
+      error stop 'stride must be >= 1 in a avgpool3d layer'
+
+    res % name = 'avgpool3d'
+
+    allocate( &
+      res % p, &
+      source=avgpool3d_layer(pool_size, stride_) &
+    )
+
+  end function avgpool3d
+
+
   module function input1d(layer_size) result(res)
     integer, intent(in) :: layer_size
     type(layer) :: res
@@ -229,6 +319,35 @@ module function maxpool2d(pool_size, stride) result(res)
   end function maxpool2d
 
 
+  module function maxpool3d(pool_size, stride) result(res)
+    integer, intent(in) :: pool_size
+    integer, intent(in), optional :: stride
+    integer :: stride_
+    type(layer) :: res
+
+    if (pool_size < 2) &
+      error stop 'pool_size must be >= 2 in a maxpool3d layer'
+
+    ! Stride defaults to pool_size if not provided
+    if (present(stride)) then
+      stride_ = stride
+    else
+      stride_ = pool_size
+    end if
+
+    if (stride_ < 1) &
+      error stop 'stride must be >= 1 in a maxpool3d layer'
+
+    res % name = 'maxpool3d'
+
+    allocate( &
+      res % p, &
+      source=maxpool3d_layer(pool_size, stride_) &
+    )
+
+  end function maxpool3d
+
+
   module function reshape2d(dim1, dim2) result(res)
     integer, intent(in) :: dim1, dim2
     type(layer) :: res
diff --git a/src/nf/nf_layer_submodule.f90 b/src/nf/nf_layer_submodule.f90
index eebedaa9..0d3bd573 100644
--- a/src/nf/nf_layer_submodule.f90
+++ b/src/nf/nf_layer_submodule.f90
@@ -1,6 +1,9 @@
 submodule(nf_layer) nf_layer_submodule
 
   use iso_fortran_env, only: stderr => error_unit
+  use nf_avgpool1d_layer, only: avgpool1d_layer
+  use nf_avgpool2d_layer, only: avgpool2d_layer
+  use nf_avgpool3d_layer, only: avgpool3d_layer
   use nf_conv1d_layer, only: conv1d_layer
   use nf_conv2d_layer, only: conv2d_layer
   use nf_dense_layer, only: dense_layer
@@ -12,6 +15,7 @@
   use nf_locally_connected1d_layer, only: locally_connected1d_layer
   use nf_maxpool1d_layer, only: maxpool1d_layer
   use nf_maxpool2d_layer, only: maxpool2d_layer
+  use nf_maxpool3d_layer, only: maxpool3d_layer
   use nf_reshape2d_layer, only: reshape2d_layer
   use nf_reshape3d_layer, only: reshape3d_layer
   use nf_linear2d_layer, only: linear2d_layer
@@ -54,6 +58,14 @@ pure module subroutine backward_1d(self, previous, gradient)
 
         ! Upstream layers permitted: input2d, input3d, conv1d, conv2d, locally_connected1d, maxpool1d, maxpool2d
         select type(prev_layer => previous % p)
+          type is(maxpool3d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool1d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool3d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
           type is(input2d_layer)
             call this_layer % backward(prev_layer % output, gradient)
           type is(locally_connected1d_layer)
@@ -141,6 +153,8 @@ pure module subroutine backward_2d(self, previous, gradient)
       select type(prev_layer => previous % p)
         type is(maxpool1d_layer)
           call this_layer % backward(prev_layer % output, gradient)
+        type is(avgpool1d_layer)
+          call this_layer % backward(prev_layer % output, gradient)
         type is(reshape2d_layer)
           call this_layer % backward(prev_layer % output, gradient)
         type is(input2d_layer)
@@ -156,6 +170,8 @@ pure module subroutine backward_2d(self, previous, gradient)
         select type(prev_layer => previous % p)
           type is(maxpool1d_layer)
             call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool1d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
           type is(reshape2d_layer)
             call this_layer % backward(prev_layer % output, gradient)
           type is(input2d_layer)
@@ -171,6 +187,25 @@ pure module subroutine backward_2d(self, previous, gradient)
         select type(prev_layer => previous % p)
           type is(maxpool1d_layer)
             call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool1d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(reshape2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(locally_connected1d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(input2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(conv1d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+        end select
+
+      type is(avgpool1d_layer)
+
+        select type(prev_layer => previous % p)
+          type is(maxpool1d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool1d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
           type is(reshape2d_layer)
             call this_layer % backward(prev_layer % output, gradient)
           type is(locally_connected1d_layer)
@@ -208,6 +243,8 @@ pure module subroutine backward_3d(self, previous, gradient)
         select type(prev_layer => previous % p)
           type is(maxpool2d_layer)
             call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
           type is(input3d_layer)
             call this_layer % backward(prev_layer % output, gradient)
           type is(conv2d_layer)
@@ -224,6 +261,24 @@ pure module subroutine backward_3d(self, previous, gradient)
             call this_layer % backward(prev_layer % output, gradient)
           type is(maxpool2d_layer)
             call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(input3d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(reshape3d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+        end select
+
+      type is(avgpool2d_layer)
+
+        ! Upstream layers permitted: conv2d, input3d, maxpool2d, reshape3d
+        select type(prev_layer => previous % p)
+          type is(conv2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(maxpool2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
+          type is(avgpool2d_layer)
+            call this_layer % backward(prev_layer % output, gradient)
           type is(input3d_layer)
             call this_layer % backward(prev_layer % output, gradient)
           type is(reshape3d_layer)
@@ -356,6 +411,36 @@ module subroutine forward(self, input)
             call this_layer % forward(prev_layer % output)
         end select
 
+      type is(avgpool1d_layer)
+
+        ! Upstream layers permitted: input1d, locally_connected1d, maxpool1d, reshape2d
+        select type(prev_layer => input % p)
+          type is(input2d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(locally_connected1d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(maxpool1d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(reshape2d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(conv1d_layer)
+            call this_layer % forward(prev_layer % output)
+        end select
+
+      type is(avgpool2d_layer)
+
+        ! Upstream layers permitted: input3d, conv2d, maxpool2d, reshape3d
+        select type(prev_layer => input % p)
+          type is(input3d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(conv2d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(maxpool2d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(reshape3d_layer)
+            call this_layer % forward(prev_layer % output)
+        end select
+
       type is(flatten_layer)
 
         ! Upstream layers permitted: input2d, input3d, conv2d, maxpool1d, maxpool2d, reshape2d, reshape3d, locally_connected2d
@@ -374,6 +459,14 @@ module subroutine forward(self, input)
             call this_layer % forward(prev_layer % output)
           type is(maxpool2d_layer)
             call this_layer % forward(prev_layer % output)
+          type is(maxpool3d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(avgpool1d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(avgpool2d_layer)
+            call this_layer % forward(prev_layer % output)
+          type is(avgpool3d_layer)
+            call this_layer % forward(prev_layer % output)
           type is(reshape2d_layer)
             call this_layer % forward(prev_layer % output)
           type is(reshape3d_layer)
@@ -481,6 +574,8 @@ pure module subroutine get_output_2d(self, output)
         allocate(output, source=this_layer % output)
       type is(maxpool1d_layer)
         allocate(output, source=this_layer % output)
+      type is(avgpool1d_layer)
+        allocate(output, source=this_layer % output)
       type is(locally_connected1d_layer)
         allocate(output, source=this_layer % output)
       type is(conv1d_layer)
@@ -518,6 +613,8 @@ pure module subroutine get_output_3d(self, output)
         allocate(output, source=this_layer % output)
       type is(maxpool2d_layer)
         allocate(output, source=this_layer % output)
+      type is(avgpool2d_layer)
+        allocate(output, source=this_layer % output)
       type is(reshape3d_layer)
         allocate(output, source=this_layer % output)
       class default
@@ -563,6 +660,14 @@ impure elemental module subroutine init(self, input)
         self % layer_shape = shape(this_layer % output)
       type is(maxpool2d_layer)
         self % layer_shape = shape(this_layer % output)
+      type is(maxpool3d_layer)
+        self % layer_shape = shape(this_layer % output)
+      type is(avgpool1d_layer)
+        self % layer_shape = shape(this_layer % output)
+      type is(avgpool2d_layer)
+        self % layer_shape = shape(this_layer % output)
+      type is(avgpool3d_layer)
+        self % layer_shape = shape(this_layer % output)
     end select
 
     self % input_layer_shape = input % layer_shape
@@ -617,6 +722,14 @@ elemental module function get_num_params(self) result(num_params)
         num_params = 0
       type is (maxpool2d_layer)
         num_params = 0
+      type is (maxpool3d_layer)
+        num_params = 0
+      type is (avgpool1d_layer)
+        num_params = 0
+      type is (avgpool2d_layer)
+        num_params = 0
+      type is (avgpool3d_layer)
+        num_params = 0
       type is (flatten_layer)
         num_params = 0
       type is (reshape2d_layer)
@@ -662,6 +775,14 @@ module function get_params(self) result(params)
         ! No parameters to get.
       type is (maxpool2d_layer)
         ! No parameters to get.
+      type is (maxpool3d_layer)
+      ! No parameters to get.
+      type is (avgpool1d_layer)
+      ! No parameters to get.
+      type is (avgpool2d_layer)
+      ! No parameters to get.
+      type is (avgpool3d_layer)
+      ! No parameters to get.
       type is (flatten_layer)
         ! No parameters to get.
       type is (reshape2d_layer)
@@ -707,6 +828,14 @@ module function get_gradients(self) result(gradients)
         ! No gradients to get.
       type is (maxpool2d_layer)
         ! No gradients to get.
+      type is (maxpool3d_layer)
+      ! No gradients to get.
+      type is (avgpool1d_layer)
+      ! No gradients to get.
+      type is (avgpool2d_layer)
+      ! No gradients to get.
+      type is (avgpool3d_layer)
+      ! No gradients to get.
       type is (flatten_layer)
         ! No gradients to get.
       type is (reshape2d_layer)
@@ -783,6 +912,26 @@ module subroutine set_params(self, params)
         ! No parameters to set.
         write(stderr, '(a)') 'Warning: calling set_params() ' &
           // 'on a zero-parameter layer; nothing to do.'
+        
+      type is (maxpool3d_layer)
+        ! No parameters to set.
+        write(stderr, '(a)') 'Warning: calling set_params() ' &
+          // 'on a zero-parameter layer; nothing to do.'
+
+      type is (avgpool1d_layer)
+          ! No parameters to set.
+          write(stderr, '(a)') 'Warning: calling set_params() ' &
+            // 'on a zero-parameter layer; nothing to do.'
+
+      type is (avgpool2d_layer)
+        ! No parameters to set.
+        write(stderr, '(a)') 'Warning: calling set_params() ' &
+          // 'on a zero-parameter layer; nothing to do.'
+        
+      type is (avgpool3d_layer)
+        ! No parameters to set.
+        write(stderr, '(a)') 'Warning: calling set_params() ' &
+          // 'on a zero-parameter layer; nothing to do.'
 
       type is (linear2d_layer)
         call this_layer % set_params(params)
diff --git a/src/nf/nf_maxpool3d_layer.f90 b/src/nf/nf_maxpool3d_layer.f90
new file mode 100644
index 00000000..b4815efe
--- /dev/null
+++ b/src/nf/nf_maxpool3d_layer.f90
@@ -0,0 +1,79 @@
+module nf_maxpool3d_layer
+
+    !! This module provides the 3D maxpooling layer.
+  
+    use nf_base_layer, only: base_layer
+    implicit none
+  
+    private
+    public :: maxpool3d_layer
+  
+    type, extends(base_layer) :: maxpool3d_layer
+  
+      integer :: channels
+      integer :: depth
+      integer :: width
+      integer :: height
+      integer :: pool_size
+      integer :: stride
+  
+      ! Locations (as input matrix indices) of the maximum values
+      ! in the depth (z), width (x), and height (y) dimensions
+      integer, allocatable :: maxloc_x(:,:,:,:)
+      integer, allocatable :: maxloc_y(:,:,:,:)
+      integer, allocatable :: maxloc_z(:,:,:,:)
+  
+      real, allocatable :: gradient(:,:,:,:)
+      real, allocatable :: output(:,:,:,:)
+  
+    contains
+  
+      procedure :: init
+      procedure :: forward
+      procedure :: backward
+  
+    end type maxpool3d_layer
+  
+    interface maxpool3d_layer
+      pure module function maxpool3d_layer_cons(pool_size, stride) result(res)
+        !! `maxpool3d` constructor function
+        integer, intent(in) :: pool_size
+          !! Depth, width, and height of the pooling window
+        integer, intent(in) :: stride
+          !! Stride of the pooling window
+        type(maxpool3d_layer) :: res
+      end function maxpool3d_layer_cons
+    end interface maxpool3d_layer
+  
+    interface
+  
+      module subroutine init(self, input_shape)
+        !! Initialize the `maxpool3d` layer instance with an input shape.
+        class(maxpool3d_layer), intent(in out) :: self
+          !! `maxpool3d_layer` instance
+        integer, intent(in) :: input_shape(:)
+          !! Array shape of the input layer
+      end subroutine init
+  
+      pure module subroutine forward(self, input)
+        !! Run a forward pass of the `maxpool3d` layer.
+        class(maxpool3d_layer), intent(in out) :: self
+          !! `maxpool3d_layer` instance
+        real, intent(in) :: input(:,:,:,:)
+          !! Input data (output of the previous layer)
+      end subroutine forward
+  
+      pure module subroutine backward(self, input, gradient)
+        !! Run a backward pass of the `maxpool3d` layer.
+        class(maxpool3d_layer), intent(in out) :: self
+          !! `maxpool3d_layer` instance
+        real, intent(in) :: input(:,:,:,:)
+          !! Input data (output of the previous layer)
+        real, intent(in) :: gradient(:,:,:,:)
+          !! Gradient from the downstream layer
+      end subroutine backward
+  
+    end interface
+  
+  end module nf_maxpool3d_layer
+  
\ No newline at end of file
diff --git a/src/nf/nf_maxpool3d_layer_submodule.f90 b/src/nf/nf_maxpool3d_layer_submodule.f90
new file mode 100644
index 00000000..e5cf2bad
--- /dev/null
+++ b/src/nf/nf_maxpool3d_layer_submodule.f90
@@ -0,0 +1,114 @@
+submodule(nf_maxpool3d_layer) nf_maxpool3d_layer_submodule
+
+  implicit none
+
+contains
+
+  pure module function maxpool3d_layer_cons(pool_size, stride) result(res)
+    implicit none
+    integer, intent(in) :: pool_size
+    integer, intent(in) :: stride
+    type(maxpool3d_layer) :: res
+    res % pool_size = pool_size
+    res % stride = stride
+  end function maxpool3d_layer_cons
+
+
+  module subroutine init(self, input_shape)
+    implicit none
+    class(maxpool3d_layer), intent(in out) :: self
+    integer, intent(in) :: input_shape(:)
+
+    self % channels = input_shape(1)
+    self % depth = input_shape(2) / self % stride
+    self % width = input_shape(3) / self % stride
+    self % height = input_shape(4) / self % stride
+
+    allocate(self % maxloc_x(self % channels, self % depth, self % width, self % height))
+    allocate(self % maxloc_y(self % channels, self % depth, self % width, self % height))
+    allocate(self % maxloc_z(self % channels, self % depth, self % width, self % height))
+    self % maxloc_x = 0
+    self % maxloc_y = 0
+    self % maxloc_z = 0
+
+    allocate(self % gradient(input_shape(1), input_shape(2), input_shape(3), input_shape(4)))
+    self % gradient = 0
+
+    allocate(self % output(self % channels, self % depth, self % width, self % height))
+    self % output = 0
+
+  end subroutine init
+
+
+  pure module subroutine forward(self, input)
+    implicit none
+    class(maxpool3d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:,:,:)
+    integer :: input_depth, input_width, input_height
+    integer :: i, j, k, n
+    integer :: ii, jj, kk
+    integer :: iend, jend, kend
+    integer :: iextent, jextent, kextent
+    integer :: maxloc_xyz(3)
+
+    input_depth = size(input, dim=2)
+    input_width = size(input, dim=3)
+    input_height = size(input, dim=4)
+
+    kextent = input_depth - mod(input_depth, self % stride)
+    iextent = input_width - mod(input_width, self % stride)
+    jextent = input_height - mod(input_height, self % stride)
+
+    do concurrent( &
+      k = 1:kextent:self % stride, &
+      i = 1:iextent:self % stride, &
+      j = 1:jextent:self % stride &
+    )
+
+      kk = k / self % stride + 1
+      ii = i / self % stride + 1
+      jj = j / self % stride + 1
+
+      kend = k + self % pool_size - 1
+      iend = i + self % pool_size - 1
+      jend = j + self % pool_size - 1
+
+      do concurrent(n = 1:self % channels)
+        maxloc_xyz = maxloc(input(n, k:kend, i:iend, j:jend))
+        self % maxloc_x(n,kk,ii,jj) = maxloc_xyz(1) + k - 1
+        self % maxloc_y(n,kk,ii,jj) = maxloc_xyz(2) + i - 1
+        self % maxloc_z(n,kk,ii,jj) = maxloc_xyz(3) + j - 1
+
+        self % output(n,kk,ii,jj) = &
+          input(n,self % maxloc_x(n,kk,ii,jj),self % maxloc_y(n,kk,ii,jj),self % maxloc_z(n,kk,ii,jj))
+      end do
+
+    end do
+
+  end subroutine forward
+
+
+  pure module subroutine backward(self, input, gradient)
+    implicit none
+    class(maxpool3d_layer), intent(in out) :: self
+    real, intent(in) :: input(:,:,:,:)
+    real, intent(in) :: gradient(:,:,:,:)
+    integer :: gradient_shape(4)
+    integer :: channels, depth, width, height
+    integer :: i, j, k, n
+
+    gradient_shape = shape(gradient)
+    channels = gradient_shape(1)
+    depth = gradient_shape(2)
+    width = gradient_shape(3)
+    height = gradient_shape(4)
+
+    do concurrent(n = 1:channels, k = 1:depth, i = 1:width, j = 1:height)
+      associate(ii => self % maxloc_x(n,k,i,j), jj => self % maxloc_y(n,k,i,j), kk => self % maxloc_z(n,k,i,j))
+        self % gradient(n,ii,jj,kk) = gradient(n,k,i,j)
+      end associate
+    end do
+
+  end subroutine backward
+
+end submodule nf_maxpool3d_layer_submodule
\ No newline at end of file