Mercurial > matrix-functions
comparison mftoolbox/logm_frechet_pade.m @ 0:8f23314345f4 draft
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 L = logm_frechet_pade(A,E) | |
| 2 %LOGM_FRECHET_PADE Frechet derivative of matrix logarithm via Pade approx. | |
| 3 % L = LOGM_FRECHET_PADE(A,E) evaluates the Frechet | |
| 4 % derivative of the matrix logarithm at A in the direction E via the | |
| 5 % inverse scaling and squaring and a Pade approximant of the function | |
| 6 % tanh(x)/x. A must have no eigenvalues on the negative real axis. | |
| 7 | |
| 8 real_data = isreal(A) && isreal(E); | |
| 9 % Form complex Schur form if A not already upper triangular. | |
| 10 use_Schur = false; | |
| 11 if ~isequal(A,triu(A)) | |
| 12 use_Schur = true; | |
| 13 [Q,T] = schur(A,'complex'); A = T; E = Q'*E*Q; | |
| 14 end | |
| 15 | |
| 16 I = eye(size(A)); | |
| 17 B = A; | |
| 18 for i = 0:inf | |
| 19 if norm(B-I,1) <= 1-1/exp(1), s = i; break, end | |
| 20 B = sqrtm(B); | |
| 21 Aroot{i+1} = B; | |
| 22 end | |
| 23 | |
| 24 % Positive zeros of p8 and q8 in r8 = p8/q8 Pade approximant. | |
| 25 load tau_r8_zeros | |
| 26 % Zeros come in \pm pairs. | |
| 27 a = complex(0, [pzero; -pzero]); | |
| 28 b = complex(0, [qzero; -qzero]); | |
| 29 | |
| 30 E = 2^s*E; | |
| 31 for i = 1:s | |
| 32 E = sylvsol(Aroot{i},Aroot{i},E); | |
| 33 end | |
| 34 | |
| 35 G = sylvsol(B,B,E); | |
| 36 | |
| 37 X = logm(B); | |
| 38 | |
| 39 for i=8:-1:1 | |
| 40 rhs = (I + X/b(i)) * G + G * (I - X/b(i)); | |
| 41 AA = I + X/a(i); BB = I - X/a(i); | |
| 42 G = sylvsol(AA, BB, rhs); | |
| 43 end | |
| 44 | |
| 45 L = 2*G; | |
| 46 if use_Schur, L = Q*L*Q'; end | |
| 47 if real_data, L = real(L); end |