Mercurial > octave

annotate scripts/linear-algebra/normest1.m @ 27919:1891570abac8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2020.
author John W. Eaton <jwe@octave.org>
date Mon, 06 Jan 2020 22:29:51 -0500
parents b442ec6dda5c
children bd51beb6205e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27919
1891570abac8 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27918
diff changeset
1 ## Copyright (C) 2016-2020 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
2 ##
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
3 ## See the file COPYRIGHT.md in the top-level directory of this distribution
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
4 ## or <https://octave.org/COPYRIGHT.html/>.
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27800
diff changeset
5 ##
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
6 ##
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
7 ## This file is part of Octave.
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
8 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23618
diff changeset
9 ## Octave is free software: you can redistribute it and/or modify it
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
10 ## under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23618
diff changeset
11 ## the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22627
diff changeset
12 ## (at your option) any later version.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
13 ##
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
14 ## Octave is distributed in the hope that it will be useful, but
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
15 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22627
diff changeset
16 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22627
diff changeset
17 ## GNU General Public License for more details.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
18 ##
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
19 ## You should have received a copy of the GNU General Public License
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
20 ## along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23618
diff changeset
21 ## <https://www.gnu.org/licenses/>.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
22
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
23 ## -*- texinfo -*-
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
24 ## @deftypefn {} {@var{nest} =} normest1 (@var{A})
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
25 ## @deftypefnx {} {@var{nest} =} normest1 (@var{A}, @var{t})
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
26 ## @deftypefnx {} {@var{nest} =} normest1 (@var{A}, @var{t}, @var{x0})
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
27 ## @deftypefnx {} {@var{nest} =} normest1 (@var{Afun}, @var{t}, @var{x0}, @var{p1}, @var{p2}, @dots{})
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
28 ## @deftypefnx {} {[@var{nest}, @var{v}] =} normest1 (@var{A}, @dots{})
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
29 ## @deftypefnx {} {[@var{nest}, @var{v}, @var{w}] =} normest1 (@var{A}, @dots{})
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
30 ## @deftypefnx {} {[@var{nest}, @var{v}, @var{w}, @var{iter}] =} normest1 (@var{A}, @dots{})
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
31 ## Estimate the 1-norm of the matrix @var{A} using a block algorithm.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
32 ##
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
33 ## @code{normest1} is best for large sparse matrices where only an estimate of
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
34 ## the norm is required. For small to medium sized matrices, consider using
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
35 ## @code{norm (@var{A}, 1)}. In addition, @code{normest1} can be used for the
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
36 ## estimate of the 1-norm of a linear operator @var{A} when matrix-vector
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
37 ## products @code{@var{A} * @var{x}} and @code{@var{A}' * @var{x}} can be
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
38 ## cheaply computed. In this case, instead of the matrix @var{A}, a function
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
39 ## @code{@var{Afun} (@var{flag}, @var{x})} is used; it must return:
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
40 ##
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
41 ## @itemize @bullet
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
42 ## @item
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
43 ## the dimension @var{n} of @var{A}, if @var{flag} is @qcode{"dim"}
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22193
diff changeset
44 ##
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
45 ## @item
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
46 ## true if @var{A} is a real operator, if @var{flag} is @qcode{"real"}
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22193
diff changeset
47 ##
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
48 ## @item
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
49 ## the result @code{@var{A} * @var{x}}, if @var{flag} is @qcode{"notransp"}
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22193
diff changeset
50 ##
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
51 ## @item
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
52 ## the result @code{@var{A}' * @var{x}}, if @var{flag} is @qcode{"transp"}
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
53 ## @end itemize
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
54 ##
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
55 ## A typical case is @var{A} defined by @code{@var{b} ^ @var{m}}, in which the
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
56 ## result @code{@var{A} * @var{x}} can be computed without even forming
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
57 ## explicitly @code{@var{b} ^ @var{m}} by:
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
58 ##
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
59 ## @example
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22193
diff changeset
60 ## @group
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
61 ## @var{y} = @var{x};
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
62 ## for @var{i} = 1:@var{m}
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
63 ## @var{y} = @var{b} * @var{y};
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
64 ## endfor
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22193
diff changeset
65 ## @end group
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
66 ## @end example
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
67 ##
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
68 ## The parameters @var{p1}, @var{p2}, @dots{} are arguments of
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
69 ## @code{@var{Afun} (@var{flag}, @var{x}, @var{p1}, @var{p2}, @dots{})}.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
70 ##
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
71 ## The default value for @var{t} is 2. The algorithm requires matrix-matrix
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
72 ## products with sizes @var{n} x @var{n} and @var{n} x @var{t}.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
73 ##
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
74 ## The initial matrix @var{x0} should have columns of unit 1-norm. The default
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
75 ## initial matrix @var{x0} has the first column
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
76 ## @code{ones (@var{n}, 1) / @var{n}} and, if @var{t} > 1, the remaining
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
77 ## columns with random elements @code{-1 / @var{n}}, @code{1 / @var{n}},
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
78 ## divided by @var{n}.
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
79 ##
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
80 ## On output, @var{nest} is the desired estimate, @var{v} and @var{w}
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
81 ## are vectors such that @code{@var{w} = @var{A} * @var{v}}, with
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
82 ## @code{norm (@var{w}, 1)} = @code{@var{c} * norm (@var{v}, 1)}. @var{iter}
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
83 ## contains in @code{@var{iter}(1)} the number of iterations (the maximum is
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
84 ## hardcoded to 5) and in @code{@var{iter}(2)} the total number of products
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
85 ## @code{@var{A} * @var{x}} or @code{@var{A}' * @var{x}} performed by the
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
86 ## algorithm.
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
87 ##
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
88 ## Algorithm Note: @code{normest1} uses random numbers during evaluation.
22299
9fc91bb2aec3 doc: grammarcheck documentation for 4.2 release.
Rik <rik@octave.org>
parents: 22193
diff changeset
89 ## Therefore, if consistent results are required, the @qcode{"state"} of the
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
90 ## random generator should be fixed before invoking @code{normest1}.
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
91 ##
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
92 ## Reference: @nospell{N. J. Higham and F. Tisseur},
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
93 ## @cite{A block algorithm for matrix 1-norm estimation, with and
22330
53e246fd8124 doc: Spellcheck documentation ahead of 4.2 release.
Rik <rik@octave.org>
parents: 22299
diff changeset
94 ## application to 1-norm @nospell{pseudospectra}},
27800
5a6a19a4e3da doc: Use Texinfo non-sentence ending periods in citations.
Rik <rik@octave.org>
parents: 26376
diff changeset
95 ## @nospell{SIAM J. Matrix Anal.@: Appl.@:},
5a6a19a4e3da doc: Use Texinfo non-sentence ending periods in citations.
Rik <rik@octave.org>
parents: 26376
diff changeset
96 ## pp.@: 1185--1201, Vol 21, No.@: 4, 2000.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
97 ##
22827
c3d3a81ad986 doc: Update documentation for norm, normest, normest1, condest.
Rik <rik@octave.org>
parents: 22627
diff changeset
98 ## @seealso{normest, norm, cond, condest}
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
99 ## @end deftypefn
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
100
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
101 ## Ideally, we would set t and X to their default values but Matlab
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
102 ## compatibility would require we set the default even when they are empty.
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
103 function [nest, v, w, iter] = normest1 (A, t = [], x0 = [], varargin)
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
104
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
105 if (nargin < 1)
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
106 print_usage ();
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
107 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
108
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
109 if (isempty (t))
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
110 t = 2;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
111 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
112
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
113 ## FIXME: t < 0 should print trace information
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
114 if (isnumeric (A) && issquare (A))
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
115 Aismat = true;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
116 Aisreal = isreal (A);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
117 n = rows (A);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
118 if (n <= 4 || t == n)
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
119 ## small input, compute directly
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
120 [nest, idx] = max (sum (abs (A), 1), [] , 2);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
121 v = zeros (n, 1);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
122 v(idx) = 1;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
123 w = A(:, idx);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
124 ## Matlab incompatible on purpose. Matlab returns iter as a row vector
22193
e35866e6a2e0 normest1: always return a column vector for IT output (patch #8837)
Carnë Draug <carandraug@octave.org>
parents: 22183
diff changeset
125 ## for this special case, but a column vector in all other cases.
e35866e6a2e0 normest1: always return a column vector for IT output (patch #8837)
Carnë Draug <carandraug@octave.org>
parents: 22183
diff changeset
126 ## This is obviously a bug in Matlab that we don't reproduce.
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
127 iter = [0; 1];
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
128 return;
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
129 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
130 elseif (is_function_handle (A))
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
131 Aismat = false;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
132 Aisreal = A ("real", [], varargin{:});
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
133 n = A ("dim", [], varargin{:});
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
134 Afun = @(x) A ("notransp", x, varargin{:});
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
135 A1fun = @(x) A ("transp", x, varargin{:});
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
136 else
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
137 error ("normest1: A must be a square matrix or a function handle");
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
138 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
139
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
140 t = min (t, n);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
141
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
142 if (isempty (x0))
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
143 X = [ones(n, 1), sign(2 * rand(n, t - 1) - 1)];
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
144 i = 2;
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
145 imax = min (t, 2^(n-1));
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
146 ## There are at most 2^(n-1) unparallel columns, see later.
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
147 while (i <= imax)
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
148 if (any (abs (X(:,i)' * X(:,1:i-1)) == n))
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
149 ## column i is parallel to a colum 1:i-1. Change it.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
150 X(:,i) = sign (2 * rand (n, 1) - 1);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
151 else
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
152 i++;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
153 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
154 endwhile
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
155 X /= n;
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
156 else
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
157 if (columns (x0) < t)
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
158 error ("normest1: X0 must have %d columns", t);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
159 endif
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
160 X = x0;
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
161 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
162
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
163 itmax = 5;
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
164 idx_hist = zeros (n, 1);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
165 nest_old = 0;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
166 idx = zeros (n, 1);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
167 S = zeros (n, t);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
168 iter = [0; 0];
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
169 converged = false;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
170 while (! converged && (iter(1) < itmax))
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
171 iter(1)++;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
172 if (Aismat)
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
173 Y = A * X;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
174 else
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
175 Y = Afun (X);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
176 endif
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
177 iter(2)++;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
178 [nest, j] = max (sum (abs (Y), 1), [], 2);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
179 if ((nest > nest_old) || (iter(1) == 2))
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
180 idx_best = idx(j);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
181 w = Y(:, j); # there is an error in Algorithm 2.4
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
182 endif
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
183 if (nest <= nest_old && iter(1) >= 2) # (1) of Algorithm 2.4
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
184 nest = nest_old;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
185 break; # while
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
186 endif
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
187 nest_old = nest;
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
188 Sold = S;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
189 S = sign (Y);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
190 S(S==0) = 1;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
191 possible_break = false;
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
192 if (Aisreal)
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
193 ## test parallel (only real case)
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
194 if (all (any (abs (Sold' * S) == n))) # (2) of Algorithm 2.4
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
195 ## all columns of S parallel to a column of Sold, exit
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
196 possible_break = true;
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
197 converged = true;
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
198 else
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
199 if (t > 1)
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
200 ## at least two columns of S are not parallel
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
201 i = 1;
23617
4ce622b7b930 Fix possible infinite loop in normest1.m (bug #51241)
Marco Caliari <marco.caliari@univr.it>
parents: 23219
diff changeset
202 ## The maximum number of unparallel columns of length n with
4ce622b7b930 Fix possible infinite loop in normest1.m (bug #51241)
Marco Caliari <marco.caliari@univr.it>
parents: 23219
diff changeset
203 ## entries {-1,1} is 2^(n-1). n of them are already in Sold.
4ce622b7b930 Fix possible infinite loop in normest1.m (bug #51241)
Marco Caliari <marco.caliari@univr.it>
parents: 23219
diff changeset
204 imax = min (t, 2 ^ (n - 1) - n);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
205 while (i <= imax)
23617
4ce622b7b930 Fix possible infinite loop in normest1.m (bug #51241)
Marco Caliari <marco.caliari@univr.it>
parents: 23219
diff changeset
206 if (any (abs (S(:,i)' * S(:,1:i-1)) == n)
4ce622b7b930 Fix possible infinite loop in normest1.m (bug #51241)
Marco Caliari <marco.caliari@univr.it>
parents: 23219
diff changeset
207 || any (abs (S(:,i)' * Sold) == n))
4ce622b7b930 Fix possible infinite loop in normest1.m (bug #51241)
Marco Caliari <marco.caliari@univr.it>
parents: 23219
diff changeset
208 ## i-th column of S is parallel to a previous column
4ce622b7b930 Fix possible infinite loop in normest1.m (bug #51241)
Marco Caliari <marco.caliari@univr.it>
parents: 23219
diff changeset
209 ## or to a column of Sold. Change it.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
210 S(:,i) = sign (2*rand (n, 1)-1);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
211 else
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
212 i++;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
213 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
214 endwhile
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
215 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
216 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
217 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
218 if (! possible_break)
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
219 if (Aismat)
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
220 Z = A' * S;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
221 else
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
222 Z = A1fun (S); # (3) of Algorithm 2.4
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
223 endif
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
224 iter(2)++;
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
225 h = max (abs (Z), [], 2);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
226 idx = (1:n)';
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
227 if (iter(1) >= 2 && (max (h, [], 1) == h(idx_best))) # (4) of Alg. 2.4
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
228 break; # while
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
229 endif
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
230 [h, idx] = sort (h, "descend");
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
231 if (t > 1)
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
232 if (all (idx_hist(idx(1:t)))) # (5) of Algorithm 2.4
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
233 break; # while
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
234 endif
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
235 idx = idx(! idx_hist(idx));
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
236 ## length(idx) could be less than t, especially if t is not << n.
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
237 ## This is not considered in point (5) of Algorithm 2.4.
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
238 tmax = min (numel (idx), t);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
239 idx = idx(1:tmax);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
240 else
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
241 tmax = 1;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
242 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
243 X = zeros (n, tmax);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
244 X(sub2ind (size (X), idx(1:tmax), (1:tmax)')) = 1;
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
245 idx_hist(idx(1:tmax)) = 1;
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
246 endif
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
247 endwhile
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
248 v = zeros (n, 1);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
249 v(idx_best) = 1;
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
250
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
251 endfunction
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
252
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
253
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
254 %!function z = afun_A (flag, x, A, n)
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
255 %! switch (flag)
22627
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
256 %! case {"dim"}
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
257 %! z = n;
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
258 %! case {"real"}
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
259 %! z = isreal (A);
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
260 %! case {"transp"}
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
261 %! z = A' * x;
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
262 %! case {"notransp"}
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
263 %! z = A * x;
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
264 %! endswitch
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
265 %!endfunction
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
266 %!function z = afun_A_P (flag, x, A, m)
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
267 %! switch (flag)
22627
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
268 %! case "dim"
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
269 %! z = length (A);
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
270 %! case "real"
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
271 %! z = isreal (A);
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
272 %! case "transp"
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
273 %! z = x; for i = 1:m, z = A' * z;, endfor
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
274 %! case "notransp"
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
275 %! z = x; for i = 1:m, z = A * z;, endfor
7b190a2f11cb maint: Use 2-space indent in definition of BIST %!functions.
Rik <rik@octave.org>
parents: 22330
diff changeset
276 %! endswitch
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
277 %!endfunction
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
278
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
279 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
280 %! A = reshape ((1:16)-8, 4, 4);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
281 %! assert (normest1 (A), norm (A, 1), eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
282
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
283 ## test t=1
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
284 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
285 %! A = rand (4); # for positive matrices always work
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
286 %! assert (normest1 (A, 1), norm (A,1), 2 * eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
287
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
288 ## test t=3
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
289 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
290 %! A = [-0.21825 0.16598 0.19388 0.75297
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
291 %! -1.47732 0.78443 -1.04254 0.42240
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
292 %! 1.39857 -0.34046 2.28617 0.68089
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
293 %! 0.31205 1.50529 -0.75804 -1.22476];
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
294 %! X = [1,1,-1;1,1,1;1,1,-1;1,-1,-1]/3;
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
295 %! assert (normest1 (A, 3, X), norm (A, 1), 2 * eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
296
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
297 ## test Afun
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
298 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
299 %! A = rand (10);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
300 %! n = length (A);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
301 %! Afun = @(flag, x) afun_A (flag, x, A, n);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
302 %! assert (normest1 (Afun), norm (A, 1), 2 * eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
303
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
304 ## test Afun with parameters
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
305 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
306 %! A = rand (10);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
307 %! assert (normest1 (@afun_A_P, [], [], A, 3), norm (A ^ 3, 1), 1000 * eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
308
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
309 ## test output
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
310 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
311 %! A = reshape (1:16,4,4);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
312 %! [nest, v, w, iter] = normest1 (A);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
313 %! assert (norm (w, 1), nest * norm (v, 1), eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
314
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
315 ## test output
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
316 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
317 %! A = rand (100);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
318 %! A(A <= 1/3) = -1;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
319 %! A(A > 1/3 & A <= 2/3) = 0;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
320 %! A(A > 2/3) = 1;
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
321 %! [nest, v, w, iter] = normest1 (A, 10);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
322 %! assert (w, A * v, eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
323
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
324 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
325 %! A = rand (5);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
326 %! nest = normest1 (A, 6);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
327 %! assert (nest, norm (A,1), eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
328
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
329 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
330 %! A = rand (5);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
331 %! nest = normest1 (A, 2, ones (5, 2) / 5);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
332 %! assert (nest, norm (A,1), eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
333
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
334 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
335 %! N = 10;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
336 %! A = ones (N);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
337 %! [nm1, v1, w1] = normest1 (A);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
338 %! [nminf, vinf, winf] = normest1 (A', 6);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
339 %! assert (nm1, N, -2*eps);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
340 %! assert (nminf, N, -2*eps);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
341 %! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
342 %! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
343
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
344 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
345 %! N = 5;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
346 %! A = hilb (N);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
347 %! [nm1, v1, w1] = normest1 (A);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
348 %! [nminf, vinf, winf] = normest1 (A', 6);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
349 %! assert (nm1, norm (A, 1), -2*eps);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
350 %! assert (nminf, norm (A, inf), -2*eps);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
351 %! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
352 %! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
353
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
354 ## Only likely to be within a factor of 10.
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
355 %!test
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
356 %! old_state = rand ("state");
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
357 %! unwind_protect
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
358 %! rand ("state", 42); # Initialize to guarantee reproducible results
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
359 %! N = 100;
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
360 %! A = rand (N);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
361 %! [nm1, v1, w1] = normest1 (A);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
362 %! [nminf, vinf, winf] = normest1 (A', 6);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
363 %! assert (nm1, norm (A, 1), -.1);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
364 %! assert (nminf, norm (A, inf), -.1);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
365 %! assert (norm (w1, 1), nm1 * norm (v1, 1), -2*eps);
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
366 %! assert (norm (winf, 1), nminf * norm (vinf, 1), -2*eps);
22183
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
367 %! unwind_protect_cleanup
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
368 %! rand ("state", old_state);
bfb1b089c230 New function normest1 as replacement for onenormest (patch #8837)
Marco Caliari <marco.caliari@univr.it>
parents:
diff changeset
369 %! end_unwind_protect
22193
e35866e6a2e0 normest1: always return a column vector for IT output (patch #8837)
Carnë Draug <carandraug@octave.org>
parents: 22183
diff changeset
370
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
371 ## Check ITER is always a column vector.
22193
e35866e6a2e0 normest1: always return a column vector for IT output (patch #8837)
Carnë Draug <carandraug@octave.org>
parents: 22183
diff changeset
372 %!test
e35866e6a2e0 normest1: always return a column vector for IT output (patch #8837)
Carnë Draug <carandraug@octave.org>
parents: 22183
diff changeset
373 %! [~, ~, ~, it] = normest1 (rand (3), 3);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
374 %! assert (iscolumn (it));
22193
e35866e6a2e0 normest1: always return a column vector for IT output (patch #8837)
Carnë Draug <carandraug@octave.org>
parents: 22183
diff changeset
375 %! [~, ~, ~, it] = normest1 (rand (50), 20);
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
376 %! assert (iscolumn (it));
22765
01aae08a0105 maint: Rename variables to match documentation in m-files.
Rik <rik@octave.org>
parents: 22755
diff changeset
377
22829
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
378 ## Test input validation
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
379 %!error normest1 ()
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
380 %!error <A must be a square matrix or a function handle> normest1 ({1})
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
381 %!error <A must be a square matrix> normest1 ([1 2])
bab05fcd38cb normest1.m: Overhaul function.
Rik <rik@octave.org>
parents: 22828
diff changeset
382 %!error <X0 must have 2 columns> normest1 (magic (5), :, [1])