Mercurial > matrix-functions
comparison mftoolbox/logm_iss.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 [X,k,m] = logm_iss(A) | |
| 2 %LOGM_ISS Matrix logarithm by inverse scaling and squaring method. | |
| 3 % X = LOGM_ISS(A) computes the logarithm of A, for a matrix with no | |
| 4 % nonpositive real eigenvalues, using the inverse scaling and squaring | |
| 5 % method with Pade approximation. | |
| 6 % Matrix square roots are computed by the product form of the | |
| 7 % Denman-Beavers teration. | |
| 8 % [X,K,M] = LOGM_ISS(A) returns the number K of square roots | |
| 9 % computed and the degree M of the Pade approximant. | |
| 10 | |
| 11 e = eig(A); | |
| 12 if any( imag(e) == 0 & real(e) <= 0 ) | |
| 13 error('A must not have any nonpositive real eigenvalues!') | |
| 14 end | |
| 15 | |
| 16 n = length(A); | |
| 17 | |
| 18 load log_pade_err_opt % mmax-by-3 matrix DATA. | |
| 19 % mvals = data(:,1); | |
| 20 xvals = data(:,2); | |
| 21 | |
| 22 X = A; | |
| 23 k = 0; p = 0; itk = 5; | |
| 24 | |
| 25 while 1 | |
| 26 | |
| 27 normdiff = norm(X-eye(n),1); | |
| 28 if normdiff <= xvals(16) | |
| 29 | |
| 30 p = p+1; | |
| 31 j1 = find(normdiff <= xvals(3:16)); | |
| 32 j1 = j1(1) + 2; | |
| 33 j2 = find(normdiff/2 <= xvals(3:16)); | |
| 34 j2 = j2(1) + 2; | |
| 35 if 2*(j1-j2)/3 < itk || p == 2, m = j1; break, end | |
| 36 | |
| 37 end | |
| 38 | |
| 39 [X,M,itk] = sqrtm_dbp(X,1); k = k+1; | |
| 40 | |
| 41 end | |
| 42 | |
| 43 X = 2^k*logm_pade_pf(X-eye(n),m); | |
| 44 if isreal(A), X = real(X); end |