|
3427
|
1 ## Copyright (C) 1995, 1998 A. Scottedward Hodel |
|
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
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. |
|
3427
|
9 ## |
|
7016
|
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. |
|
3427
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
3211
|
18 |
|
3439
|
19 ## -*- texinfo -*- |
|
|
20 ## @deftypefn {Function File} {[@var{housv}, @var{beta}, @var{zer}] =} housh (@var{x}, @var{j}, @var{z}) |
|
|
21 ## Computes householder reflection vector housv to reflect x to be |
|
3427
|
22 ## jth column of identity, i.e., (I - beta*housv*housv')x =e(j) |
|
|
23 ## inputs |
|
|
24 ## x: vector |
|
|
25 ## j: index into vector |
|
|
26 ## z: threshold for zero (usually should be the number 0) |
|
|
27 ## outputs: (see Golub and Van Loan) |
|
|
28 ## beta: If beta = 0, then no reflection need be applied (zer set to 0) |
|
|
29 ## housv: householder vector |
|
3439
|
30 ## @end deftypefn |
|
3211
|
31 |
|
3427
|
32 ## Author: A. S. Hodel |
|
3439
|
33 ## Created: August 1995 |
|
3427
|
34 |
|
|
35 function [housv, beta, zer] = housh (x, j, z) |
|
3211
|
36 |
|
3427
|
37 ## check for valid inputs |
|
4030
|
38 if (!isvector (x) && ! isscalar (x)) |
|
3427
|
39 error ("housh: first input must be a vector") |
|
4030
|
40 elseif (! isscalar(j)) |
|
3427
|
41 error ("housh: second argment must be an integer scalar") |
|
3211
|
42 else |
|
|
43 housv = x; |
|
3427
|
44 m = max (abs (housv)); |
|
|
45 if (m != 0.0) |
|
|
46 housv = housv / m; |
|
|
47 alpha = norm (housv); |
|
3426
|
48 if (alpha > z) |
|
3427
|
49 beta = 1.0 / (alpha * (alpha + abs (housv(j)))); |
|
|
50 sg = sign (housv(j)); |
|
|
51 if (sg == 0) |
|
3211
|
52 sg = 1; |
|
|
53 endif |
|
|
54 housv(j) = housv(j) + alpha*sg; |
|
|
55 else |
|
|
56 beta = 0.0; |
|
|
57 endif |
|
|
58 else |
|
|
59 beta = 0.0; |
|
|
60 endif |
|
|
61 zer = (beta == 0); |
|
|
62 endif |
|
3427
|
63 |
|
3211
|
64 endfunction |
