|
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 |
|
|
6 ## under the terms of the GNU General Public License as published by the |
|
|
7 ## Free Software Foundation; either version 2, or (at your option) any |
|
|
8 ## later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
13 ## for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
3211
|
18 |
|
3427
|
19 ## function [housv,beta,zer] = housh(x,j,z) |
|
|
20 ## computes householder reflection vector housv to reflect x to be |
|
|
21 ## jth column of identity, i.e., (I - beta*housv*housv')x =e(j) |
|
|
22 ## inputs |
|
|
23 ## x: vector |
|
|
24 ## j: index into vector |
|
|
25 ## z: threshold for zero (usually should be the number 0) |
|
|
26 ## outputs: (see Golub and Van Loan) |
|
|
27 ## beta: If beta = 0, then no reflection need be applied (zer set to 0) |
|
|
28 ## housv: householder vector |
|
|
29 ## mar 6,1987 : rev dec 17,1988 |
|
|
30 ## rev sep 19,1991 (blas) |
|
|
31 ## translated from FORTRAN Aug 1995 |
|
3211
|
32 |
|
3427
|
33 ## Author: A. S. Hodel |
|
|
34 |
|
|
35 function [housv, beta, zer] = housh (x, j, z) |
|
3211
|
36 |
|
3427
|
37 ## check for valid inputs |
|
|
38 if (!is_vector (x) && ! is_scalar (x)) |
|
|
39 error ("housh: first input must be a vector") |
|
|
40 elseif (! is_scalar(j)) |
|
|
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 |
