diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt
index fc3f823d9..de9eb4f99 100644
--- a/example/linalg/CMakeLists.txt
+++ b/example/linalg/CMakeLists.txt
@@ -16,6 +16,7 @@ ADD_EXAMPLE(inverse_operator)
 ADD_EXAMPLE(inverse_function)
 ADD_EXAMPLE(inverse_inplace)
 ADD_EXAMPLE(inverse_subroutine)
+ADD_EXAMPLE(pseudoinverse)
 ADD_EXAMPLE(outer_product)
 ADD_EXAMPLE(eig)
 ADD_EXAMPLE(eigh)
diff --git a/example/linalg/example_pseudoinverse.f90 b/example/linalg/example_pseudoinverse.f90
new file mode 100644
index 000000000..6cf01642b
--- /dev/null
+++ b/example/linalg/example_pseudoinverse.f90
@@ -0,0 +1,32 @@
+! Matrix pseudo-inversion example: function, subroutine, and operator interfaces
+program example_pseudoinverse
+  use stdlib_linalg, only: pinv, pseudoinvert, operator(.pinv.), linalg_state_type
+  implicit none(type,external)
+
+  real :: A(15,5), Am1(5,15)
+  type(linalg_state_type) :: state
+  integer :: i, j
+  real, parameter :: tol = sqrt(epsilon(0.0))
+  
+  ! Generate random matrix A (15x15)
+  call random_number(A)
+
+  ! Pseudo-inverse: Function interfcae
+  Am1 = pinv(A, err=state)
+  print *, 'Max error (function)  : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
+  
+  ! User threshold
+  Am1 = pinv(A, rtol=0.001, err=state)
+  print *, 'Max error (rtol=0.001): ', maxval(abs(A-matmul(A, matmul(Am1,A))))  
+
+  ! Pseudo-inverse: Subroutine interface
+  call pseudoinvert(A, Am1, err=state)
+
+  print *, 'Max error (subroutine): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
+
+  ! Operator interface
+  Am1 = .pinv.A
+  
+  print *, 'Max error (operator)  : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
+
+end program example_pseudoinverse