diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md
index f8cf15227..cfd632456 100644
--- a/doc/specs/stdlib_math.md
+++ b/doc/specs/stdlib_math.md
@@ -61,6 +61,41 @@ Here inputs are of type `real` and kind `sp`
 {!example/math/example_clip_real.f90!}
 ```
 
+<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
+### `swap` subroutine
+
+#### Description
+
+Swaps the values in `lhs` and `rhs`.
+
+#### Syntax
+
+`call` [[stdlib_math(module):swap(interface)]] ` (lhs, rhs)`
+
+#### Status
+
+Experimental
+
+#### Class
+
+Elemental subroutine.
+
+#### Argument(s)
+
+`lhs`: scalar or array of any of the intrinsic types `integer`, `real`, `complex`, `logical`, `character`, `string_type`, `bitset` type. This argument is `intent(inout)`.
+`rhs`: scalar or array of any of the intrinsic types `integer`, `real`, `complex`, `logical`, `character`, `string_type`, `bitset` type. This argument is `intent(inout)`.
+
+##### Note
+All arguments must have same `type` and same `kind`.
+
+**WARNING**: For fix size characters with different length, the `swap` subroutine will truncate the longest amongst `lhs` and `rhs`. To avoid truncation it is possible to pass a subsection of the string.
+
+#### Examples
+
+```fortran
+{!example/math/example_math_swap.f90!}
+```
+
 ### `gcd` function
 
 #### Description
diff --git a/example/math/CMakeLists.txt b/example/math/CMakeLists.txt
index 0a1411bb9..8b834f348 100644
--- a/example/math/CMakeLists.txt
+++ b/example/math/CMakeLists.txt
@@ -15,4 +15,5 @@ ADD_EXAMPLE(math_argpi)
 ADD_EXAMPLE(math_deg2rad)
 ADD_EXAMPLE(math_rad2deg)
 ADD_EXAMPLE(math_is_close)
+ADD_EXAMPLE(math_swap)
 ADD_EXAMPLE(meshgrid)
diff --git a/example/math/example_math_swap.f90 b/example/math/example_math_swap.f90
new file mode 100644
index 000000000..7c2388b05
--- /dev/null
+++ b/example/math/example_math_swap.f90
@@ -0,0 +1,54 @@
+program example_math_swap
+    use stdlib_math, only: swap
+    implicit none
+    
+    block
+        integer :: x, y
+        x = 9
+        y = 18
+        call swap(x,y)
+    end block
+
+    block
+        real :: x, y
+        x = 4.0
+        y = 8.0
+        call swap(x,y)
+    end block
+
+    block
+        real :: x(3), y(3)
+        x = [1.0,2.0,3.0]
+        y = [4.0,5.0,6.0]
+        call swap(x,y)
+    end block
+
+    block
+        character(4) :: x
+        character(6) :: y
+        x = 'abcd'
+        y = 'efghij'
+        call swap(x,y)      ! x=efgh,  y=abcd
+
+        x = 'abcd'
+        y = 'efghij'
+        call swap(x,y(1:4)) ! x=efgh,  y=abcdij
+    end block
+
+    block
+        use stdlib_string_type
+        type(string_type) :: x, y
+        x = 'abcde'
+        y = 'fghij'
+        call swap(x,y)
+    end block
+
+    block
+        use stdlib_bitsets
+        type(bitset_64) :: x, y
+        call x%from_string('0000')
+        call y%from_string('1111')
+        call swap(x,y)
+    end block
+
+end program example_math_swap  
\ No newline at end of file
diff --git a/src/stdlib_math.fypp b/src/stdlib_math.fypp
index 92a946cd9..1b1abb363 100644
--- a/src/stdlib_math.fypp
+++ b/src/stdlib_math.fypp
@@ -1,14 +1,15 @@
 #:include "common.fypp"
 #:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES
 #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
-
+#:set BITSET_KINDS_TYPES = list(zip(BITSET_KINDS, BITSET_TYPES))
 module stdlib_math
     use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp
     use stdlib_optval, only: optval
+    use stdlib_bitsets, only: bitset_64, bitset_large
 
     implicit none
     private
-    public :: clip, gcd, linspace, logspace
+    public :: clip, swap, gcd, linspace, logspace
     public :: EULERS_NUMBER_SP, EULERS_NUMBER_DP
 #:if WITH_QP
     public :: EULERS_NUMBER_QP
@@ -42,6 +43,22 @@ module stdlib_math
     #:endfor
     end interface clip
 
+    !> Swap the values of the lhs and rhs arguments
+    !> ([Specification](../page/specs/stdlib_math.html#swap_subroutine))
+    !>
+    !> Version: experimental
+    interface swap
+      #:for k1, t1 in INT_KINDS_TYPES + REAL_KINDS_TYPES + BITSET_KINDS_TYPES
+      module procedure :: swap_${k1}$
+      #:endfor
+      #:for k1, t1 in CMPLX_KINDS_TYPES
+      module procedure :: swap_c${k1}$
+      #:endfor
+      module procedure :: swap_bool
+      module procedure :: swap_str
+      module procedure :: swap_stt
+    end interface
+
     !> Returns the greatest common divisor of two integers
     !> ([Specification](../page/specs/stdlib_math.html#gcd))
     !>
@@ -509,5 +526,42 @@ contains
     end function gcd_${k1}$
 
     #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES + REAL_KINDS_TYPES + BITSET_KINDS_TYPES
+    elemental subroutine swap_${k1}$(lhs, rhs)
+        ${t1}$, intent(inout) :: lhs, rhs
+        ${t1}$ :: temp
+        temp = lhs; lhs = rhs; rhs = temp
+    end subroutine
+
+    #:endfor
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    elemental subroutine swap_c${k1}$(lhs, rhs)
+        ${t1}$, intent(inout) :: lhs, rhs
+        ${t1}$ :: temp
+        temp = lhs; lhs = rhs; rhs = temp
+    end subroutine
+
+    #:endfor
+
+    elemental subroutine swap_bool(lhs, rhs)
+        logical, intent(inout) :: lhs, rhs
+        logical :: temp
+        temp = lhs; lhs = rhs; rhs = temp
+    end subroutine
+
+    elemental subroutine swap_str(lhs,rhs)
+        character(*), intent(inout) :: lhs, rhs
+        character(len=max(len(lhs), len(rhs))) :: temp
+        temp = lhs ; lhs = rhs ; rhs = temp
+    end subroutine
+
+    elemental subroutine swap_stt(lhs,rhs)
+        use stdlib_string_type, only: string_type
+        type(string_type), intent(inout) :: lhs, rhs
+        type(string_type) :: temp
+        temp = lhs ; lhs = rhs ; rhs = temp
+    end subroutine
     
 end module stdlib_math
diff --git a/test/math/test_stdlib_math.fypp b/test/math/test_stdlib_math.fypp
index 1f64622e4..1365756b9 100644
--- a/test/math/test_stdlib_math.fypp
+++ b/test/math/test_stdlib_math.fypp
@@ -4,7 +4,7 @@
 
 module test_stdlib_math
     use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
-    use stdlib_math, only: clip, arg, argd, argpi, arange, is_close, all_close, diff, &
+    use stdlib_math, only: clip, swap, arg, argd, argpi, arange, is_close, all_close, diff, &
                            arange, deg2rad, rad2deg
     use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp
     implicit none
@@ -38,6 +38,16 @@ contains
             new_unittest("clip-real-quad", test_clip_rqp), &
             new_unittest("clip-real-quad-bounds", test_clip_rqp_bounds) &
             
+            !> Tests swap
+            #:for k1, t1 in INT_KINDS_TYPES + REAL_KINDS_TYPES
+            , new_unittest("swap_${k1}$", test_swap_${k1}$) &
+            #:endfor
+            #:for k1, t1 in CMPLX_KINDS_TYPES
+            , new_unittest("swap_c${k1}$", test_swap_c${k1}$) &
+            #:endfor
+            , new_unittest("swap_str", test_swap_str) &
+            , new_unittest("swap_stt", test_swap_stt) &
+
             !> Tests for arg/argd/argpi
             #:for k1 in CMPLX_KINDS
             , new_unittest("arg-cmplx-${k1}$", test_arg_${k1}$) &
@@ -246,6 +256,166 @@ contains
 
     end subroutine test_clip_rqp_bounds
 
+    #:for k1, t1 in INT_KINDS_TYPES + REAL_KINDS_TYPES
+    subroutine test_swap_${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        ${t1}$ :: x(3), y(3)
+        
+        x = [${t1}$ :: 1, 2, 3]
+        y = [${t1}$ :: 4, 5, 6]
+
+        call swap(x,y)
+        
+        call check(error, all( x == [${t1}$ :: 4, 5, 6] ) )
+        if (allocated(error)) return
+        call check(error, all( y == [${t1}$ :: 1, 2, 3] ) )
+        if (allocated(error)) return
+
+        ! check self swap
+        call swap(x,x)
+
+        call check(error, all( x == [${t1}$ :: 4, 5, 6] ) )
+        if (allocated(error)) return
+    end subroutine test_swap_${k1}$
+    #:endfor
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    subroutine test_swap_c${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        ${t1}$ :: x(3), y(3)
+        
+        x = cmplx( [1, 2, 3] , [4, 5, 6] )
+        y = cmplx( [4, 5, 6] , [1, 2, 3] )
+
+        call swap(x,y)
+        
+        call check(error, all( x == cmplx( [4, 5, 6] , [1, 2, 3] ) ) )
+        if (allocated(error)) return
+        call check(error, all( y == cmplx( [1, 2, 3] , [4, 5, 6] ) ) )
+        if (allocated(error)) return
+
+        ! check self swap
+        call swap(x,x)
+        
+        call check(error, all( x == cmplx( [4, 5, 6] , [1, 2, 3] ) ) )
+        if (allocated(error)) return
+    end subroutine test_swap_c${k1}$
+    #:endfor
+
+    subroutine test_swap_str(error)
+        type(error_type), allocatable, intent(out) :: error
+        block
+            character(5) :: x(2), y(2)
+            
+            x = ['abcde','fghij'] 
+            y = ['fghij','abcde']
+
+            call swap(x,y)
+            
+            call check(error, all( x == ['fghij','abcde'] ) )
+            if (allocated(error)) return
+            call check(error, all( y == ['abcde','fghij']  ) )
+            if (allocated(error)) return
+
+            ! check self swap
+            call swap(x,x)
+            
+            call check(error, all( x == ['fghij','abcde'] ) )
+            if (allocated(error)) return
+        end block
+
+        block
+            character(4) :: x
+            character(6) :: y
+            
+            x = 'abcd'
+            y = 'efghij'
+            call swap(x,y)
+            
+            call check(error, x == 'efgh' )
+            if (allocated(error)) return
+            call check(error, y(1:6) == 'abcd  ' )
+            if (allocated(error)) return
+
+            x = 'abcd'
+            y = 'efghij'
+            call swap(x,y(1:4))
+            
+            call check(error, x == 'efgh' )
+            if (allocated(error)) return
+            call check(error, y == 'abcdij' )
+            if (allocated(error)) return
+        end block
+    end subroutine test_swap_str
+
+    subroutine test_swap_stt(error)
+        use stdlib_string_type
+        type(error_type), allocatable, intent(out) :: error
+        type(string_type) :: x(2), y(2)
+        
+        x = ['abcde','fghij'] 
+        y = ['fghij','abcde']
+
+        call swap(x,y)
+        
+        call check(error, all( x == ['fghij','abcde'] ) )
+        if (allocated(error)) return
+        call check(error, all( y == ['abcde','fghij']  ) )
+        if (allocated(error)) return
+
+        ! check self swap
+        call swap(x,x)
+        
+        call check(error, all( x == ['fghij','abcde'] ) )
+        if (allocated(error)) return
+    end subroutine test_swap_stt
+
+    subroutine test_swap_bitset_64(error)
+        use stdlib_bitsets
+        type(error_type), allocatable, intent(out) :: error
+        type(bitset_64) :: x, y, u, v
+        
+        x = [.true.,.false.,.true.,.false.] 
+        u = x
+        y = [.false.,.true.,.false.,.true.]
+        v = y
+        call swap(x,y)
+        
+        call check(error, x == v )
+        if (allocated(error)) return
+        call check(error, y == u  )
+        if (allocated(error)) return
+
+        ! check self swap
+        call swap(x,x)
+        
+        call check(error, x == v )
+        if (allocated(error)) return
+    end subroutine test_swap_bitset_64
+
+    subroutine test_swap_bitset_large(error)
+        use stdlib_bitsets
+        type(error_type), allocatable, intent(out) :: error
+        type(bitset_large) :: x, y, u, v
+        
+        x = [.true.,.false.,.true.,.false.] 
+        u = x
+        y = [.false.,.true.,.false.,.true.]
+        v = y
+        call swap(x,y)
+        
+        call check(error, x == v )
+        if (allocated(error)) return
+        call check(error, y == u )
+        if (allocated(error)) return
+
+        ! check self swap
+        call swap(x,x)
+        
+        call check(error, x == v )
+        if (allocated(error)) return
+    end subroutine test_swap_bitset_large
+
     #:for k1 in CMPLX_KINDS
     subroutine test_arg_${k1}$(error)
         type(error_type), allocatable, intent(out) :: error