Mercurial > octave-nkf
view scripts/linear-algebra/housh.m @ 7125:f084ba47812b
[project @ 2007-11-08 02:29:23 by jwe]
author | jwe |
---|---|
date | Thu, 08 Nov 2007 02:29:24 +0000 |
parents | a1dbe9d80eee |
children | 76341ffda11e |
line wrap: on
line source
## Copyright (C) 1995, 1998, 1999, 2000, 2002, 2005, 2007 ## A. Scottedward Hodel ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{housv}, @var{beta}, @var{zer}] =} housh (@var{x}, @var{j}, @var{z}) ## Computes householder reflection vector housv to reflect x to be ## jth column of identity, i.e., (I - beta*housv*housv')x =e(j) ## inputs ## x: vector ## j: index into vector ## z: threshold for zero (usually should be the number 0) ## outputs: (see Golub and Van Loan) ## beta: If beta = 0, then no reflection need be applied (zer set to 0) ## housv: householder vector ## @end deftypefn ## Author: A. S. Hodel ## Created: August 1995 function [housv, beta, zer] = housh (x, j, z) if (nargin != 3) print_usage (); endif ## Check for valid inputs. if (! isvector (x) && ! isscalar (x)) error ("housh: first input must be a vector") elseif (! isscalar(j)) error ("housh: second argment must be an integer scalar") else housv = x; m = max (abs (housv)); if (m != 0.0) housv = housv / m; alpha = norm (housv); if (alpha > z) beta = 1.0 / (alpha * (alpha + abs (housv(j)))); sg = sign (housv(j)); if (sg == 0) sg = 1; endif housv(j) = housv(j) + alpha*sg; else beta = 0.0; endif else beta = 0.0; endif zer = (beta == 0); endif endfunction