Mercurial > octave
annotate scripts/linear-algebra/logm.m @ 29770:4778b21b1386 stable
expm.m, logm.m: Use function "isdiag" to detect if input is a diagonal matrix (bug #60738).
* scripts/linear-algebra/expm.m, scripts/linear-algebra/logm.m: Use function
"isdiag" to detect if input is a diagonal matrix. Add tests.
* scripts/linear-algebra/isdiag.m: Add test.
author | Markus Mützel <markus.muetzel@gmx.de> |
---|---|
date | Thu, 17 Jun 2021 18:23:26 +0200 |
parents | fe06c183cc4e |
children | 42dc5cf93f83 20380c9bed30 |
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 ## |
29358
0a5b15007766
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
27985
diff
changeset
|
3 ## Copyright (C) 2008-2021 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
|
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/>. |
4334 | 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 |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
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. |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
14 ## |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
15 ## 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
|
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. |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
19 ## |
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
20 ## 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
|
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 ######################################################################## |
4334 | 25 |
26 ## -*- texinfo -*- | |
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
27 ## @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
|
28 ## @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
|
29 ## @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
|
30 ## 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
|
31 ## |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
32 ## 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
|
33 ## |
10827
228cd18455a6
logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents:
10825
diff
changeset
|
34 ## @example |
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
35 ## 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
|
36 ## @end example |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
37 ## |
20160
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
38 ## 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
|
39 ## 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
|
40 ## |
03b9d17a2d95
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
41 ## 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
|
42 ## computed. |
12584
7ef7e20057fa
Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
43 ## @seealso{expm, sqrtm} |
4334 | 44 ## @end deftypefn |
45 | |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
46 ## 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
|
47 ## (SIAM, 2008.) |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
48 ## |
4334 | 49 |
16768 | 50 ## Author: N. J. Higham |
51 ## Author: Richard T. Guy <guyrt7@wfu.edu> | |
27985
9f9ac219896d
maint: Remove remaining "Author:" instances from code base.
Rik <rik@octave.org>
parents:
27978
diff
changeset
|
52 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
53 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
|
54 |
11458 | 55 if (nargin == 0 || nargin > 2) |
6046 | 56 print_usage (); |
4334 | 57 endif |
58 | |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
59 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
|
60 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
|
61 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
62 |
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
|
63 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
|
64 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
|
65 return; |
29770
4778b21b1386
expm.m, logm.m: Use function "isdiag" to detect if input is a diagonal matrix (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29769
diff
changeset
|
66 elseif (isdiag (A)) |
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
|
67 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
|
68 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
|
69 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
|
70 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
71 [u, s] = schur (A); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
72 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
73 if (isreal (A)) |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
74 [u, s] = rsf2csf (u, s); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
75 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
76 |
13097
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
77 eigv = diag (s); |
29769
fe06c183cc4e
logm.m: Allow tolerance in check for real negative values in complex vector (bug #60738).
Steven <steven.waldrip@gmail.com>
parents:
29738
diff
changeset
|
78 tol = rows (A) * eps (max (abs (eigv))); |
fe06c183cc4e
logm.m: Allow tolerance in check for real negative values in complex vector (bug #60738).
Steven <steven.waldrip@gmail.com>
parents:
29738
diff
changeset
|
79 real_neg_eigv = (real (eigv) < -tol) & (imag (eigv) <= tol); |
29738
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
80 if (any (real_neg_eigv)) |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
81 warning ("Octave:logm:non-principal", |
11458 | 82 "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
|
83 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
84 |
29738
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
85 real_eig = ! any (real_neg_eigv); |
13097
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
86 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
87 k = 0; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
88 ## 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
|
89 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
|
90 p = 0; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
91 m = 7; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
92 while (k < opt_iters) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
93 tau = norm (s - eye (size (s)),1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
94 if (tau <= theta (7)) |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
95 p += 1; |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
96 j(1) = find (tau <= theta, 1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
97 j(2) = find (tau / 2 <= theta, 1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
98 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
|
99 m = j(1); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
100 break |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
101 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
102 endif |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
103 k += 1; |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
104 s = sqrtm (s); |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
105 endwhile |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
106 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
107 if (k >= opt_iters) |
11458 | 108 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
|
109 endif |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
110 |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
111 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
|
112 |
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
|
113 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
|
114 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
|
115 endif |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
116 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
117 s = 2^k * u * s * u'; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
118 |
13097
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
119 ## 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
|
120 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
|
121 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
|
122 endif |
52e4aa30d5b2
logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents:
12584
diff
changeset
|
123 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
124 if (nargout == 2) |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
125 iters = k; |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
126 endif |
4334 | 127 |
128 endfunction | |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
129 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
130 ################## ANCILLARY FUNCTIONS ################################ |
11033
d9c8916bb9dd
Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents:
10827
diff
changeset
|
131 ###### 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
|
132 ###### Reference: |
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
133 ###### 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
|
134 ###### (SIAM, 2008.). |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
135 ####################################################################### |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
136 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
137 ##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
|
138 ## 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
|
139 ## 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
|
140 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
141 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
|
142 |
11458 | 143 [nodes, wts] = gauss_legendre (m); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
144 ## Convert from [-1,1] to [0,1]. |
11458 | 145 nodes = (nodes+1)/2; |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20160
diff
changeset
|
146 wts /= 2; |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
147 |
11471
994e2a93a8e2
Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents:
11465
diff
changeset
|
148 n = length (A); |
11458 | 149 s = zeros (n); |
150 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
|
151 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
|
152 endfor |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
153 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
154 endfunction |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
155 |
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 ## 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
|
158 ## [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
|
159 ## for N-point Gauss-Legendre quadrature. |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
160 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
161 ## Reference: |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
162 ## 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
|
163 ## 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
|
164 |
11458 | 165 function [x, w] = gauss_legendre (n) |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
166 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
167 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
|
168 v = i./sqrt ((2*i).^2-1); |
11458 | 169 [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
|
170 x = diag (D); |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
171 w = 2*(V(1,:)'.^2); |
21758
ffad2baa90f7
maint: Use newlines to make code more readable.
Rik <rik@octave.org>
parents:
20852
diff
changeset
|
172 |
10825
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
173 endfunction |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
174 |
cace99cb01ab
rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents:
7017
diff
changeset
|
175 |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
176 %!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
|
177 %!test |
e5facc6eec13
Silence warning messages in %!test code.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
178 %! warning ("off", "Octave:logm:non-principal", "local"); |
e5facc6eec13
Silence warning messages in %!test code.
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
179 %! 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
|
180 %!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
|
181 %!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
|
182 %!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
|
183 %!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
|
184 %!assert (logm (expm ([0 1i; -1i 0])), [0 1i; -1i 0], 10 * eps) |
29738
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
185 %!test <*60738> |
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
186 %! A = [0.2510, 1.2808, -1.2252; ... |
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
187 %! 0.2015, 1.0766, 0.5630; ... |
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
188 %! -1.9769, -1.0922, -0.5831]; |
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
189 %! warning ("off", "Octave:logm:non-principal", "local"); |
5928dd4f7a02
logm.m: Fix check for real negative values in complex vector (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29358
diff
changeset
|
190 %! assert (expm (logm (A)), A, 40*eps); |
29770
4778b21b1386
expm.m, logm.m: Use function "isdiag" to detect if input is a diagonal matrix (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29769
diff
changeset
|
191 %!assert (expm (logm (diag (ones (1, 3)))), diag (ones (1, 3))); |
4778b21b1386
expm.m, logm.m: Use function "isdiag" to detect if input is a diagonal matrix (bug #60738).
Markus Mützel <markus.muetzel@gmx.de>
parents:
29769
diff
changeset
|
192 %!assert (expm (logm (zeros (3))), zeros (3)); |
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
|
193 |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
194 ## Test input validation |
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
195 %!error logm () |
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
196 %!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
|
197 %!error <logm: A must be a square matrix> logm ([1 0;0 1; 2 2]) |