Mercurial > matrix-functions
diff funm_files/myfunm.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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/funm_files/myfunm.m Wed May 06 14:56:53 2015 +0200 @@ -0,0 +1,96 @@ +function [F,n_swaps,n_calls,terms,ind,T] = myfunm(A,fun,delta,tol,prnt,m) +%MYFUNM Evaluate general matrix function. +% F = FUNM(A,FUN), for a square matrix argument A, evaluates the +% function FUN at the square matrix A. +% FUN(X,K) must return the K'th derivative of +% the function represented by FUN evaluated at the vector X. +% The MATLAB functions COS, SIN, EXP, LOG can be passed as FUN, +% i.e. FUNM(A,@COS), FUNM(A,@SIN), FUNM(@EXP), FUNM(A,@LOG). +% For matrix square roots use SQRTM(A) instead. +% For matrix exponentials, either of EXPM(A) and FUNM(A,@EXPM) +% may be the faster or the more accurate, depending on A. +% +% F = FUNM(A,FUN,DELTA,TOL,PRNT,M) specifies a tolerance +% DELTA used in determining the blocking (default 0.1), +% and a tolerance TOL used in a convergence test for evaluating the +% Taylor series (default EPS). +% If PRNT is nonzero then information describing the +% behaviour of the algorithm is printed. +% M, if supplied, defines a blocking. +% +% [F,N_SWAPS,N_CALLS,TERMS,IND,T] = FUNM(A,FUN,...) also returns +% N_SWAPS: the total number of swaps in the Schur re-ordering. +% N_CALLS: the total number of calls to ZTREXC for the re-ordering. +% TERMS(I): the number of Taylor series terms used when evaluating +% the I'th atomic triangular block. +% IND: a cell array specifying the blocking: the (I,J) block of +% the re-ordered Schur factor T is T(IND{I},IND{j}), +% T: the re-ordered Schur form. + +if isequal(fun,@cos) || isequal(fun,'cos'), fun = @fun_cos; end +if isequal(fun,@sin) || isequal(fun,'sin'), fun = @fun_sin; end +if isequal(fun,@exp) || isequal(fun,'exp'), fun = @fun_exp; end + +if nargin < 3 || isempty(delta), delta = 0.1; end +if nargin < 4 || isempty(tol), tol = eps; end +if nargin < 5 || isempty(prnt), prnt = 0; end +if nargin < 6, m = []; end + +n = length(A); + +% First form complex Schur form (if A not already upper triangular). +if isequal(A,triu(A)) + T = A; U = eye(n); +else + [U,T] = schur(A,'complex'); +end + +if isequal(T,tril(T)) % Handle special case of diagonal T. + F = U*diag(feval(fun,diag(T)))*U'; + n_swaps = 0; n_calls = 0; terms = 0; ind = {1:n}; + return +end + +% Determine reordering of Schur form into block form. +if isempty(m), m = blocking(T,delta,abs(prnt)>=3); end + +if prnt, fprintf('delta (blocking) = %9.2e, tol (TS) = %9.2e\n', delta, tol), +end + +[M,ind,n_swaps] = swapping(m); +n_calls = size(M,1); +if n_calls > 0 % If there are swaps to do... + [U,T] = swap(U,T,M); % MEX file +end + +m = length(ind); + +% Calculate F(T) +F = zeros(n); + +for col=1:m + j = ind{col}; + [F(j,j), n_terms] = funm_atom(T(j,j),fun,tol,abs(prnt)*(prnt ~= 1)); + terms(col) = n_terms; + + for row=col-1:-1:1 + i = ind{row}; + if length(i) == 1 && length(j) == 1 + % Scalar case. + k = i+1:j-1; + temp = T(i,j)*(F(i,i) - F(j,j)) + F(i,k)*T(k,j) - T(i,k)*F(k,j); + F(i,j) = temp/(T(i,i)-T(j,j)); + else + k = cat(2,ind{row+1:col-1}); + 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); + F(i,j) = sylv_tri(T(i,i),-T(j,j),rhs); + end + end +end + +F = U*F*U'; + +% As in FUNM: +if isreal(A) && norm(imag(F),1) <= 10*n*eps*norm(F,1) + F = real(F); +end