|
8920
|
1 ## Copyright (C) 1995, 1998, 1999, 2000, 2002, 2005, 2007, 2009 |
|
7017
|
2 ## A. Scottedward Hodel |
|
3427
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
3427
|
10 ## |
|
7016
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
3427
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
3211
|
19 |
|
3439
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {[@var{housv}, @var{beta}, @var{zer}] =} housh (@var{x}, @var{j}, @var{z}) |
|
7201
|
22 ## Compute Householder reflection vector @var{housv} to reflect @var{x} |
|
8494
|
23 ## to be the j-th column of identity, i.e., |
|
7201
|
24 ## |
|
|
25 ## @example |
|
|
26 ## @group |
|
|
27 ## (I - beta*housv*housv')x = norm(x)*e(j) if x(1) < 0, |
|
|
28 ## (I - beta*housv*housv')x = -norm(x)*e(j) if x(1) >= 0 |
|
|
29 ## @end group |
|
|
30 ## @end example |
|
|
31 ## |
|
|
32 ## @noindent |
|
|
33 ## Inputs |
|
|
34 ## |
|
|
35 ## @table @var |
|
|
36 ## @item x |
|
|
37 ## vector |
|
|
38 ## @item j |
|
|
39 ## index into vector |
|
|
40 ## @item z |
|
|
41 ## threshold for zero (usually should be the number 0) |
|
|
42 ## @end table |
|
|
43 ## |
|
|
44 ## @noindent |
|
|
45 ## Outputs (see Golub and Van Loan): |
|
|
46 ## |
|
|
47 ## @table @var |
|
|
48 ## @item beta |
|
|
49 ## If beta = 0, then no reflection need be applied (zer set to 0) |
|
|
50 ## @item housv |
|
|
51 ## householder vector |
|
|
52 ## @end table |
|
3439
|
53 ## @end deftypefn |
|
3211
|
54 |
|
3427
|
55 ## Author: A. S. Hodel |
|
3439
|
56 ## Created: August 1995 |
|
3427
|
57 |
|
|
58 function [housv, beta, zer] = housh (x, j, z) |
|
3211
|
59 |
|
7125
|
60 if (nargin != 3) |
|
|
61 print_usage (); |
|
|
62 endif |
|
|
63 |
|
|
64 ## Check for valid inputs. |
|
|
65 if (! isvector (x) && ! isscalar (x)) |
|
8664
|
66 error ("housh: first input must be a vector"); |
|
4030
|
67 elseif (! isscalar(j)) |
|
8664
|
68 error ("housh: second argment must be an integer scalar"); |
|
3211
|
69 else |
|
|
70 housv = x; |
|
3427
|
71 m = max (abs (housv)); |
|
|
72 if (m != 0.0) |
|
|
73 housv = housv / m; |
|
|
74 alpha = norm (housv); |
|
3426
|
75 if (alpha > z) |
|
3427
|
76 beta = 1.0 / (alpha * (alpha + abs (housv(j)))); |
|
|
77 sg = sign (housv(j)); |
|
|
78 if (sg == 0) |
|
3211
|
79 sg = 1; |
|
|
80 endif |
|
|
81 housv(j) = housv(j) + alpha*sg; |
|
|
82 else |
|
|
83 beta = 0.0; |
|
|
84 endif |
|
|
85 else |
|
|
86 beta = 0.0; |
|
|
87 endif |
|
|
88 zer = (beta == 0); |
|
|
89 endif |
|
3427
|
90 |
|
3211
|
91 endfunction |
