Mercurial > matrix-functions

annotate toolbox/matsignt.m @ 2:c124219d7bfa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
1 function S = matsignt(T)
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
2 %MATSIGNT Matrix sign function of a triangular matrix.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
3 % S = MATSIGN(T) computes the matrix sign function S of the
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
4 % upper triangular matrix T using a recurrence.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
5
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
6 % Adapted from FUNM. Called by SIGNM.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
7
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
8 if norm(tril(T,-1),1), error('Matrix must be upper triangular.'), end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
9
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
10 n = max(size(T));
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
11
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
12 S = diag( sign( diag(real(T)) ) );
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
13 tol = 0;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
14 for p = 1:n-1
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
15 for i = 1:n-p
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
16
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
17 j = i+p;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
18 d = T(j,j) - T(i,i);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
19
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
20 if S(i,i) ~= -S(j,j) % Solve via S^2 = I if we can.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
21
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
22 % Get S(i,j) from S^2 = I.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
23 k = i+1:j-1;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
24 RHS = 0;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
25 if k, RHS = RHS - S(i,k)*S(k,j); end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
26 S(i,j) = RHS / (S(i,i)+S(j,j));
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
27
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
28 else
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
29
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
30 % Get S(i,j) from S*T = T*S.
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
31 s = T(i,j)*(S(j,j)-S(i,i));
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
32 if p > 1
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
33 k = i+1:j-1;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
34 s = s + T(i,k)*S(k,j) - S(i,k)*T(k,j);
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
35 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
36 S(i,j) = s/d;
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
37
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
38 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
39
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
40 end
8f23314345f4 Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff changeset
41 end