Mercurial > octave
annotate scripts/linear-algebra/logm.m @ 27919:1891570abac8
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 |
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) 2008-2020 The Octave Project Developers |
27918
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
2 ## |
b442ec6dda5c
use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents:
26376
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:
26376
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:
26376
diff
changeset
|
5 ## |
4334 | 6 ## |
7 ## This file is part of Octave. | |
8 ## | |
24534
194eb4bd202b
maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
9 ## Octave is free software: you can redistribute it and/or modify it |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
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:
23220
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:
22323
diff
changeset
|
12 ## (at your option) any later version. |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
13 ## |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
14 ## Octave is distributed in the hope that it will be useful, but |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
15 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
22755
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
16 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
3a2b891d0b33
maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents:
22323
diff
changeset
|
17 ## GNU General Public License for more details. |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
18 ## |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
19 ## You should have received a copy of the GNU General Public License |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
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:
23220
diff
changeset
|
21 ## <https://www.gnu.org/licenses/>. |
4334 | 22 |
23 ## -*- texinfo -*- | |
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
24 ## @deftypefn {} {@var{s} =} logm (@var{A}) |
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
25 ## @deftypefnx {} {@var{s} =} logm (@var{A}, @var{opt_iters}) |
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
26 ## @deftypefnx {} {[@var{s}, @var{iters}] =} logm (@dots{}) |
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
27 ## Compute the matrix logarithm of the square matrix @var{A}. |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
28 ## |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
29 ## The implementation utilizes a Pad@'e approximant and the identity |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
30 ## |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
31 ## @example |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
32 ## logm (@var{A}) = 2^k * logm (@var{A}^(1 / 2^k)) |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
33 ## @end example |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
34 ## |
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
35 ## The optional input @var{opt_iters} is the maximum number of square roots |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
36 ## to compute and defaults to 100. |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
37 ## |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
38 ## The optional output @var{iters} is the number of square roots actually |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
39 ## computed. |
12584
7ef7e20057fa
Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
40 ## @seealso{expm, sqrtm} |
4334 | 41 ## @end deftypefn |
42 | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
43 ## Reference: N. J. Higham, Functions of Matrices: Theory and Computation |
11033
d9c8916bb9dd
Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents:
10827
diff
changeset
|
44 ## (SIAM, 2008.) |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
45 ## |
4334 | 46 |
16768 | 47 ## Author: N. J. Higham |
48 ## Author: Richard T. Guy <guyrt7@wfu.edu> | |
49 ## Author: Marco Caliari <marco.caliari@univr.it> | |
50 | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
51 function [s, iters] = logm (A, opt_iters = 100) |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
52 |
11458 | 53 if (nargin == 0 || nargin > 2) |
6046 | 54 print_usage (); |
4334 | 55 endif |
56 | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
57 if (! issquare (A)) |
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11471
diff
changeset
|
58 error ("logm: A must be a square matrix"); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
59 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
60 |
11476
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
61 if (isscalar (A)) |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
62 s = log (A); |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
63 return; |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
64 elseif (strfind (typeinfo (A), "diagonal matrix")) |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
65 s = diag (log (diag (A))); |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
66 return; |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
67 endif |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
68 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
69 [u, s] = schur (A); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
70 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
71 if (isreal (A)) |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
72 [u, s] = rsf2csf (u, s); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
73 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
74 |
13097
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
75 eigv = diag (s); |
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
76 if (any (eigv < 0)) |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
77 warning ("Octave:logm:non-principal", |
11458 | 78 "logm: principal matrix logarithm is not defined for matrices with negative eigenvalues; computing non-principal logarithm"); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
79 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
80 |
13097
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
81 real_eig = all (eigv >= 0); |
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
82 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
83 k = 0; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
84 ## Algorithm 11.9 in "Function of matrices", by N. Higham |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
85 theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1]; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
86 p = 0; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
87 m = 7; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
88 while (k < opt_iters) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
89 tau = norm (s - eye (size (s)),1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
90 if (tau <= theta (7)) |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
91 p += 1; |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
92 j(1) = find (tau <= theta, 1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
93 j(2) = find (tau / 2 <= theta, 1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
94 if (j(1) - j(2) <= 1 || p == 2) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
95 m = j(1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
96 break |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
97 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
98 endif |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
99 k += 1; |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
100 s = sqrtm (s); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
101 endwhile |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
102 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
103 if (k >= opt_iters) |
11458 | 104 warning ("logm: maximum number of square roots exceeded; results may still be accurate"); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
105 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
106 |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
107 s -= eye (size (s)); |
11476
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
108 |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
109 if (m > 1) |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
110 s = logm_pade_pf (s, m); |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
111 endif |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
112 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
113 s = 2^k * u * s * u'; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
114 |
13097
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
115 ## Remove small complex values (O(eps)) which may have entered calculation |
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14491
diff
changeset
|
116 if (real_eig && isreal (A)) |
13097
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
117 s = real (s); |
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
118 endif |
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
119 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
120 if (nargout == 2) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
121 iters = k; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
122 endif |
4334 | 123 |
124 endfunction | |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
125 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
126 ################## ANCILLARY FUNCTIONS ################################ |
11033
d9c8916bb9dd
Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents:
10827
diff
changeset
|
127 ###### Taken from the mfttoolbox (GPL 3) by D. Higham. |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
128 ###### Reference: |
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
129 ###### D. Higham, Functions of Matrices: Theory and Computation |
11033
d9c8916bb9dd
Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents:
10827
diff
changeset
|
130 ###### (SIAM, 2008.). |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
131 ####################################################################### |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
132 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
133 ##LOGM_PADE_PF Evaluate Pade approximant to matrix log by partial fractions. |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
134 ## Y = LOGM_PADE_PF(A,M) evaluates the [M/M] Pade approximation to |
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
135 ## LOG(EYE(SIZE(A))+A) using a partial fraction expansion. |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
136 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
137 function s = logm_pade_pf (A, m) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
138 |
11458 | 139 [nodes, wts] = gauss_legendre (m); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
140 ## Convert from [-1,1] to [0,1]. |
11458 | 141 nodes = (nodes+1)/2; |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
142 wts /= 2; |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
143 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
144 n = length (A); |
11458 | 145 s = zeros (n); |
146 for j = 1:m | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
147 s += wts(j)*(A/(eye (n) + nodes(j)*A)); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
148 endfor |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
149 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
150 endfunction |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
151 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
152 ###################################################################### |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
153 ## GAUSS_LEGENDRE Nodes and weights for Gauss-Legendre quadrature. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
154 ## [X,W] = GAUSS_LEGENDRE(N) computes the nodes X and weights W |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
155 ## for N-point Gauss-Legendre quadrature. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
156 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
157 ## Reference: |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
158 ## G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
159 ## rules, Math. Comp., 23(106):221-230, 1969. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
160 |
11458 | 161 function [x, w] = gauss_legendre (n) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
162 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
163 i = 1:n-1; |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
164 v = i./sqrt ((2*i).^2-1); |
11458 | 165 [V, D] = eig (diag (v, -1) + diag (v, 1)); |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
166 x = diag (D); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
167 w = 2*(V(1,:)'.^2); |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
168 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
169 endfunction |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
170 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
171 |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
172 %!assert (norm (logm ([1 -1;0 1]) - [0 -1; 0 0]) < 1e-5) |
19670
e5facc6eec13
Silence warning messages in %!test code.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
173 %!test |
e5facc6eec13
Silence warning messages in %!test code.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
174 %! warning ("off", "Octave:logm:non-principal", "local"); |
e5facc6eec13
Silence warning messages in %!test code.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
175 %! assert (norm (expm (logm ([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5); |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
176 %!assert (logm ([1 -1 -1;0 1 -1; 0 0 1]), [0 -1 -1.5; 0 0 -1; 0 0 0], 1e-5) |
11476
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
177 %!assert (logm (10), log (10)) |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
178 %!assert (full (logm (eye (3))), logm (full (eye (3)))) |
ff7e0776ba0f
logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents:
11472
diff
changeset
|
179 %!assert (full (logm (10*eye (3))), logm (full (10*eye (3))), 8*eps) |
14481
d2bffa78730e
Fix logm for complex matrix with real eigenvalues (bug #34893).
Marco Caliari <marco.caliari@univr.it>
parents:
14327
diff
changeset
|
180 %!assert (logm (expm ([0 1i; -1i 0])), [0 1i; -1i 0], 10 * eps) |
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11471
diff
changeset
|
181 |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
182 ## Test input validation |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
183 %!error logm () |
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
184 %!error logm (1, 2, 3) |
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
185 %!error <logm: A must be a square matrix> logm ([1 0;0 1; 2 2]) |