Mercurial > matrix-functions

comparison mftoolbox/rootpm_schur_newton.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function [X,Y] = rootpm_schur_newton(A,p)
2 %ROOTPM_SCHUR_NEWTON Matrix pth root by Schur-Newton method.
3 % [X,Y] = ROOTPM_SCHUR_NEWTON(A,p) computes the principal pth root
4 % X of the real matrix A using the Schur-Newton algorithm.
5 % It also returns Y = INV(X).
6
7 if norm(imag(A),1), error('A must be real.'), end
8 A = real(A); % Discard any zero imaginary part.
9 [Q,R] = schur(A,'real'); % Quasitriangular R.
10
11 e = eig(R);
12 if any (e(find(e == real(e))) < 0 )
13 error('A has a negative real eigenvalue: principal pth root not defined')
14 end
15
16 f = factor(p);
17 k0 = length(find(f == 2)); % Number of factors 2.
18 q = p/2^(k0);
19 k1 = k0;
20 if q > 1
21
22 emax = max(abs(e)); emin = min(abs(e));
23 if emax > emin % Avoid log(0).
24 k1 = max(k1, ceil( log2( log2(emax/emin) ) ));
25 end
26
27 max_arg = norm(angle(e),inf);
28 if max_arg > pi/8
29 k3 = 1;
30 if max_arg > pi/2
31 k3 = 3;
32 elseif max_arg > pi/4
33 k3 = 2;
34 end
35 k1 = max(k1,k3);
36 end
37
38 end
39
40 for i = 1:k1, R = sqrtm_real(R); end
41
42
43 if q ~= 1
44 pw = 2^(-k1); emax = emax^pw; emin = emin^pw;
45 if ~any(imag(e))
46 % Real eigenvalues.
47 if emax > emin
48 alpha = emax/emin;
49 c = ( (alpha^(1/q)*emax-emin)/( (alpha^(1/q)-1)*(q+1) ) )^(1/q);
50 else
51 c = emin^(1/q);
52 end
53 else
54 % Complex eigenvalues.
55 c = (( emax+emin)/2 )^(1/q);
56 end
57
58 X = rootpm_newton(R,q,c);
59 for i = 1:k1-k0
60 X = X*X;
61 end
62 else
63 X = R;
64 end
65 Y = Q*(X\Q'); % Return inverse pth root, too.
66 X = Q*X*Q';