1 ## Copyright (C) 1995-2012 A. Scottedward Hodel
3 ## This file is part of Octave.
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 3 of the License, or (at
8 ## your option) any later version.
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
20 ## @deftypefn {Function File} {[@var{housv}, @var{beta}, @var{zer}] =} housh (@var{x}, @var{j}, @var{z})
21 ## Compute Householder reflection vector @var{housv} to reflect @var{x}
22 ## to be the j-th column of identity, i.e.,
26 ## (I - beta*housv*housv')x = norm(x)*e(j) if x(j) < 0,
27 ## (I - beta*housv*housv')x = -norm(x)*e(j) if x(j) >= 0
42 ## threshold for zero (usually should be the number 0)
46 ## Outputs (see Golub and Van Loan):
50 ## If beta = 0, then no reflection need be applied (zer set to 0)
57 ## Author: A. S. Hodel
58 ## Created: August 1995
60 function [housv, beta, zer] = housh (x, j, z)
66 ## Check for valid inputs.
67 if (! isvector (x) && ! isscalar (x))
68 error ("housh: first input must be a vector");
69 elseif (! isscalar(j))
70 error ("housh: second argment must be an integer scalar");
73 m = max (abs (housv));
78 beta = 1.0 / (alpha * (alpha + abs (housv(j))));
83 housv(j) = housv(j) + alpha*sg;
98 %! [hv, b, z] = housh(x, j, 0);
99 %! r = (eye(3) - b*hv*hv') * x;
100 %! d = - norm(x) * [0 0 1]';
101 %! assert(r, d, 2e-8);
102 %! assert(z, 0, 2e-8);
107 %! [hv, b, z] = housh(x, j, 0);
108 %! r = (eye(3) - b*hv*hv') * x;
109 %! d = norm(x) * [0 1 0]';
110 %! assert(r, d, 2e-8);
111 %! assert(z, 0, 2e-8);
116 %! [hv, b, z] = housh(x, j, 10);
117 %! r = (eye(3) - b*hv*hv') * x;
118 %! d = norm(x) * [1 0 0]';
119 %! assert(r, d, 2e-8);
120 %! assert(z, 1, 2e-8);
125 %! [hv, b, z] = housh(x, j, 0);
126 %! r = (eye(4) - b*hv*hv') * x;
127 %! d = - norm(x) * [0 1 0 0]';
128 %! assert(r, d, 2e-8);
129 %! assert(z, 0, 2e-8);