|
3211
|
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 |
|
3284
|
17 # Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
3211
|
18 |
|
|
19 function [housv,beta,zer] = housh(x,j,z) |
|
|
20 # function [housv,beta,zer] = housh(x,j,z) |
|
|
21 # computes householder reflection vector housv to reflect x to be |
|
|
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 |
|
|
30 # mar 6,1987 : rev dec 17,1988 |
|
|
31 # rev sep 19,1991 (blas) |
|
|
32 # translated from FORTRAN Aug 1995 |
|
|
33 # A. S. Hodel |
|
|
34 |
|
|
35 # $Revision: 1.1.1.1 $ |
|
|
36 # $Log$ |
|
|
37 |
|
|
38 # check for legal inputs |
|
|
39 if( !is_vector(x) && !is_scalar(x)) |
|
|
40 error("housh: first input must be a vector") |
|
|
41 elseif( !is_scalar(j) ) |
|
|
42 error("housh: second argment must be an integer scalar") |
|
|
43 else |
|
|
44 housv = x; |
|
|
45 m = max(abs(housv)); |
|
|
46 if (m ~= 0.0) |
|
|
47 housv = housv/m; |
|
|
48 alpha = norm(housv); |
|
|
49 if (alpha > z) |
|
|
50 beta = 1.0/(alpha*(alpha+abs(housv(j)))); |
|
|
51 sg = sign(housv(j)); |
|
|
52 if( sg == 0) |
|
|
53 sg = 1; |
|
|
54 endif |
|
|
55 housv(j) = housv(j) + alpha*sg; |
|
|
56 else |
|
|
57 beta = 0.0; |
|
|
58 endif |
|
|
59 else |
|
|
60 beta = 0.0; |
|
|
61 endif |
|
|
62 zer = (beta == 0); |
|
|
63 endif |
|
|
64 endfunction |
