Mercurial > octave
annotate scripts/special-matrix/hilb.m @ 28520:627da618dcc4 stable
hilb.m: Speed-up function using broadcasting.
* scripts/special-matrix/hilb.m: Speed-up function using broadcasting. For more
information read https://octave.discourse.group/t/boosting-hilb/48 .
author | Kai T. Ohlhus <k.ohlhus@gmail.com> |
---|---|
date | Wed, 01 Jul 2020 16:03:41 +0900 |
parents | bd51beb6205e |
children | d8318c12d903 0a5b15007766 |
rev | line source |
---|---|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 ######################################################################## |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 ## |
27919
1891570abac8
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
27918
diff
changeset
|
3 ## Copyright (C) 1993-2020 The Octave Project Developers |
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
27898
diff
changeset
|
4 ## |
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 ## See the file COPYRIGHT.md in the top-level directory of this |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 ## distribution or <https://octave.org/copyright/>. |
2313 | 7 ## |
8 ## This file is part of Octave. | |
9 ## | |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
10 ## Octave is free software: you can redistribute it and/or modify it |
2313 | 11 ## 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:
23220
diff
changeset
|
12 ## the Free Software Foundation, either version 3 of the License, or |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
13 ## (at your option) any later version. |
2313 | 14 ## |
15 ## Octave is distributed in the hope that it will be useful, but | |
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
18 ## GNU General Public License for more details. |
2313 | 19 ## |
20 ## You should have received a copy of the GNU General Public License | |
7016 | 21 ## 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:
23220
diff
changeset
|
22 ## <https://www.gnu.org/licenses/>. |
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 ## |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 ######################################################################## |
245 | 25 |
3369 | 26 ## -*- texinfo -*- |
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20735
diff
changeset
|
27 ## @deftypefn {} {} hilb (@var{n}) |
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
28 ## Return the Hilbert matrix of order @var{n}. |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
29 ## |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
30 ## The @math{i,j} element of a Hilbert matrix is defined as |
3369 | 31 ## @tex |
32 ## $$ | |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
33 ## H(i, j) = {1 \over (i + j - 1)} |
3369 | 34 ## $$ |
35 ## @end tex | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7411
diff
changeset
|
36 ## @ifnottex |
3426 | 37 ## |
3369 | 38 ## @example |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
39 ## H(i, j) = 1 / (i + j - 1) |
3369 | 40 ## @end example |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
9211
diff
changeset
|
41 ## |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7411
diff
changeset
|
42 ## @end ifnottex |
12639
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
43 ## |
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
44 ## Hilbert matrices are close to being singular which make them difficult to |
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
45 ## invert with numerical routines. Comparing the condition number of a random |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
46 ## matrix 5x5 matrix with that of a Hilbert matrix of order 5 reveals just how |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
47 ## difficult the problem is. |
12639
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
48 ## |
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
49 ## @example |
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
50 ## @group |
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
51 ## cond (rand (5)) |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
52 ## @result{} 14.392 |
12639
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
53 ## cond (hilb (5)) |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
54 ## @result{} 4.7661e+05 |
12639
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
55 ## @end group |
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
56 ## @end example |
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
57 ## |
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
58 ## @seealso{invhilb} |
3369 | 59 ## @end deftypefn |
4 | 60 |
2311 | 61 function retval = hilb (n) |
4 | 62 |
63 if (nargin != 1) | |
6046 | 64 print_usage (); |
13872
779e15b69738
hilb.m: 10% speedup by using in-place accumulation.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
65 elseif (! isscalar (n)) |
779e15b69738
hilb.m: 10% speedup by using in-place accumulation.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
66 error ("hilb: N must be a scalar integer"); |
4 | 67 endif |
68 | |
28520
627da618dcc4
hilb.m: Speed-up function using broadcasting.
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
27923
diff
changeset
|
69 ## Very elegant solution by N. Higham |
627da618dcc4
hilb.m: Speed-up function using broadcasting.
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
27923
diff
changeset
|
70 ## https://nhigham.com/2020/06/30/what-is-the-hilbert-matrix/ |
627da618dcc4
hilb.m: Speed-up function using broadcasting.
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
27923
diff
changeset
|
71 j = 1:n; |
627da618dcc4
hilb.m: Speed-up function using broadcasting.
Kai T. Ohlhus <k.ohlhus@gmail.com>
parents:
27923
diff
changeset
|
72 retval = 1 ./ (j' + j - 1); |
4 | 73 |
74 endfunction | |
7411 | 75 |
13872
779e15b69738
hilb.m: 10% speedup by using in-place accumulation.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
76 |
779e15b69738
hilb.m: 10% speedup by using in-place accumulation.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
77 %!assert (hilb (2), [1, 1/2; 1/2, 1/3]) |
779e15b69738
hilb.m: 10% speedup by using in-place accumulation.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
78 %!assert (hilb (3), [1, 1/2, 1/3; 1/2, 1/3, 1/4; 1/3, 1/4, 1/5]) |
7411 | 79 |
13872
779e15b69738
hilb.m: 10% speedup by using in-place accumulation.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
80 %!error hilb () |
779e15b69738
hilb.m: 10% speedup by using in-place accumulation.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
81 %!error hilb (1, 2) |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
82 %!error <N must be a scalar integer> hilb (ones (2)) |