function x = hanbliro(A,b)
%   Процедура Хансена-Блика-Рона для интервальных линейных систем.
%   HANBLIRO(A,b) выдаёт внешнюю интервальную оценку множества решений 
%   интервальной системы линейных алгебраических уравнений Ax = b.  
% 
%   Входные аргументы 
%       A - интервальная квадратная матрица решаемой системы 
%       b - интервальный вектор правой части решаемой системы
% 
%   Применима для интервальных линейных систем с H-матрицами, иначе
%   результат полагается равным [-inf,+inf] по каждой компоненте.
%   Реализация основана на работе 
%     Neumaier A. A simple derivation of Hansen-Bliek-Rohn-Ning-Kearfott
%     enclosure for linear interval equations // Reliable Computing. -
%     1999. - Vol. 5, No. 2. - P. 131-136.
%
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   обработка входных данных, проверка их корректности
  
[m, n] = size(A);
if ( m~=n )
    error('Матрица системы должна быть квадратной')
end 
[k, l] = size(b);
if ( k~=n || l>1 )
    error('Матрица и правая часть должен иметь одинаковые размеры')
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   проверяем режим округления и устанавливаем его равным 'nearest'
  
e = 1e-30;
if 1+e==1-e 
    rndold = 0;
else
    rndold = getround;
    setround(0);
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   организуем компарант матрицы системы
C = compmat(A);
  
%   выделяем из матрицы диагональ, обращаем компарант и т.д.
dA = diag(A);
E = inv(C);
v = abs(E*ones(n,1));
setround(-1);
u = C*v;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   выполняем собственно процедуру Хансена-Блика-Рона
  
if all( min(u) > 0 ) 
    %   A является H-матрицей
    dC = diag(C); 
    C = C*E - eye(n); 
    setround(1);
    w = zeros(1,n);
    for i = 1:n
        w = max(w,(-C(i,:))/u(i));
    end; 
    E = E + v*w; 
    u = E*mag(b); 
    d = diag(E); 
    alpha = dC + (-1)./d;  
    beta = u./d - mag(b); 
    x = (b + midrad(0,beta))./(dA + midrad(0,alpha));
else 
    %   A не является H-матрицей
    x = infsup(-Inf,Inf)*ones(n,1);
end  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   возвращаем старый режим округления
  
setround(rndold);
  
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  