function [diagonal,m] = hmat(a, epsilon, maxit)

%           [diagonal,m] = hmat(a,epsilon)
%     Given an nxn matrix A, find, if possible, a diagonal matrix 
%     D such that AD is strictly row diagonally dominant.
%     This is achieved using Algorithm H found in 
%              "AN ITERATIVE CRITERION FOR H-MATRICES"  
%       by B. Li, L. Li, H. Niki,M. Harada, and M. Tsatsomeros. 
%
%   INPUT 
%         a=square matrix
%         epsilon=the parameter of re-distribution (default=0.001)
%         maxit=maximum number of iterations allowed (default=100)
%   OUTPUT 
%         m=number of iterations performed
%         diagonal=the matrix that achieves strict diagonal dominance
%                    (if it exists)
%         (Written by M. Tsatsomeros, Nov. 20, 1996)


n=size(a,1); diagonal=eye(n); m=1; one=ones(1,n); stoppage=0;
if (nargin==1) epsilon=0.001; maxit=100; end
if (nargin==2) maxit=100; end

% check first if diagonal entries are nonzero
if (1-all(diag(a)))
  stoppage=1; diagonal=[]; m=m-1; 'Input is NOT an H-matrix',
end  

while (stoppage==0 & m<maxit+1)
  for i=1:n
        r(i)=sum(abs(a(i,1:n)))-abs(a(i,i));
        if (abs(a(i,i))>r(i))
           d(i)=(r(i)+epsilon)/(abs(a(i,i))+epsilon);
        else
           d(i)=1;
        end
  end   
  if (d==one)
        stoppage=1; diagonal=[]; 'Input is NOT an H-matrix',
  elseif (d<one)
        stoppage=1; 'Input IS an H-matrix', 
  else
        for i=1:n
            diagonal(i,i)=diagonal(i,i)*d(i); 
        end
        a=a*diag(d); m=m+1;              
  end
end
if (m==maxit+1 & stoppage==0) 
   diagonal=[]; m=m-1;
   'Inconclusive: decrease "epsilon" or increase "maxit"', 
end