update packages

[CreaPhase.git] / octave_packages / m / linear-algebra / housh.m

diff --git a/octave_packages/m/linear-algebra/housh.m b/octave_packages/m/linear-algebra/housh.m
new file mode 100644 (file)
index 0000000..54246c7
--- /dev/null
+++ b/octave_packages/m/linear-algebra/housh.m
@@ -0,0 +1,133 @@
+## Copyright (C) 1995-2012 A. Scottedward Hodel
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{housv}, @var{beta}, @var{zer}] =} housh (@var{x}, @var{j}, @var{z})
+## Compute Householder reflection vector @var{housv} to reflect @var{x}
+## to be the j-th column of identity, i.e.,
+##
+## @example
+## @group
+## (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
+##
+## @noindent
+## Inputs
+##
+## @table @var
+## @item x
+## vector
+##
+## @item j
+## index into vector
+##
+## @item z
+## threshold for zero  (usually should be the number 0)
+## @end table
+##
+## @noindent
+## Outputs (see Golub and Van Loan):
+##
+## @table @var
+## @item beta
+## If beta = 0, then no reflection need be applied (zer set to 0)
+##
+## @item housv
+## householder vector
+## @end table
+## @end deftypefn
+
+## Author: A. S. Hodel
+## Created: August 1995
+
+function [housv, beta, zer] = housh (x, j, z)
+
+  if (nargin != 3)
+    print_usage ();
+  endif
+
+  ## Check for valid inputs.
+  if (! isvector (x) && ! isscalar (x))
+    error ("housh: first input must be a vector");
+  elseif (! isscalar(j))
+    error ("housh: second argment must be an integer scalar");
+  else
+    housv = x;
+    m = max (abs (housv));
+    if (m != 0.0)
+      housv = housv / m;
+      alpha = norm (housv);
+      if (alpha > z)
+        beta = 1.0 / (alpha * (alpha + abs (housv(j))));
+        sg = sign (housv(j));
+        if (sg == 0)
+          sg = 1;
+        endif
+        housv(j) = housv(j) + alpha*sg;
+      else
+        beta = 0.0;
+      endif
+    else
+      beta = 0.0;
+    endif
+    zer = (beta == 0);
+  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();
+