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