Mercurial > hg > octave-lyh
changeset 13048:c5c94b63931f
codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref
author | Roman Belov <romblv@gmail.com> |
---|---|
date | Sat, 03 Sep 2011 18:40:46 +0400 |
parents | 69a4609e61e2 |
children | b37d8e5aedf3 |
files | scripts/linear-algebra/cross.m scripts/linear-algebra/housh.m scripts/linear-algebra/planerot.m scripts/linear-algebra/qzhess.m scripts/linear-algebra/rref.m |
diffstat | 5 files changed, 165 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- 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(); +