function [X, iter] = logm_isst(T, prnt)
%LOGM_ISST   Log of triangular matrix by Schur-Pade method with scaling.
%         X = LOGM_ISST(A) computes the logarithm of an upper triangular
%         matrix A, for a matrix with no nonpositive real eigenvalues,
%         using the inverse scaling and squaring method with Pade
%         approximation.  TOL is an error tolerance.
%         [X, ITER] = LOGM_ISST(A, PRNT) returns the number ITER of square
%         roots computed and prints this information if PRNT is nonzero.

% References:
% S. H. Cheng, N. J. Higham, C. S. Kenney, and A. J. Laub, Approximating the
%    logarithm of a matrix to specified accuracy, SIAM J. Matrix Anal. Appl.,
%    22(4):1112-1125, 2001.
% N. J. Higham, Evaluating Pade approximants of the matrix logarithm,
%    SIAM J. Matrix Anal. Appl., 22(4):1126-1135, 2001.

if nargin < 2, prnt = 0; end
n = length(T);

if any( imag(diag(T)) == 0 & real(diag(T)) <= 0 )
   error('A must not have nonpositive real eigenvalues!')
end

if n == 1, X = log(T); iter = 0; return, end

R = T;
maxlogiter = 50;

for iter = 0:maxlogiter

    phi  = norm(T-eye(n),'fro');

    if phi <= 0.25
       if prnt, fprintf('LOGM_ISST computed %g square roots. \n', iter), end
       break
    end
    if iter == maxlogiter, error('Too many square roots in LOGM_ISST.\n'), end

    % Compute upper triangular square root R of T, a column at a time.
    for j=1:n
        R(j,j) = sqrt(T(j,j));
        for i=j-1:-1:1
            R(i,j) = (T(i,j) - R(i,i+1:j-1)*R(i+1:j-1,j))/(R(i,i) + R(j,j));
        end
    end
    T = R;
end

X = 2^(iter)*logm_pf(T-eye(n),8);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S = logm_pf(A,m)
%LOGM_PF   Pade approximation to matrix log by partial fraction expansion.
%          Y = LOGM_PF(A,m) approximates LOG(I+A).

[nodes,wts] = gauss_legendre(m);
% Convert from [-1,1] to [0,1].
nodes = (nodes + 1)/2;
wts = wts/2;

n = length(A);
S = zeros(n);

for j=1:m
    S = S + wts(j)*(A/(eye(n) + nodes(j)*A));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,w] = gauss_legendre(n)
%GAUSS_LEGENDRE  Nodes and weights for Gauss-Legendre quadrature.

% Reference:
% G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature
% rules, Math. Comp., 23(106):221-230, 1969.

i = 1:n-1;
v = i./sqrt((2*i).^2-1);
[V,D] = eig( diag(v,-1)+diag(v,1) );
x = diag(D);
w = 2*(V(1,:)'.^2);