annotate scripts/linear-algebra/logm.m @ 14363:f3d52523cde1

Use Octave coding conventions in all m-file %!test blocks * wavread.m, acosd.m, acot.m, acotd.m, acoth.m, acsc.m, acscd.m, acsch.m, asec.m, asecd.m, asech.m, asind.m, atand.m, cosd.m, cot.m, cotd.m, coth.m, csc.m, cscd.m, csch.m, sec.m, secd.m, sech.m, sind.m, tand.m, accumarray.m, accumdim.m, bitcmp.m, bitget.m, bitset.m, blkdiag.m, cart2pol.m, cart2sph.m, celldisp.m, chop.m, circshift.m, colon.m, common_size.m, cplxpair.m, cumtrapz.m, curl.m, dblquad.m, deal.m, divergence.m, flipdim.m, fliplr.m, flipud.m, genvarname.m, gradient.m, idivide.m, int2str.m, interp1.m, interp1q.m, interp2.m, interp3.m, interpft.m, interpn.m, isa.m, isdir.m, isequal.m, isequalwithequalnans.m, issquare.m, logspace.m, nargchk.m, narginchk.m, nargoutchk.m, nextpow2.m, nthargout.m, num2str.m, pol2cart.m, polyarea.m, postpad.m, prepad.m, profile.m, profshow.m, quadgk.m, quadv.m, randi.m, rat.m, repmat.m, rot90.m, rotdim.m, shift.m, shiftdim.m, sph2cart.m, structfun.m, trapz.m, triplequad.m, convhull.m, dsearch.m, dsearchn.m, griddata3.m, griddatan.m, rectint.m, tsearchn.m, __makeinfo__.m, doc.m, get_first_help_sentence.m, help.m, type.m, unimplemented.m, which.m, imread.m, imwrite.m, dlmwrite.m, fileread.m, is_valid_file_id.m, strread.m, textread.m, textscan.m, commutation_matrix.m, cond.m, condest.m, cross.m, duplication_matrix.m, expm.m, housh.m, isdefinite.m, ishermitian.m, issymmetric.m, logm.m, normest.m, null.m, onenormest.m, orth.m, planerot.m, qzhess.m, rank.m, rref.m, trace.m, vech.m, ans.m, bincoeff.m, bug_report.m, bzip2.m, comma.m, compare_versions.m, computer.m, edit.m, fileparts.m, fullfile.m, getfield.m, gzip.m, info.m, inputname.m, isappdata.m, isdeployed.m, ismac.m, ispc.m, isunix.m, list_primes.m, ls.m, mexext.m, namelengthmax.m, news.m, orderfields.m, paren.m, recycle.m, rmappdata.m, semicolon.m, setappdata.m, setfield.m, substruct.m, symvar.m, ver.m, version.m, warning_ids.m, xor.m, fminbnd.m, fsolve.m, fzero.m, lsqnonneg.m, optimset.m, pqpnonneg.m, sqp.m, matlabroot.m, __gnuplot_drawnow__.m, __plt_get_axis_arg__.m, ancestor.m, cla.m, clf.m, close.m, colorbar.m, colstyle.m, comet3.m, contourc.m, figure.m, gca.m, gcbf.m, gcbo.m, gcf.m, ginput.m, graphics_toolkit.m, gtext.m, hggroup.m, hist.m, hold.m, isfigure.m, ishghandle.m, ishold.m, isocolors.m, isonormals.m, isosurface.m, isprop.m, legend.m, line.m, loglog.m, loglogerr.m, meshgrid.m, ndgrid.m, newplot.m, orient.m, patch.m, plot3.m, plotyy.m, __print_parse_opts__.m, quiver3.m, refreshdata.m, ribbon.m, semilogx.m, semilogxerr.m, semilogy.m, stem.m, stem3.m, subplot.m, title.m, uigetfile.m, view.m, whitebg.m, compan.m, conv.m, deconv.m, mkpp.m, mpoles.m, pchip.m, poly.m, polyaffine.m, polyder.m, polyfit.m, polygcd.m, polyint.m, polyout.m, polyval.m, polyvalm.m, ppder.m, ppint.m, ppjumps.m, ppval.m, residue.m, roots.m, spline.m, intersect.m, ismember.m, powerset.m, setdiff.m, setxor.m, union.m, unique.m, autoreg_matrix.m, bartlett.m, blackman.m, detrend.m, fftconv.m, fftfilt.m, fftshift.m, freqz.m, hamming.m, hanning.m, ifftshift.m, sinc.m, sinetone.m, sinewave.m, unwrap.m, bicg.m, bicgstab.m, gmres.m, gplot.m, nonzeros.m, pcg.m, pcr.m, spaugment.m, spconvert.m, spdiags.m, speye.m, spfun.m, spones.m, sprand.m, sprandsym.m, spstats.m, spy.m, svds.m, treelayout.m, bessel.m, beta.m, betaln.m, factor.m, factorial.m, isprime.m, lcm.m, legendre.m, nchoosek.m, nthroot.m, perms.m, pow2.m, primes.m, reallog.m, realpow.m, realsqrt.m, hadamard.m, hankel.m, hilb.m, invhilb.m, magic.m, rosser.m, vander.m, __finish__.m, center.m, cloglog.m, corr.m, cov.m, gls.m, histc.m, iqr.m, kendall.m, kurtosis.m, logit.m, mahalanobis.m, mean.m, meansq.m, median.m, mode.m, moment.m, ols.m, ppplot.m, prctile.m, probit.m, quantile.m, range.m, ranks.m, run_count.m, runlength.m, skewness.m, spearman.m, statistics.m, std.m, table.m, var.m, zscore.m, betacdf.m, betainv.m, betapdf.m, betarnd.m, binocdf.m, binoinv.m, binopdf.m, binornd.m, cauchy_cdf.m, cauchy_inv.m, cauchy_pdf.m, cauchy_rnd.m, chi2cdf.m, chi2inv.m, chi2pdf.m, chi2rnd.m, discrete_cdf.m, discrete_inv.m, discrete_pdf.m, discrete_rnd.m, empirical_cdf.m, empirical_inv.m, empirical_pdf.m, empirical_rnd.m, expcdf.m, expinv.m, exppdf.m, exprnd.m, fcdf.m, finv.m, fpdf.m, frnd.m, gamcdf.m, gaminv.m, gampdf.m, gamrnd.m, geocdf.m, geoinv.m, geopdf.m, geornd.m, hygecdf.m, hygeinv.m, hygepdf.m, hygernd.m, kolmogorov_smirnov_cdf.m, laplace_cdf.m, laplace_inv.m, laplace_pdf.m, laplace_rnd.m, logistic_cdf.m, logistic_inv.m, logistic_pdf.m, logistic_rnd.m, logncdf.m, logninv.m, lognpdf.m, lognrnd.m, nbincdf.m, nbininv.m, nbinpdf.m, nbinrnd.m, normcdf.m, norminv.m, normpdf.m, normrnd.m, poisscdf.m, poissinv.m, poisspdf.m, poissrnd.m, stdnormal_cdf.m, stdnormal_inv.m, stdnormal_pdf.m, stdnormal_rnd.m, tcdf.m, tinv.m, tpdf.m, trnd.m, unidcdf.m, unidinv.m, unidpdf.m, unidrnd.m, unifcdf.m, unifinv.m, unifpdf.m, unifrnd.m, wblcdf.m, wblinv.m, wblpdf.m, wblrnd.m, kolmogorov_smirnov_test.m, kruskal_wallis_test.m, base2dec.m, bin2dec.m, blanks.m, cstrcat.m, deblank.m, dec2base.m, dec2bin.m, dec2hex.m, findstr.m, hex2dec.m, index.m, isletter.m, mat2str.m, rindex.m, str2num.m, strcat.m, strjust.m, strmatch.m, strsplit.m, strtok.m, strtrim.m, strtrunc.m, substr.m, validatestring.m, demo.m, example.m, fail.m, speed.m, addtodate.m, asctime.m, clock.m, ctime.m, date.m, datenum.m, datetick.m, datevec.m, eomday.m, etime.m, is_leap_year.m, now.m: Use Octave coding conventions in all m-file %!test blocks
author Rik <octave@nomad.inbox5.com>
date Mon, 13 Feb 2012 07:29:44 -0800
parents 4d917a6a858b
children 6c3441f3146b
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
14138
72c96de7a403 maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents: 13097
diff changeset
1 ## Copyright (C) 2008-2012 N.J. Higham
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
2 ## Copyright (C) 2010 Richard T. Guy <guyrt7@wfu.edu>
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
3 ## Copyright (C) 2010 Marco Caliari <marco.caliari@univr.it>
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
4 ##
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
5 ## This file is part of Octave.
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
6 ##
10827
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
7 ## Octave is free software; you can redistribute it and/or modify it
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
8 ## under the terms of the GNU General Public License as published by
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
9 ## the Free Software Foundation; either version 3 of the License, or (at
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
10 ## your option) any later version.
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
11 ##
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
12 ## 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
13 ## WITHOUT ANY WARRANTY; without even the implied warranty of
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
14 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
15 ## General Public License for more details.
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
16 ##
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
17 ## 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
18 ## along with Octave; see the file COPYING. If not, see
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
19 ## <http://www.gnu.org/licenses/>.
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
20
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
21 ## -*- texinfo -*-
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
22 ## @deftypefn {Function File} {@var{s} =} logm (@var{A})
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
23 ## @deftypefnx {Function File} {@var{s} =} logm (@var{A}, @var{opt_iters})
10827
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
24 ## @deftypefnx {Function File} {[@var{s}, @var{iters}] =} logm (@dots{})
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
25 ## Compute the matrix logarithm of the square matrix @var{A}. The
10827
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
26 ## 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
27 ##
10827
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
28 ## @example
14327
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
29 ## 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
30 ## @end example
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
31 ##
10827
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
32 ## The optional argument @var{opt_iters} is the maximum number of square roots
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
33 ## to compute and defaults to 100. The optional output @var{iters} is the
228cd18455a6 logm.m: Improve documentation string. Add GPL header. Add additional test block.
Rik <octave@nomad.inbox5.com>
parents: 10825
diff changeset
34 ## number of square roots actually computed.
12584
7ef7e20057fa Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
35 ## @seealso{expm, sqrtm}
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
36 ## @end deftypefn
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
37
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
38 ## 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
39 ## (SIAM, 2008.)
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
40 ##
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
41
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
42 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
43
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
44 if (nargin == 0 || nargin > 2)
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5307
diff changeset
45 print_usage ();
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
46 endif
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
47
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
48 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
49 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
50 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
51
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
52 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
53 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
54 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
55 elseif (strfind (typeinfo (A), "diagonal matrix"))
ff7e0776ba0f logm.m: handle scalar and diagonal matrix arguments specially; call log_pade_pf only if m > 1
John W. Eaton <jwe@octave.org>
parents: 11472
diff changeset
56 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
57 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
58 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
59
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
60 [u, s] = schur (A);
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
61
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
62 if (isreal (A))
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
63 [u, s] = rsf2csf (u, s);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
64 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
65
13097
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
66 eigv = diag (s);
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
67 if (any (eigv < 0))
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
68 warning ("Octave:logm:non-principal",
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
69 "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
70 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
71
13097
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
72 real_eig = all (eigv >= 0);
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
73
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
74 k = 0;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
75 ## 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
76 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
77 p = 0;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
78 m = 7;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
79 while (k < opt_iters)
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
80 tau = norm (s - eye (size (s)),1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
81 if (tau <= theta (7))
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
82 p = p + 1;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
83 j(1) = find (tau <= theta, 1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
84 j(2) = find (tau / 2 <= theta, 1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
85 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
86 m = j(1);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
87 break
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
88 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
89 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
90 k = k + 1;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
91 s = sqrtm (s);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
92 endwhile
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
93
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
94 if (k >= opt_iters)
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
95 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
96 endif
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
97
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
98 s = s - eye (size (s));
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
99
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
100 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
101 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
102 endif
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
103
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
104 s = 2^k * u * s * u';
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
105
13097
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
106 ## Remove small complex values (O(eps)) which may have entered calculation
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
107 if (real_eig)
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
108 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
109 endif
52e4aa30d5b2 logm.m: Return real matrix when all eigenvalues are real (Bug #32121).
Rik <octave@nomad.inbox5.com>
parents: 12584
diff changeset
110
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
111 if (nargout == 2)
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
112 iters = k;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
113 endif
4334
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
114
764229f9a5c8 [project @ 2003-02-19 06:24:02 by jwe]
jwe
parents:
diff changeset
115 endfunction
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 ################## ANCILLARY FUNCTIONS ################################
11033
d9c8916bb9dd Untabify a few remaining .m scripts.
Rik <octave@nomad.inbox5.com>
parents: 10827
diff changeset
118 ###### 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
119 ###### Reference:
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
120 ###### 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
121 ###### (SIAM, 2008.).
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
122 #######################################################################
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
123
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
124 ##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
125 ## 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
126 ## 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
127
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
128 function s = logm_pade_pf (A, m)
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
129 [nodes, wts] = gauss_legendre (m);
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
130 ## Convert from [-1,1] to [0,1].
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
131 nodes = (nodes+1)/2;
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
132 wts = wts/2;
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
133
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 11465
diff changeset
134 n = length (A);
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
135 s = zeros (n);
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
136 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
137 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
138 endfor
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
139 endfunction
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
140
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
141 ######################################################################
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
142 ## 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
143 ## [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
144 ## for N-point Gauss-Legendre quadrature.
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
145
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
146 ## Reference:
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
147 ## 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
148 ## 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
149
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
150 function [x, w] = gauss_legendre (n)
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
151 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
152 v = i./sqrt ((2*i).^2-1);
11458
93a039fe681e logm: style fixes
John W. Eaton <jwe@octave.org>
parents: 11033
diff changeset
153 [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
154 x = diag (D);
10825
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
155 w = 2*(V(1,:)'.^2);
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
156 endfunction
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
157
cace99cb01ab rewrite logm (M. Caliari, R.T. Guy)
Jaroslav Hajek <highegg@gmail.com>
parents: 7017
diff changeset
158
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
159 %!assert (norm (logm ([1 -1;0 1]) - [0 -1; 0 0]) < 1e-5)
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
160 %!assert (norm (expm (logm ([-1 2 ; 4 -1])) - [-1 2 ; 4 -1]) < 1e-5)
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
161 %!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
162 %!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
163 %!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
164 %!assert (full (logm (10*eye (3))), logm (full (10*eye (3))), 8*eps)
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
165
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
166 %% Test input validation
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
167 %!error logm ()
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
168 %!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
169 %!error <logm: A must be a square matrix> logm ([1 0;0 1; 2 2])
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
170