diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md
index 6aa76fa88..cab16279c 100644
--- a/doc/specs/stdlib_linalg.md
+++ b/doc/specs/stdlib_linalg.md
@@ -168,3 +168,41 @@ program demo_trace
     print *, trace(A) ! 1 + 5 + 9
 end program demo_trace
 ```
+
+## `outer_product` - Computes the outer product of two vectors
+
+### Status
+
+Experimental
+
+### Description
+
+Computes the outer product of two vectors
+
+### Syntax
+
+`d = [[stdlib_linalg(module):outer_product(interface)]](u, v)`
+
+### Arguments
+
+`u`: Shall be a rank-1 array
+
+`v`: Shall be a rank-1 array
+
+### Return value
+
+Returns a rank-2 array equal to `u v^T` (where `u, v` are considered column vectors). The shape of the returned array is `[size(u), size(v)]`.
+
+### Example
+
+```fortran
+program demo_outer_product
+    use stdlib_linalg, only: outer_product
+    implicit none
+    real, allocatable :: A(:,:), u(:), v(:)
+    u = [1., 2., 3. ]
+    v = [3., 4.]
+    A = outer_product(u,v)
+    !A = reshape([3., 6., 9., 4., 8., 12.], [3,2])
+end program demo_outer_product
+```
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fa7f6e639..ef1d7963c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -9,6 +9,7 @@ set(fppFiles
     stdlib_io.fypp
     stdlib_linalg.fypp
     stdlib_linalg_diag.fypp
+    stdlib_linalg_outer_product.fypp
     stdlib_optval.fypp
     stdlib_sorting.fypp
     stdlib_sorting_ord_sort.fypp
diff --git a/src/Makefile.manual b/src/Makefile.manual
index e55bb80fd..f869cbd7f 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -6,6 +6,7 @@ SRCFYPP =\
         stdlib_io.fypp \
         stdlib_linalg.fypp \
         stdlib_linalg_diag.fypp \
+        stdlib_linalg_outer_product.fypp \
         stdlib_optval.fypp \
         stdlib_quadrature.fypp \
         stdlib_quadrature_trapz.fypp \
@@ -128,3 +129,4 @@ stdlib_stats_distribution_PRNG.o: \
 stdlib_string_type.o: stdlib_ascii.o stdlib_kinds.o
 stdlib_strings.o: stdlib_ascii.o stdlib_string_type.o
 stdlib_math.o: stdlib_kinds.o
+stdlib_linalg_outer_product.o: stdlib_linalg.o
diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp
index 51c2cdd54..5e0388c0b 100644
--- a/src/stdlib_linalg.fypp
+++ b/src/stdlib_linalg.fypp
@@ -11,6 +11,7 @@ module stdlib_linalg
   public :: diag
   public :: eye
   public :: trace
+  public :: outer_product
 
   interface diag
     !! version: experimental
@@ -52,6 +53,7 @@ module stdlib_linalg
     #:endfor
   end interface
 
+
   ! Matrix trace
   interface trace
     !! version: experimental
@@ -63,6 +65,21 @@ module stdlib_linalg
     #:endfor
   end interface
 
+
+  ! Outer product (of two vectors)
+  interface outer_product
+    !! version: experimental
+    !!
+    !! Computes the outer product of two vectors, returning a rank-2 array
+    !! ([Specification](../page/specs/stdlib_linalg.html#description_3))
+    #:for k1, t1 in RCI_KINDS_TYPES
+      pure module function outer_product_${t1[0]}$${k1}$(u, v) result(res)
+        ${t1}$, intent(in) :: u(:), v(:)
+        ${t1}$ :: res(size(u),size(v))
+      end function outer_product_${t1[0]}$${k1}$
+    #:endfor
+  end interface outer_product
+
 contains
 
     function eye(n) result(res)
diff --git a/src/stdlib_linalg_outer_product.fypp b/src/stdlib_linalg_outer_product.fypp
new file mode 100644
index 000000000..26c726435
--- /dev/null
+++ b/src/stdlib_linalg_outer_product.fypp
@@ -0,0 +1,20 @@
+#:include "common.fypp"
+#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
+submodule (stdlib_linalg) stdlib_linalg_outer_product
+
+  implicit none
+
+contains
+
+  #:for k1, t1 in RCI_KINDS_TYPES
+    pure module function outer_product_${t1[0]}$${k1}$(u, v) result(res)
+      ${t1}$, intent(in) :: u(:), v(:)
+      ${t1}$ :: res(size(u),size(v))
+      integer :: col
+      do col = 1, size(v)
+        res(:,col) = v(col) * u
+      end do
+    end function outer_product_${t1[0]}$${k1}$
+  #:endfor
+
+end submodule
diff --git a/src/tests/linalg/test_linalg.f90 b/src/tests/linalg/test_linalg.f90
index 4ad178d5b..cc8d0db68 100644
--- a/src/tests/linalg/test_linalg.f90
+++ b/src/tests/linalg/test_linalg.f90
@@ -2,7 +2,7 @@ program test_linalg
   
   use stdlib_error, only: check
   use stdlib_kinds, only: sp, dp, qp, int8, int16, int32, int64
-  use stdlib_linalg, only: diag, eye, trace
+  use stdlib_linalg, only: diag, eye, trace, outer_product
   
   implicit none
   
@@ -56,6 +56,22 @@ program test_linalg
   call test_trace_int32
   call test_trace_int64
 
+  !
+  ! outer product
+  !
+  call test_outer_product_rsp
+  call test_outer_product_rdp
+  call test_outer_product_rqp
+
+  call test_outer_product_csp
+  call test_outer_product_cdp
+  call test_outer_product_cqp
+
+  call test_outer_product_int8
+  call test_outer_product_int16
+  call test_outer_product_int32
+  call test_outer_product_int64
+
 
 contains
 
@@ -75,7 +91,7 @@ subroutine test_eye
     cye = eye(7)
     call check(abs(trace(cye) - cmplx(7.0_sp,0.0_sp,kind=sp)) < sptol, &
           msg="abs(trace(cye) - cmplx(7.0_sp,0.0_sp,kind=sp)) < sptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_eye
 
   subroutine test_diag_rsp
     integer, parameter :: n = 3
@@ -90,7 +106,7 @@ subroutine test_diag_rsp
 
     call check(all(diag(3*a) == 3*v), &
             msg="all(diag(3*a) == 3*v) failed.",warn=warn)
-  end subroutine
+  end subroutine test_diag_rsp
 
   subroutine test_diag_rsp_k
     integer, parameter :: n = 4
@@ -118,7 +134,7 @@ subroutine test_diag_rsp_k
     end do
     call check(size(diag(a,n+1)) == 0, &
         msg="size(diag(a,n+1)) == 0 failed.",warn=warn)
-  end subroutine
+  end subroutine test_diag_rsp_k
 
   subroutine test_diag_rdp
     integer, parameter :: n = 3
@@ -133,7 +149,7 @@ subroutine test_diag_rdp
 
     call check(all(diag(3*a) == 3*v), &
             msg="all(diag(3*a) == 3*v) failed.",warn=warn)
-  end subroutine
+  end subroutine test_diag_rdp
 
   subroutine test_diag_rqp
     integer, parameter :: n = 3
@@ -148,7 +164,7 @@ subroutine test_diag_rqp
 
     call check(all(diag(3*a) == 3*v), &
             msg="all(diag(3*a) == 3*v) failed.", warn=warn)
-  end subroutine
+  end subroutine test_diag_rqp
 
   subroutine test_diag_csp
     integer, parameter :: n = 3
@@ -165,7 +181,7 @@ subroutine test_diag_csp
           msg="all(abs(real(diag(a)) - [(i,i=1,n)]) < sptol)", warn=warn)
     call check(all(abs(aimag(diag(a)) - [(1,i=1,n)]) < sptol), &
           msg="all(abs(aimag(diag(a)) - [(1,i=1,n)]) < sptol)", warn=warn)
-  end subroutine
+  end subroutine test_diag_csp
 
   subroutine test_diag_cdp
     integer, parameter :: n = 3
@@ -175,7 +191,7 @@ subroutine test_diag_cdp
     a = diag([i_],-2) + diag([i_],2)
     call check(a(3,1) == i_ .and. a(1,3) == i_, &
           msg="a(3,1) == i_ .and. a(1,3) == i_ failed.",warn=warn)
-  end subroutine
+  end subroutine test_diag_cdp
 
   subroutine test_diag_cqp
     integer, parameter :: n = 3
@@ -185,7 +201,7 @@ subroutine test_diag_cqp
     a = diag([i_,i_],-1) + diag([i_,i_],1)
     call check(all(diag(a,-1) == i_) .and. all(diag(a,1) == i_), &
           msg="all(diag(a,-1) == i_) .and. all(diag(a,1) == i_) failed.",warn=warn)
-  end subroutine
+  end subroutine test_diag_cqp
 
   subroutine test_diag_int8
     integer, parameter :: n = 3
@@ -199,7 +215,7 @@ subroutine test_diag_int8
           msg="all(diag(a) == pack(a,mask)) failed.", warn=warn)
     call check(all(diag(diag(a)) == merge(a,0_int8,mask)), &
           msg="all(diag(diag(a)) == merge(a,0_int8,mask)) failed.", warn=warn)
-  end subroutine
+  end subroutine test_diag_int8
   subroutine test_diag_int16
     integer, parameter :: n = 4
     integer(int16), allocatable :: a(:,:)
@@ -212,7 +228,7 @@ subroutine test_diag_int16
           msg="all(diag(a) == pack(a,mask))", warn=warn)
     call check(all(diag(diag(a)) == merge(a,0_int16,mask)), &
           msg="all(diag(diag(a)) == merge(a,0_int16,mask)) failed.", warn=warn)
-  end subroutine
+  end subroutine test_diag_int16
   subroutine test_diag_int32
     integer, parameter :: n = 3
     integer(int32) :: a(n,n)
@@ -226,7 +242,7 @@ subroutine test_diag_int32
           msg="all(diag([1,1],-1) == a) failed.", warn=warn)
     call check(all(diag([1,1],1) == transpose(a)), &
           msg="all(diag([1,1],1) == transpose(a)) failed.", warn=warn)
-  end subroutine
+  end subroutine test_diag_int32
   subroutine test_diag_int64
     integer, parameter :: n = 4
     integer(int64) :: a(n,n), c(0:2*n-1)
@@ -257,7 +273,7 @@ subroutine test_diag_int64
     end do
     call check(all(diag(a,-2) == diag(a,2)), &
           msg="all(diag(a,-2) == diag(a,2))", warn=warn)
-  end subroutine
+  end subroutine test_diag_int64
 
 
 
@@ -270,7 +286,7 @@ subroutine test_trace_rsp
     a = reshape([(i,i=1,n**2)],[n,n])
     call check(abs(trace(a) - sum(diag(a))) < sptol, &
           msg="abs(trace(a) - sum(diag(a))) < sptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_rsp
 
   subroutine test_trace_rsp_nonsquare
     integer, parameter :: n = 4
@@ -287,7 +303,7 @@ subroutine test_trace_rsp_nonsquare
 
     call check(abs(trace(a) - ans) < sptol, &
           msg="abs(trace(a) - ans) < sptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_rsp_nonsquare
 
   subroutine test_trace_rdp
     integer, parameter :: n = 4
@@ -297,7 +313,7 @@ subroutine test_trace_rdp
     a = reshape([(i,i=1,n**2)],[n,n])
     call check(abs(trace(a) - sum(diag(a))) < dptol, &
           msg="abs(trace(a) - sum(diag(a))) < dptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_rdp
 
   subroutine test_trace_rdp_nonsquare
     integer, parameter :: n = 4
@@ -314,7 +330,7 @@ subroutine test_trace_rdp_nonsquare
 
     call check(abs(trace(a) - ans) < dptol, &
           msg="abs(trace(a) - ans) < dptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_rdp_nonsquare
 
   subroutine test_trace_rqp
     integer, parameter :: n = 3
@@ -324,7 +340,7 @@ subroutine test_trace_rqp
     a = reshape([(i,i=1,n**2)],[n,n])
     call check(abs(trace(a) - sum(diag(a))) < qptol, &
           msg="abs(trace(a) - sum(diag(a))) < qptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_rqp
 
 
   subroutine test_trace_csp
@@ -345,7 +361,7 @@ subroutine test_trace_csp
     ! tr(A + B) = tr(A) + tr(B)
     call check(abs(trace(a+b) - (trace(a) + trace(b))) < sptol, &
           msg="abs(trace(a+b) - (trace(a) + trace(b))) < sptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_csp
 
   subroutine test_trace_cdp
     integer, parameter :: n = 3
@@ -359,7 +375,7 @@ subroutine test_trace_cdp
 
     call check(abs(trace(a) - ans) < dptol, &
           msg="abs(trace(a) - ans) < dptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_cdp
 
   subroutine test_trace_cqp
     integer, parameter :: n = 3
@@ -369,7 +385,7 @@ subroutine test_trace_cqp
     a = 3*eye(n) + 4*eye(n)*i_ ! pythagorean triple
     call check(abs(trace(a)) - 3*5.0_qp < qptol, &
           msg="abs(trace(a)) - 3*5.0_qp < qptol failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_cqp
 
 
   subroutine test_trace_int8
@@ -380,7 +396,7 @@ subroutine test_trace_int8
     a = reshape([(i**2,i=1,n**2)],[n,n])
     call check(trace(a) == (1 + 25 + 81), &
           msg="trace(a) == (1 + 25 + 81) failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_int8
 
   subroutine test_trace_int16
     integer, parameter :: n = 3
@@ -390,7 +406,7 @@ subroutine test_trace_int16
     a = reshape([(i**3,i=1,n**2)],[n,n])
     call check(trace(a) == (1 + 125 + 729), &
           msg="trace(a) == (1 + 125 + 729) failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_int16
 
   subroutine test_trace_int32
     integer, parameter :: n = 3
@@ -400,7 +416,7 @@ subroutine test_trace_int32
     a = reshape([(i**4,i=1,n**2)],[n,n])
     call check(trace(a) == (1 + 625 + 6561), &
           msg="trace(a) == (1 + 625 + 6561) failed.",warn=warn)
-  end subroutine
+  end subroutine test_trace_int32
 
   subroutine test_trace_int64
     integer, parameter :: n = 5
@@ -424,7 +440,129 @@ subroutine test_trace_int64
     call check(trace(h) == sum(c(0:nd:2)), &
           msg="trace(h) == sum(c(0:nd:2)) failed.",warn=warn)
 
-  end subroutine
+  end subroutine test_trace_int64
+
+
+  subroutine test_outer_product_rsp
+    integer, parameter :: n = 2
+    real(sp) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_rsp" 
+    u = [1.,2.]
+    v = [1.,3.]
+    expected = reshape([1.,2.,3.,6.],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) < sptol), &
+         msg="all(abs(diff) < sptol) failed.",warn=warn)
+  end subroutine test_outer_product_rsp
+
+  subroutine test_outer_product_rdp
+    integer, parameter :: n = 2
+    real(dp) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_rdp" 
+    u = [1.,2.]
+    v = [1.,3.]
+    expected = reshape([1.,2.,3.,6.],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) < dptol), &
+         msg="all(abs(diff) < dptol) failed.",warn=warn)
+  end subroutine test_outer_product_rdp
+
+  subroutine test_outer_product_rqp
+    integer, parameter :: n = 2
+    real(qp) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_rqp" 
+    u = [1.,2.]
+    v = [1.,3.]
+    expected = reshape([1.,2.,3.,6.],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) < qptol), &
+         msg="all(abs(diff) < qptol) failed.",warn=warn)
+  end subroutine test_outer_product_rqp
+
+  subroutine test_outer_product_csp
+    integer, parameter :: n = 2
+    complex(sp) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_csp" 
+    u = [cmplx(1.,1.),cmplx(2.,0.)]
+    v = [cmplx(1.,0.),cmplx(3.,1.)]
+    expected = reshape([cmplx(1.,1.),cmplx(2.,0.),cmplx(2.,4.),cmplx(6.,2.)],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) < sptol), &
+         msg="all(abs(diff) < sptol) failed.",warn=warn)
+  end subroutine test_outer_product_csp
+
+  subroutine test_outer_product_cdp
+    integer, parameter :: n = 2
+    complex(dp) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_cdp" 
+    u = [cmplx(1.,1.),cmplx(2.,0.)]
+    v = [cmplx(1.,0.),cmplx(3.,1.)]
+    expected = reshape([cmplx(1.,1.),cmplx(2.,0.),cmplx(2.,4.),cmplx(6.,2.)],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) < dptol), &
+         msg="all(abs(diff) < dptol) failed.",warn=warn)
+  end subroutine test_outer_product_cdp
+
+  subroutine test_outer_product_cqp
+    integer, parameter :: n = 2
+    complex(qp) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_cqp" 
+    u = [cmplx(1.,1.),cmplx(2.,0.)]
+    v = [cmplx(1.,0.),cmplx(3.,1.)]
+    expected = reshape([cmplx(1.,1.),cmplx(2.,0.),cmplx(2.,4.),cmplx(6.,2.)],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) < qptol), &
+         msg="all(abs(diff) < qptol) failed.",warn=warn)
+  end subroutine test_outer_product_cqp
+    
+  subroutine test_outer_product_int8
+    integer, parameter :: n = 2
+    integer(int8) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_int8" 
+    u = [1,2]
+    v = [1,3]
+    expected = reshape([1,2,3,6],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) == 0), &
+         msg="all(abs(diff) == 0) failed.",warn=warn)
+  end subroutine test_outer_product_int8
+
+  subroutine test_outer_product_int16
+    integer, parameter :: n = 2
+    integer(int16) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_int16" 
+    u = [1,2]
+    v = [1,3]
+    expected = reshape([1,2,3,6],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) == 0), &
+         msg="all(abs(diff) == 0) failed.",warn=warn)
+  end subroutine test_outer_product_int16
+
+  subroutine test_outer_product_int32
+    integer, parameter :: n = 2
+    integer(int32) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_int32" 
+    u = [1,2]
+    v = [1,3]
+    expected = reshape([1,2,3,6],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) == 0), &
+         msg="all(abs(diff) == 0) failed.",warn=warn)
+  end subroutine test_outer_product_int32
+
+  subroutine test_outer_product_int64
+    integer, parameter :: n = 2
+    integer(int64) :: u(n), v(n), expected(n,n), diff(n,n)
+    write(*,*) "test_outer_product_int64" 
+    u = [1,2]
+    v = [1,3]
+    expected = reshape([1,2,3,6],[n,n])
+    diff = expected - outer_product(u,v)
+    call check(all(abs(diff) == 0), &
+         msg="all(abs(diff) == 0) failed.",warn=warn)
+  end subroutine test_outer_product_int64
+
 
   pure recursive function catalan_number(n) result(value)
     integer, intent(in) :: n