Mercurial > octave-nkf

comparison scripts/linear-algebra/logm.m @ 20266:83792dd9bcc1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use in-place operators in m-files where possible. * scripts/audio/@audioplayer/set.m, scripts/audio/@audiorecorder/set.m, scripts/audio/mu2lin.m, scripts/elfun/cosd.m, scripts/general/del2.m, scripts/general/profexplore.m, scripts/general/quadl.m, scripts/general/rat.m, scripts/general/rotdim.m, scripts/help/get_first_help_sentence.m, scripts/help/private/__strip_html_tags__.m, scripts/image/cubehelix.m, scripts/io/textread.m, scripts/linear-algebra/duplication_matrix.m, scripts/linear-algebra/housh.m, scripts/linear-algebra/krylov.m, scripts/linear-algebra/logm.m, scripts/linear-algebra/normest.m, scripts/linear-algebra/onenormest.m, scripts/optimization/fminsearch.m, scripts/optimization/lsqnonneg.m, scripts/optimization/qp.m, scripts/plot/appearance/annotation.m, scripts/plot/appearance/axis.m, scripts/plot/appearance/legend.m, scripts/plot/appearance/specular.m, scripts/plot/draw/colorbar.m, scripts/plot/draw/hist.m, scripts/plot/draw/plotmatrix.m, scripts/plot/draw/private/__stem__.m, scripts/plot/util/__actual_axis_position__.m, scripts/plot/util/__gnuplot_drawnow__.m, scripts/plot/util/findobj.m, scripts/plot/util/print.m, scripts/plot/util/private/__go_draw_axes__.m, scripts/plot/util/private/__print_parse_opts__.m, scripts/plot/util/rotate.m, scripts/polynomial/pchip.m, scripts/polynomial/polyaffine.m, scripts/polynomial/polyder.m, scripts/polynomial/private/__splinefit__.m, scripts/polynomial/residue.m, scripts/signal/arch_fit.m, scripts/signal/arch_rnd.m, scripts/signal/bartlett.m, scripts/signal/blackman.m, scripts/signal/freqz.m, scripts/signal/hamming.m, scripts/signal/hanning.m, scripts/signal/spectral_adf.m, scripts/signal/spectral_xdf.m, scripts/signal/stft.m, scripts/sparse/bicgstab.m, scripts/sparse/cgs.m, scripts/sparse/private/__sprand_impl__.m, scripts/sparse/qmr.m, scripts/sparse/sprandsym.m, scripts/sparse/svds.m, scripts/specfun/legendre.m, scripts/special-matrix/gallery.m, scripts/statistics/base/gls.m, scripts/statistics/models/logistic_regression.m, scripts/statistics/tests/kruskal_wallis_test.m, scripts/statistics/tests/manova.m, scripts/statistics/tests/wilcoxon_test.m, scripts/time/datevec.m: Use in-place operators in m-files where possible.
author Rik <rik@octave.org>
date Tue, 26 May 2015 21:07:42 -0700
parents 03b9d17a2d95
children
comparison
equal deleted inserted replaced
20265:e914b5399c67 20266:83792dd9bcc1
84 p = 0; 84 p = 0;
85 m = 7; 85 m = 7;
86 while (k < opt_iters) 86 while (k < opt_iters)
87 tau = norm (s - eye (size (s)),1); 87 tau = norm (s - eye (size (s)),1);
88 if (tau <= theta (7)) 88 if (tau <= theta (7))
89 p = p + 1; 89 p += 1;
90 j(1) = find (tau <= theta, 1); 90 j(1) = find (tau <= theta, 1);
91 j(2) = find (tau / 2 <= theta, 1); 91 j(2) = find (tau / 2 <= theta, 1);
92 if (j(1) - j(2) <= 1 || p == 2) 92 if (j(1) - j(2) <= 1 || p == 2)
93 m = j(1); 93 m = j(1);
94 break 94 break
95 endif 95 endif
96 endif 96 endif
97 k = k + 1; 97 k += 1;
98 s = sqrtm (s); 98 s = sqrtm (s);
99 endwhile 99 endwhile
100 100
101 if (k >= opt_iters) 101 if (k >= opt_iters)
102 warning ("logm: maximum number of square roots exceeded; results may still be accurate"); 102 warning ("logm: maximum number of square roots exceeded; results may still be accurate");
103 endif 103 endif
104 104
105 s = s - eye (size (s)); 105 s -= eye (size (s));
106 106
107 if (m > 1) 107 if (m > 1)
108 s = logm_pade_pf (s, m); 108 s = logm_pade_pf (s, m);
109 endif 109 endif
110 110
134 134
135 function s = logm_pade_pf (A, m) 135 function s = logm_pade_pf (A, m)
136 [nodes, wts] = gauss_legendre (m); 136 [nodes, wts] = gauss_legendre (m);
137 ## Convert from [-1,1] to [0,1]. 137 ## Convert from [-1,1] to [0,1].
138 nodes = (nodes+1)/2; 138 nodes = (nodes+1)/2;
139 wts = wts/2; 139 wts /= 2;
140 140
141 n = length (A); 141 n = length (A);
142 s = zeros (n); 142 s = zeros (n);
143 for j = 1:m 143 for j = 1:m
144 s += wts(j)*(A/(eye (n) + nodes(j)*A)); 144 s += wts(j)*(A/(eye (n) + nodes(j)*A));