Mercurial > matrix-functions
comparison toolbox/matsignt.m @ 2:c124219d7bfa draft
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
| author | Antonio Pino Robles <data.script93@gmail.com> |
|---|---|
| date | Thu, 07 May 2015 18:36:24 +0200 |
| parents | 8f23314345f4 |
| children |
comparison
equal
deleted
inserted
replaced
| 1:e471a92d17be | 2:c124219d7bfa |
|---|---|
| 1 function S = matsignt(T) | |
| 2 %MATSIGNT Matrix sign function of a triangular matrix. | |
| 3 % S = MATSIGN(T) computes the matrix sign function S of the | |
| 4 % upper triangular matrix T using a recurrence. | |
| 5 | |
| 6 % Adapted from FUNM. Called by SIGNM. | |
| 7 | |
| 8 if norm(tril(T,-1),1), error('Matrix must be upper triangular.'), end | |
| 9 | |
| 10 n = max(size(T)); | |
| 11 | |
| 12 S = diag( sign( diag(real(T)) ) ); | |
| 13 tol = 0; | |
| 14 for p = 1:n-1 | |
| 15 for i = 1:n-p | |
| 16 | |
| 17 j = i+p; | |
| 18 d = T(j,j) - T(i,i); | |
| 19 | |
| 20 if S(i,i) ~= -S(j,j) % Solve via S^2 = I if we can. | |
| 21 | |
| 22 % Get S(i,j) from S^2 = I. | |
| 23 k = i+1:j-1; | |
| 24 RHS = 0; | |
| 25 if k, RHS = RHS - S(i,k)*S(k,j); end | |
| 26 S(i,j) = RHS / (S(i,i)+S(j,j)); | |
| 27 | |
| 28 else | |
| 29 | |
| 30 % Get S(i,j) from S*T = T*S. | |
| 31 s = T(i,j)*(S(j,j)-S(i,i)); | |
| 32 if p > 1 | |
| 33 k = i+1:j-1; | |
| 34 s = s + T(i,k)*S(k,j) - S(i,k)*T(k,j); | |
| 35 end | |
| 36 S(i,j) = s/d; | |
| 37 | |
| 38 end | |
| 39 | |
| 40 end | |
| 41 end |