octave: scripts/linear-algebra/logm.m annotate

annotate scripts/linear-algebra/logm.m @ 31551:fd29c7a50a78 stable

maint: use commas, semicolons consistently with Octave conventions. * makeValidName.m: Remove %!test and move BIST %!asserts to column 1. * base64decode.m, base64encode.m, which.m, logm.m, uniquetol.m, perms.m: Delete semicolon (';') at end of %!assert BIST. * lin2mu.m, interp2.m, interpn.m, lsqnonneg.m, pqpnonneg.m, uniquetol.m, betainc.m, normalize.m: Add semicolon (';') to end of assert statement within %!test BIST. * __memoize__.m, tar_is_bsd.m, publish.m: Add semicolon (';') to line with keyword "persistent". * stft.m: Use comma (',') after "case" keyword when code immediately follows. * gallery.m: Align commas used in case statements in massive switch block. Remove unnecessary parentheses around a numeric case argument. * ranks.m: Remove semicolon (';') from case statemnt argument.

author	Rik <rik@octave.org>
date	Sat, 26 Nov 2022 06:32:08 -0800
parents	a40c0b7aa376
children	597f3ee61a48

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 ##
30564 796f54d4ddbf update Octave Project Developers copyright for the new year John W. Eaton <jwe@octave.org> parents: 30383 diff changeset	3 ## Copyright (C) 2008-2022 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 764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	7 ##
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	8 ## This file is part of Octave.
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	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 764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	25
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	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 764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	44 ## @end deftypefn
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	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 764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	49
16768 e2de3c8882be copyright notice fixes John W. Eaton <jwe@octave.org> parents: 14868 diff changeset	50 ## Author: N. J. Higham
e2de3c8882be copyright notice fixes John W. Eaton <jwe@octave.org> parents: 14868 diff changeset	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
28789 28de41192f3c Eliminate unneeded verification of nargin, nargout in m-files. Rik <rik@octave.org> parents: 27985 diff changeset	55 if (nargin == 0)
6046 34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe] jwe parents: 5307 diff changeset	56 print_usage ();
4334 764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	57 endif
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	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);
30383 8afdeac24ba4 maint: Use space after function name and before opening parenthesis. Rik <rik@octave.org> parents: 29806 diff changeset	78 n = rows (A);
29778 c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	79 tol = n * eps (max (abs (eigv)));
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	80 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	81 if (any (real_neg_eigv))
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	82 warning ("Octave:logm:non-principal",
11458 93a039fe681e logm: style fixes John W. Eaton <jwe@octave.org> parents: 11033 diff changeset	83 "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	84 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	85
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	86 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	87
29778 c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	88 if (max (abs (triu (s,1))(:)) < tol)
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	89 ## Will run for Hermitian matrices as Schur decomposition is diagonal.
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	90 ## This way is faster and more accurate but only works on a diagonal matrix.
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	91 logeigv = log (eigv);
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	92 logeigv(isinf (logeigv)) = -log (realmax ());
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	93 s = u * diag (logeigv) * u';
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	94 iters = 0;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	95 else
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	96 k = 0;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	97 ## Algorithm 11.9 in "Function of matrices", by N. Higham
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	98 theta = [0, 0, 1.61e-2, 5.38e-2, 1.13e-1, 1.86e-1, 2.6429608311114350e-1];
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	99 p = 0;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	100 m = 7;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	101 while (k < opt_iters)
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	102 tau = norm (s - eye (n), 1);
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	103 if (tau <= theta (7))
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	104 p += 1;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	105 j(1) = find (tau <= theta, 1);
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	106 j(2) = find (tau / 2 <= theta, 1);
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	107 if (j(1) - j(2) <= 1 \|\| p == 2)
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	108 m = j(1);
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	109 break;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	110 endif
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	111 endif
29778 c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	112 k += 1;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	113 s = sqrtm (s);
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	114 endwhile
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	115
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	116 if (k >= opt_iters)
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	117 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	118 endif
29778 c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	119
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	120 s -= eye (n);
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	121
29778 c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	122 if (m > 1)
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	123 s = logm_pade_pf (s, m);
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	124 endif
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	125
29778 c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	126 s = 2^k * u * s * u';
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	127
29778 c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	128 if (nargout == 2)
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	129 iters = k;
c44e54bb018b logm.m: Improve performance and accuracy for Hermitian input (bug #60738). Steven Waldrip <steven.waldrip@gmail.com> parents: 29771 diff changeset	130 endif
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	131 endif
13097 52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121). Rik <octave@nomad.inbox5.com> parents: 12584 diff changeset	132 ## 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	133 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	134 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	135 endif
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121). Rik <octave@nomad.inbox5.com> parents: 12584 diff changeset	136
4334 764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe] jwe parents: diff changeset	137 endfunction
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	138
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	139 ################## ANCILLARY FUNCTIONS ################################
11033 d9c8916bb9dd Untabify a few remaining .m scripts. Rik <octave@nomad.inbox5.com> parents: 10827 diff changeset	140 ###### 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	141 ###### Reference:
c792872f8942 all script files: untabify and strip trailing whitespace John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	142 ###### 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	143 ###### (SIAM, 2008.).
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	144 #######################################################################
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	145
31253 a40c0b7aa376 maint: changes to follow Octave coding conventions. Rik <rik@octave.org> parents: 30564 diff changeset	146 ## 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	147 ## 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	148 ## 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	149
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 11465 diff changeset	150 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	151
11458 93a039fe681e logm: style fixes John W. Eaton <jwe@octave.org> parents: 11033 diff changeset	152 [nodes, wts] = gauss_legendre (m);
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	153 ## Convert from [-1,1] to [0,1].
11458 93a039fe681e logm: style fixes John W. Eaton <jwe@octave.org> parents: 11033 diff changeset	154 nodes = (nodes+1)/2;
20231 83792dd9bcc1 Use in-place operators in m-files where possible. Rik <rik@octave.org> parents: 20160 diff changeset	155 wts /= 2;
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	156
11471 994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files. Rik <octave@nomad.inbox5.com> parents: 11465 diff changeset	157 n = length (A);
11458 93a039fe681e logm: style fixes John W. Eaton <jwe@octave.org> parents: 11033 diff changeset	158 s = zeros (n);
93a039fe681e logm: style fixes John W. Eaton <jwe@octave.org> parents: 11033 diff changeset	159 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	160 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	161 endfor
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 endfunction
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	164
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	165 ######################################################################
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	166 ## 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	167 ## [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	168 ## for N-point Gauss-Legendre quadrature.
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	169
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	170 ## Reference:
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	171 ## 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	172 ## 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	173
11458 93a039fe681e logm: style fixes John W. Eaton <jwe@octave.org> parents: 11033 diff changeset	174 function [x, w] = gauss_legendre (n)
21758 ffad2baa90f7 maint: Use newlines to make code more readable. Rik <rik@octave.org> parents: 20852 diff changeset	175
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	176 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	177 v = i./sqrt ((2*i).^2-1);
11458 93a039fe681e logm: style fixes John W. Eaton <jwe@octave.org> parents: 11033 diff changeset	178 [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	179 x = diag (D);
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	180 w = 2*(V(1,:)'.^2);
21758 ffad2baa90f7 maint: Use newlines to make code more readable. Rik <rik@octave.org> parents: 20852 diff changeset	181
10825 cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	182 endfunction
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	183
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy) Jaroslav Hajek <highegg@gmail.com> parents: 7017 diff changeset	184
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14327 diff changeset	185 %!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	186 %!test
e5facc6eec13 Silence warning messages in %!test code. Rik <rik@octave.org> parents: 17744 diff changeset	187 %! warning ("off", "Octave:logm:non-principal", "local");
e5facc6eec13 Silence warning messages in %!test code. Rik <rik@octave.org> parents: 17744 diff changeset	188 %! 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	189 %!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	190 %!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	191 %!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	192 %!assert (full (logm (10eye (3))), logm (full (10eye (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	193 %!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	194 %!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	195 %! 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	196 %! 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	197 %! -1.9769, -1.0922, -0.5831];
29805 20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	198 %! if (ismac ())
20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	199 %! ## The math libraries on macOS seem to require larger tolerances
20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	200 %! tol = 60*eps;
20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	201 %! else
20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	202 %! tol = 40*eps;
20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	203 %! endif
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	204 %! warning ("off", "Octave:logm:non-principal", "local");
29805 20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	205 %! assert (expm (logm (A)), A, tol);
31551 fd29c7a50a78 maint: use commas, semicolons consistently with Octave conventions. Rik <rik@octave.org> parents: 31253 diff changeset	206 %!assert (expm (logm (eye (3))), eye (3))
fd29c7a50a78 maint: use commas, semicolons consistently with Octave conventions. Rik <rik@octave.org> parents: 31253 diff changeset	207 %!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	208
19833 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19697 diff changeset	209 ## Test input validation
28896 90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin. Rik <rik@octave.org> parents: 28789 diff changeset	210 %!error <Invalid call> logm ()
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14327 diff changeset	211 %!error <logm: A must be a square matrix> logm ([1 0;0 1; 2 2])

Mercurial > octave

annotate scripts/linear-algebra/logm.m @ 31551:fd29c7a50a78 stable