octave-antonio: scripts/linear-algebra/housh.m annotate

annotate scripts/linear-algebra/housh.m @ 20170:af2b7695f1c4 draft default tip @

gallery.m: clean and vectorize part of qmult. * scripts/special-matrix/gallery.m: make qmult() follow Octave guidelines, vectorize last for loop in qmult().

author	Antonio Pino Robles <data.script93@gmail.com>
date	Thu, 28 May 2015 18:32:47 +0200
parents	4197fc428c7d
children	03b9d17a2d95

rev	line source
19697 4197fc428c7d maint: Update copyright notices for 2015. John W. Eaton <jwe@octave.org> parents: 19040 diff changeset	1 ## Copyright (C) 1995-2015 A. Scottedward Hodel
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	2 ##
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	3 ## This file is part of Octave.
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	4 ##
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	5 ## Octave is free software; you can redistribute it and/or modify it
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	6 ## under the terms of the GNU General Public License as published by
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	7 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	8 ## your option) any later version.
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	9 ##
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	10 ## Octave is distributed in the hope that it will be useful, but
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	13 ## General Public License for more details.
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	14 ##
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	15 ## You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	16 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5307 diff changeset	17 ## <http://www.gnu.org/licenses/>.
3211 440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	18
3439 3234a698073a [project @ 2000-01-14 09:51:14 by jwe] jwe parents: 3427 diff changeset	19 ## -- texinfo --
3234a698073a [project @ 2000-01-14 09:51:14 by jwe] jwe parents: 3427 diff changeset	20 ## @deftypefn {Function File} {[@var{housv}, @var{beta}, @var{zer}] =} housh (@var{x}, @var{j}, @var{z})
7201 76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	21 ## Compute Householder reflection vector @var{housv} to reflect @var{x}
8494 836618fee9d6 [docs] add hyphen Brian Gough <bjg@gnu.org> parents: 7201 diff changeset	22 ## to be the j-th column of identity, i.e.,
7201 76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	23 ##
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	24 ## @example
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	25 ## @group
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14363 diff changeset	26 ## (I - betahousvhousv')x = norm (x)*e(j) if x(j) < 0,
5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14363 diff changeset	27 ## (I - betahousvhousv')x = -norm (x)*e(j) if x(j) >= 0
7201 76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	28 ## @end group
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	29 ## @end example
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	30 ##
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	31 ## @noindent
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	32 ## Inputs
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	33 ##
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	34 ## @table @var
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	35 ## @item x
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	36 ## vector
10821 693e22af08ae Grammarcheck documentation of m-files Rik <octave@nomad.inbox5.com> parents: 8920 diff changeset	37 ##
7201 76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	38 ## @item j
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	39 ## index into vector
10821 693e22af08ae Grammarcheck documentation of m-files Rik <octave@nomad.inbox5.com> parents: 8920 diff changeset	40 ##
7201 76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	41 ## @item z
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	42 ## threshold for zero (usually should be the number 0)
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	43 ## @end table
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	44 ##
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	45 ## @noindent
19040 0850b5212619 doc: Add @nospell macro around proper names in documentation. Rik <rik@octave.org> parents: 17744 diff changeset	46 ## Outputs (see @nospell{Golub and Van Loan}):
7201 76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	47 ##
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	48 ## @table @var
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	49 ## @item beta
16826 a4969508008e doc: Periodic spellcheck of the documentation. Rik <rik@octave.org> parents: 14868 diff changeset	50 ## If beta = 0, then no reflection need be applied (@nospell{zer} set to 0)
10821 693e22af08ae Grammarcheck documentation of m-files Rik <octave@nomad.inbox5.com> parents: 8920 diff changeset	51 ##
7201 76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	52 ## @item housv
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	53 ## householder vector
76341ffda11e [project @ 2007-11-27 18:23:48 by jwe] jwe parents: 7125 diff changeset	54 ## @end table
3439 3234a698073a [project @ 2000-01-14 09:51:14 by jwe] jwe parents: 3427 diff changeset	55 ## @end deftypefn
3211 440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	56
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	57 ## Author: A. S. Hodel
3439 3234a698073a [project @ 2000-01-14 09:51:14 by jwe] jwe parents: 3427 diff changeset	58 ## Created: August 1995
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	59
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	60 function [housv, beta, zer] = housh (x, j, z)
3211 440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	61
7125 f084ba47812b [project @ 2007-11-08 02:29:23 by jwe] jwe parents: 7017 diff changeset	62 if (nargin != 3)
f084ba47812b [project @ 2007-11-08 02:29:23 by jwe] jwe parents: 7017 diff changeset	63 print_usage ();
f084ba47812b [project @ 2007-11-08 02:29:23 by jwe] jwe parents: 7017 diff changeset	64 endif
f084ba47812b [project @ 2007-11-08 02:29:23 by jwe] jwe parents: 7017 diff changeset	65
f084ba47812b [project @ 2007-11-08 02:29:23 by jwe] jwe parents: 7017 diff changeset	66 ## Check for valid inputs.
f084ba47812b [project @ 2007-11-08 02:29:23 by jwe] jwe parents: 7017 diff changeset	67 if (! isvector (x) && ! isscalar (x))
8664 e07e93c04080 style fixes John W. Eaton <jwe@octave.org> parents: 8494 diff changeset	68 error ("housh: first input must be a vector");
14868 5d3a684236b0 maint: Use Octave coding conventions for cuddling parentheses in scripts directory Rik <octave@nomad.inbox5.com> parents: 14363 diff changeset	69 elseif (! isscalar (j))
8664 e07e93c04080 style fixes John W. Eaton <jwe@octave.org> parents: 8494 diff changeset	70 error ("housh: second argment must be an integer scalar");
3211 440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	71 else
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	72 housv = x;
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	73 m = max (abs (housv));
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	74 if (m != 0.0)
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	75 housv = housv / m;
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	76 alpha = norm (housv);
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3284 diff changeset	77 if (alpha > z)
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	78 beta = 1.0 / (alpha * (alpha + abs (housv(j))));
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	79 sg = sign (housv(j));
e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	80 if (sg == 0)
3211 440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	81 sg = 1;
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	82 endif
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	83 housv(j) = housv(j) + alpha*sg;
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	84 else
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	85 beta = 0.0;
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	86 endif
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	87 else
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	88 beta = 0.0;
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	89 endif
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	90 zer = (beta == 0);
440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	91 endif
3427 e098ebb77023 [project @ 2000-01-13 09:25:53 by jwe] jwe parents: 3426 diff changeset	92
3211 440b2b28e74a [project @ 1998-11-05 04:16:22 by jwe] jwe parents: diff changeset	93 endfunction
13048 c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	94
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	95
13048 c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	96 %!test
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	97 %! x = [1 2 3]';
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	98 %! j = 3;
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	99 %! [hv, b, z] = housh (x, j, 0);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	100 %! r = (eye (3) - bhvhv') * x;
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	101 %! d = - norm (x) * [0 0 1]';
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	102 %! assert (r, d, 2e-8);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	103 %! assert (z, 0, 2e-8);
13048 c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	104
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	105 %!test
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	106 %! x = [7 -3 1]';
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	107 %! j = 2;
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	108 %! [hv, b, z] = housh (x, j, 0);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	109 %! r = (eye (3) - bhvhv') * x;
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	110 %! d = norm (x) * [0 1 0]';
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	111 %! assert (r, d, 2e-8);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	112 %! assert (z, 0, 2e-8);
13048 c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	113
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	114 %!test
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	115 %! x = [1 0 0]';
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	116 %! j = 1;
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	117 %! [hv, b, z] = housh (x, j, 10);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	118 %! r = (eye (3) - bhvhv') * x;
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	119 %! d = norm (x) * [1 0 0]';
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	120 %! assert (r, d, 2e-8);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	121 %! assert (z, 1, 2e-8);
13048 c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	122
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	123 %!test
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	124 %! x = [5 0 4 1]';
c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	125 %! j = 2;
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	126 %! [hv, b, z] = housh (x, j, 0);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	127 %! r = (eye (4) - bhvhv') * x;
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	128 %! d = - norm (x) * [0 1 0 0]';
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	129 %! assert (r, d, 2e-8);
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	130 %! assert (z, 0, 2e-8);
13048 c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	131
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	132 %!error housh ([0])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	133 %!error housh ()
13048 c5c94b63931f codesprint: linear algebra tests: cross, housh, planerot, qzhess, rref Roman Belov <romblv@gmail.com> parents: 11523 diff changeset	134

Mercurial > octave-antonio

annotate scripts/linear-algebra/housh.m @ 20170:af2b7695f1c4 draft default tip @