Mercurial > matrix-functions

comparison funm_files/myfunm.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 [F,n_swaps,n_calls,terms,ind,T] = myfunm(A,fun,delta,tol,prnt,m)
2 %MYFUNM Evaluate general matrix function.
3 % F = FUNM(A,FUN), for a square matrix argument A, evaluates the
4 % function FUN at the square matrix A.
5 % FUN(X,K) must return the K'th derivative of
6 % the function represented by FUN evaluated at the vector X.
7 % The MATLAB functions COS, SIN, EXP, LOG can be passed as FUN,
8 % i.e. FUNM(A,@COS), FUNM(A,@SIN), FUNM(@EXP), FUNM(A,@LOG).
9 % For matrix square roots use SQRTM(A) instead.
10 % For matrix exponentials, either of EXPM(A) and FUNM(A,@EXPM)
11 % may be the faster or the more accurate, depending on A.
12 %
13 % F = FUNM(A,FUN,DELTA,TOL,PRNT,M) specifies a tolerance
14 % DELTA used in determining the blocking (default 0.1),
15 % and a tolerance TOL used in a convergence test for evaluating the
16 % Taylor series (default EPS).
17 % If PRNT is nonzero then information describing the
18 % behaviour of the algorithm is printed.
19 % M, if supplied, defines a blocking.
20 %
21 % [F,N_SWAPS,N_CALLS,TERMS,IND,T] = FUNM(A,FUN,...) also returns
22 % N_SWAPS: the total number of swaps in the Schur re-ordering.
23 % N_CALLS: the total number of calls to ZTREXC for the re-ordering.
24 % TERMS(I): the number of Taylor series terms used when evaluating
25 % the I'th atomic triangular block.
26 % IND: a cell array specifying the blocking: the (I,J) block of
27 % the re-ordered Schur factor T is T(IND{I},IND{j}),
28 % T: the re-ordered Schur form.
29
30 if isequal(fun,@cos) || isequal(fun,'cos'), fun = @fun_cos; end
31 if isequal(fun,@sin) || isequal(fun,'sin'), fun = @fun_sin; end
32 if isequal(fun,@exp) || isequal(fun,'exp'), fun = @fun_exp; end
33
34 if nargin < 3 || isempty(delta), delta = 0.1; end
35 if nargin < 4 || isempty(tol), tol = eps; end
36 if nargin < 5 || isempty(prnt), prnt = 0; end
37 if nargin < 6, m = []; end
38
39 n = length(A);
40
41 % First form complex Schur form (if A not already upper triangular).
42 if isequal(A,triu(A))
43 T = A; U = eye(n);
44 else
45 [U,T] = schur(A,'complex');
46 end
47
48 if isequal(T,tril(T)) % Handle special case of diagonal T.
49 F = U*diag(feval(fun,diag(T)))*U';
50 n_swaps = 0; n_calls = 0; terms = 0; ind = {1:n};
51 return
52 end
53
54 % Determine reordering of Schur form into block form.
55 if isempty(m), m = blocking(T,delta,abs(prnt)>=3); end
56
57 if prnt, fprintf('delta (blocking) = %9.2e, tol (TS) = %9.2e\n', delta, tol),
58 end
59
60 [M,ind,n_swaps] = swapping(m);
61 n_calls = size(M,1);
62 if n_calls > 0 % If there are swaps to do...
63 [U,T] = swap(U,T,M); % MEX file
64 end
65
66 m = length(ind);
67
68 % Calculate F(T)
69 F = zeros(n);
70
71 for col=1:m
72 j = ind{col};
73 [F(j,j), n_terms] = funm_atom(T(j,j),fun,tol,abs(prnt)*(prnt ~= 1));
74 terms(col) = n_terms;
75
76 for row=col-1:-1:1
77 i = ind{row};
78 if length(i) == 1 && length(j) == 1
79 % Scalar case.
80 k = i+1:j-1;
81 temp = T(i,j)*(F(i,i) - F(j,j)) + F(i,k)*T(k,j) - T(i,k)*F(k,j);
82 F(i,j) = temp/(T(i,i)-T(j,j));
83 else
84 k = cat(2,ind{row+1:col-1});
85 rhs = F(i,i)*T(i,j) - T(i,j)*F(j,j) + F(i,k)*T(k,j) - T(i,k)*F(k,j);
86 F(i,j) = sylv_tri(T(i,i),-T(j,j),rhs);
87 end
88 end
89 end
90
91 F = U*F*U';
92
93 % As in FUNM:
94 if isreal(A) && norm(imag(F),1) <= 10*n*eps*norm(F,1)
95 F = real(F);
96 end