--- a/src/DLD-FUNCTIONS/pinv.cc
+++ b/src/DLD-FUNCTIONS/pinv.cc
@@ -169,3 +169,23 @@
 
   return retval;
 }
+
+/*
+%!shared a, b, hitol, d, u, x, y, tolmult
+%! a = reshape (rand*[1:16], 4, 4);   ## Rank 2 matrix
+%! b = pinv (a);
+%! hitol = sqrt(eps);
+%! d = diag ([rand, rand, hitol, hitol]);
+%! u = rand (4);                      ## Could be singular by freak accident
+%! x = inv (u)*d*u;
+%! y = pinv (x, sqrt(eps));
+%! tolmult = 15;
+%!assert(a*b*a, a, tolmult*eps);
+%!assert(b*a*b, b, tolmult*eps);
+%!assert((b*a)', b*a, tolmult*eps);
+%!assert((a*b)', a*b, tolmult*eps);
+%!assert(x*y*x, x, tolmult*hitol);
+%!assert(y*x*y, y, tolmult*hitol);
+%!assert((x*y)', x*y, tolmult*hitol);
+%!assert((y*x)', y*x, tolmult*hitol);
+*/