Mercurial > octave-antonio
annotate scripts/special-matrix/invhilb.m @ 20162:2645f9ef8c88 stable
doc: Update more docstrings to have one sentence summary as first line.
Reviewed specfun, special-matrix, testfun, and time script directories.
* scripts/specfun/expint.m, scripts/specfun/isprime.m,
scripts/specfun/legendre.m, scripts/specfun/primes.m,
scripts/specfun/reallog.m, scripts/specfun/realsqrt.m,
scripts/special-matrix/gallery.m, scripts/special-matrix/hadamard.m,
scripts/special-matrix/hankel.m, scripts/special-matrix/hilb.m,
scripts/special-matrix/invhilb.m, scripts/special-matrix/magic.m,
scripts/special-matrix/pascal.m, scripts/special-matrix/rosser.m,
scripts/special-matrix/toeplitz.m, scripts/special-matrix/vander.m,
scripts/special-matrix/wilkinson.m, scripts/testfun/assert.m,
scripts/testfun/demo.m, scripts/testfun/example.m, scripts/testfun/fail.m,
scripts/testfun/rundemos.m, scripts/testfun/runtests.m,
scripts/testfun/speed.m, scripts/time/asctime.m, scripts/time/calendar.m,
scripts/time/clock.m, scripts/time/ctime.m, scripts/time/datenum.m,
scripts/time/datestr.m, scripts/time/datevec.m, scripts/time/etime.m,
scripts/time/is_leap_year.m, scripts/time/now.m, scripts/time/weekday.m:
Update more docstrings to have one sentence summary as first line.
author | Rik <rik@octave.org> |
---|---|
date | Sun, 03 May 2015 17:00:11 -0700 |
parents | 4197fc428c7d |
children |
rev | line source |
---|---|
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19040
diff
changeset
|
1 ## Copyright (C) 1993-2015 Dirk Laurie |
2313 | 2 ## |
3 ## This file is part of Octave. | |
4 ## | |
5 ## Octave is free software; you can redistribute it and/or modify it | |
6 ## under the terms of the GNU General Public License as published by | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
2313 | 9 ## |
10 ## Octave is distributed in the hope that it will be useful, but | |
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 ## General Public License for more details. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
245 | 18 |
3369 | 19 ## -*- texinfo -*- |
20 ## @deftypefn {Function File} {} invhilb (@var{n}) | |
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
21 ## Return the inverse of 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
|
22 ## |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
23 ## This can be computed exactly using |
3889 | 24 ## @tex |
25 ## $$\eqalign{ | |
26 ## A_{ij} &= -1^{i+j} (i+j-1) | |
27 ## \left( \matrix{n+i-1 \cr n-j } \right) | |
28 ## \left( \matrix{n+j-1 \cr n-i } \right) | |
29 ## \left( \matrix{i+j-2 \cr i-2 } \right)^2 \cr | |
30 ## &= { p(i)p(j) \over (i+j-1) } | |
31 ## }$$ | |
32 ## where | |
33 ## $$ | |
34 ## p(k) = -1^k \left( \matrix{ k+n-1 \cr k-1 } \right) | |
35 ## \left( \matrix{ n \cr k } \right) | |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
36 ## $$ |
3889 | 37 ## @end tex |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7411
diff
changeset
|
38 ## @ifnottex |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
39 ## |
3889 | 40 ## @example |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
41 ## @group |
3889 | 42 ## |
43 ## (i+j) /n+i-1\ /n+j-1\ /i+j-2\ 2 | |
44 ## A(i,j) = -1 (i+j-1)( )( ) ( ) | |
45 ## \ n-j / \ n-i / \ i-2 / | |
46 ## | |
47 ## = p(i) p(j) / (i+j-1) | |
48 ## | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
49 ## @end group |
3889 | 50 ## @end example |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
51 ## |
10846
a4f482e66b65
Grammarcheck more of the documentation.
Rik <octave@nomad.inbox5.com>
parents:
10821
diff
changeset
|
52 ## @noindent |
3889 | 53 ## where |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
54 ## |
3889 | 55 ## @example |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
56 ## @group |
3889 | 57 ## k /k+n-1\ /n\ |
58 ## p(k) = -1 ( ) ( ) | |
59 ## \ k-1 / \k/ | |
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
60 ## @end group |
3889 | 61 ## @end example |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10791
diff
changeset
|
62 ## |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7411
diff
changeset
|
63 ## @end ifnottex |
20162
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
64 ## The validity of this formula can easily be checked by expanding the binomial |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
65 ## coefficients in both formulas as factorials. It can be derived more |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
66 ## directly via the theory of Cauchy matrices. See @nospell{J. W. Demmel}, |
2645f9ef8c88
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
67 ## @cite{Applied Numerical Linear Algebra}, p. 92. |
3889 | 68 ## |
69 ## Compare this with the numerical calculation of @code{inverse (hilb (n))}, | |
3369 | 70 ## which suffers from the ill-conditioning of the Hilbert matrix, and the |
71 ## finite precision of your computer's floating point arithmetic. | |
12639
4d777e05d47c
doc: Review and update documentation for "Matrix Manipulation" chapter.
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
72 ## @seealso{hilb} |
3369 | 73 ## @end deftypefn |
4 | 74 |
5132 | 75 ## Author: Dirk Laurie <dlaurie@na-net.ornl.gov> |
2314 | 76 |
2311 | 77 function retval = invhilb (n) |
4 | 78 |
79 if (nargin != 1) | |
6046 | 80 print_usage (); |
13892
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
81 elseif (! isscalar (n)) |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
82 error ("invhilb: N must be a scalar integer"); |
4 | 83 endif |
84 | |
13892
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
85 ## The point about the second formula above is that when vectorized, |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
86 ## p(k) is evaluated for k=1:n which involves O(n) calls to bincoeff |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
87 ## instead of O(n^2). |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
88 ## |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
89 ## We evaluate the expression as (-1)^(i+j)*(p(i)*p(j))/(i+j-1) except |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
90 ## when p(i)*p(j) would overflow. In cases where p(i)*p(j) is an exact |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
91 ## machine number, the result is also exact. Otherwise we calculate |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
92 ## (-1)^(i+j)*p(i)*(p(j)/(i+j-1)). |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
93 ## |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
94 ## The Octave bincoeff routine uses transcendental functions (gammaln |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
95 ## and exp) rather than multiplications, for the sake of speed. |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
96 ## However, it rounds the answer to the nearest integer, which |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
97 ## justifies the claim about exactness made above. |
3889 | 98 |
13892
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
99 retval = zeros (n); |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
100 k = [1:n]; |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
101 p = k .* bincoeff (k+n-1, k-1) .* bincoeff (n, k); |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
102 p(2:2:n) = -p(2:2:n); |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
103 if (n < 203) |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
104 for l = 1:n |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
105 retval(l,:) = (p(l) * p) ./ [l:l+n-1]; |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
106 endfor |
4 | 107 else |
13892
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
108 for l = 1:n |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
109 retval(l,:) = p(l) * (p ./ [l:l+n-1]); |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
110 endfor |
4 | 111 endif |
112 | |
113 endfunction | |
7411 | 114 |
13892
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
115 |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
116 %!assert (invhilb (1), 1) |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
117 %!assert (invhilb (2), [4, -6; -6, 12]) |
7411 | 118 %!test |
13892
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
119 %! result4 = [16 , -120 , 240 , -140; |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
120 %! -120, 1200 , -2700, 1680; |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
121 %! 240 , -2700, 6480 , -4200; |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
122 %! -140, 1680 , -4200, 2800]; |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
123 %! assert (invhilb (4), result4); |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
124 %!assert (invhilb (7) * hilb (7), eye (7), sqrt (eps)) |
7411 | 125 |
13892
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
126 %!error invhilb () |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
127 %!error invhilb (1, 2) |
dd01f0bfd78d
invhilb.m: update coding style.
Rik <octave@nomad.inbox5.com>
parents:
12639
diff
changeset
|
128 %!error <N must be a scalar integer> invhilb ([1, 2]) |
7411 | 129 |