--- a/scripts/linear-algebra/cross.m
+++ b/scripts/linear-algebra/cross.m
@@ -90,3 +90,26 @@
   endif
 
 endfunction
+
+%!test
+%! x = [1 0 0];
+%! y = [0 1 0];
+%! r = [0 0 1];
+%! assert(cross(x, y), r, 2e-8);
+
+%!test
+%! x = [1 2 3];
+%! y = [4 5 6];
+%! r = [(2*6-3*5) (3*4-1*6) (1*5-2*4)];
+%! assert(cross(x, y), r, 2e-8);
+
+%!test
+%! x = [1 0 0; 0 1 0; 0 0 1];
+%! y = [0 1 0; 0 0 1; 1 0 0];
+%! r = [0 0 1; 1 0 0; 0 1 0];
+%! assert(cross(x, y, 2), r, 2e-8);
+%! assert(cross(x, y, 1), -r, 2e-8);
+
+%!error cross(0,0);
+%!error cross();
+

--- a/scripts/linear-algebra/housh.m
+++ b/scripts/linear-algebra/housh.m
@@ -23,8 +23,8 @@
 ##
 ## @example
 ## @group
-## (I - beta*housv*housv')x =  norm(x)*e(j) if x(1) < 0,
-## (I - beta*housv*housv')x = -norm(x)*e(j) if x(1) >= 0
+## (I - beta*housv*housv')x =  norm(x)*e(j) if x(j) < 0,
+## (I - beta*housv*housv')x = -norm(x)*e(j) if x(j) >= 0
 ## @end group
 ## @end example
 ##
@@ -91,3 +91,43 @@
   endif
 
 endfunction
+
+%!test
+%! x = [1 2 3]';
+%! j = 3;
+%! [hv, b, z] = housh(x, j, 0);
+%! r = (eye(3) - b*hv*hv') * x;
+%! d = - norm(x) * [0 0 1]';
+%! assert(r, d, 2e-8);
+%! assert(z, 0, 2e-8);
+
+%!test
+%! x = [7 -3 1]';
+%! j = 2;
+%! [hv, b, z] = housh(x, j, 0);
+%! r = (eye(3) - b*hv*hv') * x;
+%! d = norm(x) * [0 1 0]';
+%! assert(r, d, 2e-8);
+%! assert(z, 0, 2e-8);
+
+%!test
+%! x = [1 0 0]';
+%! j = 1;
+%! [hv, b, z] = housh(x, j, 10);
+%! r = (eye(3) - b*hv*hv') * x;
+%! d = norm(x) * [1 0 0]';
+%! assert(r, d, 2e-8);
+%! assert(z, 1, 2e-8);
+
+%!test
+%! x = [5 0 4 1]';
+%! j = 2;
+%! [hv, b, z] = housh(x, j, 0);
+%! r = (eye(4) - b*hv*hv') * x;
+%! d = - norm(x) * [0 1 0 0]';
+%! assert(r, d, 2e-8);
+%! assert(z, 0, 2e-8);
+
+%!error housh([0]);
+%!error housh();
+

--- a/scripts/linear-algebra/planerot.m
+++ b/scripts/linear-algebra/planerot.m
@@ -34,3 +34,14 @@
   G = givens (x(1), x(2));
   y = G * x(:);
 endfunction
+
+%!test
+%! x = [3 4];
+%! [g y] = planerot(x);
+%! assert(g - [x(1) x(2); -x(2) x(1)] / sqrt(x(1)^2 + x(2)^2), zeros(2), 2e-8);
+%! assert(y(2), 0, 2e-8);
+
+%!error planerot([0]);
+%!error planerot([0 0 0]);
+%!error planerot();
+

--- a/scripts/linear-algebra/qzhess.m
+++ b/scripts/linear-algebra/qzhess.m
@@ -90,3 +90,52 @@
   endfor
 
 endfunction
+
+%!test
+%! a = [1 2 1 3;
+%!      2 5 3 2;
+%!      5 5 1 0;
+%!      4 0 3 2];
+%! b = [0 4 2 1;
+%!      2 3 1 1;
+%!      1 0 2 1;
+%!      2 5 3 2];
+%! mask = [0 0 0 0;
+%!         0 0 0 0;
+%!         1 0 0 0;
+%!         1 1 0 0];
+%! [aa, bb, q, z] = qzhess(a, b);
+%! assert(inv(q) - q', zeros(4), 2e-8);
+%! assert(inv(z) - z', zeros(4), 2e-8);
+%! assert(q * a * z, aa, 2e-8);
+%! assert(aa .* mask, zeros(4), 2e-8);
+%! assert(q * b * z, bb, 2e-8);
+%! assert(bb .* mask, zeros(4), 2e-8);
+
+%!test
+%! a = [1 2 3 4 5;
+%!      3 2 3 1 0;
+%!      4 3 2 1 1;
+%!      0 1 0 1 0;
+%!      3 2 1 0 5];
+%! b = [5 0 4 0 1;
+%!      1 1 1 2 5;
+%!      0 3 2 1 0;
+%!      4 3 0 3 5;
+%!      2 1 2 1 3];
+%! mask = [0 0 0 0 0;
+%!         0 0 0 0 0;
+%!         1 0 0 0 0;
+%!         1 1 0 0 0;
+%!         1 1 1 0 0];
+%! [aa, bb, q, z] = qzhess(a, b);
+%! assert(inv(q) - q', zeros(5), 2e-8);
+%! assert(inv(z) - z', zeros(5), 2e-8);
+%! assert(q * a * z, aa, 2e-8);
+%! assert(aa .* mask, zeros(5), 2e-8);
+%! assert(q * b * z, bb, 2e-8);
+%! assert(bb .* mask, zeros(5), 2e-8);
+
+%!error qzhess([0]);
+%!error qzhess();
+

--- a/scripts/linear-algebra/rref.m
+++ b/scripts/linear-algebra/rref.m
@@ -86,3 +86,43 @@
   k = find (used);
 
 endfunction
+
+%!test
+%! a = [1];
+%! [r k] = rref(a);
+%! assert(r, [1], 2e-8);
+%! assert(k, [1], 2e-8);
+
+%!test
+%! a = [1 3; 4 5];
+%! [r k] = rref(a);
+%! assert(rank(a), rank(r), 2e-8);
+%! assert(r, eye(2), 2e-8);
+%! assert(k == [1, 2] || k == [2, 1]);
+
+
+%!test
+%! a = [1 3; 4 5; 7 9];
+%! [r k] = rref(a);
+%! assert(rank(a), rank(r), 2e-8);
+%! assert(r, eye(3)(:,1:2), 2e-8);
+%! assert(k, [1 2], 2e-8);
+
+%!test
+%! a = [1 2 3; 2 4 6; 7 2 0];
+%! [r k] = rref(a);
+%! assert(rank(a), rank(r), 2e-8);
+%! assert(r, [1 0 (3-7/2); 0 1 (7/4); 0 0 0], 2e-8);
+%! assert(k, [1 2], 2e-8);
+
+%!test
+%! a = [1 2 1; 2 4 2.01; 2 4 2.1];
+%! tol = 0.02;
+%! [r k] = rref(a, tol);
+%! assert(rank(a, tol), rank(r, tol), 2e-8);
+%! tol = 0.2;
+%! [r k] = rref(a, tol);
+%! assert(rank(a, tol), rank(r, tol), 2e-8);
+
+%!error rref();
+