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

annotate scripts/linear-algebra/logm.m @ 31253:a40c0b7aa376

maint: changes to follow Octave coding conventions. * NEWS.8.md: Wrap lines to 72 chars. * LSODE-opts.in: Use two spaces after sentence ending period. * LSODE.cc: Use minimum of two spaces between code and start of comment. * MemoizedFunction.m: Change copyright date to 2022 since this is the year it was accepted into core. Don't wrap error() lines to 80 chars. Use newlines to improve readability of switch statements. Use minimum of two spaces between code and start of comment. * del2.m, integral.m, interp1.m, interp2.m, griddata.m, inpolygon.m, waitbar.m, cubehelix.m, ind2x.m, importdata.m, textread.m, logm.m, lighting.m, shading.m, xticklabels.m, yticklabels.m, zticklabels.m, colorbar.m, meshc.m, print.m, __gnuplot_draw_axes__.m, struct2hdl.m, ppval.m, ismember.m, iqr.m: Use a space between comment character '#' and start of comment. Use hyphen for adjectives describing dimensions such as "1-D". * vectorize.m, ode23s.m: Use is_function_handle() instead of "isa (x, "function_handle")" for clarity and performance. * clearAllMemoizedCaches.m: Change copyright date to 2022 since this is the year it was accepted into core. Remove input validation which is done by interpreter. Use two newlines between end of code and start of BIST tests. * memoize.m: Change copyright date to 2022 since this is the year it was accepted into core. Re-wrap documentation to 80 chars. Use is_function_handle() instead of "isa (x, "function_handle")" for clarity and performance. Use two newlines between end of code and start of BIST tests. Use semicolon for assert statements within %!test block. Re-write BIST tests for input validation. * __memoize__.m: Change copyright date to 2022 since this is the year it was accepted into core. Use spaces in for statements to improve readability. * unique.m: Add FIXME note to commented BIST test * dec2bin.m: Remove stray newline at end of file. * triplequad.m: Reduce doubly-commented BIST syntax using "#%!#" to "#%!". * delaunayn.m: Use input variable names in error() statements. Use minimum of two spaces between code and start of comment. Use hyphen for describing dimensions. Use two newlines between end of code and start of BIST tests. Update BIST tests to pass.

author	Rik <rik@octave.org>
date	Mon, 03 Oct 2022 18:06:55 -0700
parents	796f54d4ddbf
children	fd29c7a50a78

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);
20380c9bed30 logm.m: Allow larger tolerance for test on macOS. Markus Mützel <markus.muetzel@gmx.de> parents: 29770 diff changeset	206 %!assert (expm (logm (eye (3))), eye (3));
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	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 @ 31253:a40c0b7aa376