function y = hbrgs(A,b)
//  Внешнее интервальное оценивание множества решений для интервальной 
//  системы линейных алгебраических уравнений с H-матрицами с помощью 
//  процедуры Хансена-Блика-Рона и последующее уточнение полученной 
//  оценки с помощью интервального итерационного метода Гаусса-Зейделя. 
//  Если матрица системы не является H-матрицей, результат полагается 
//  равным [-inf,+inf] по каждой компоненте.
//
//  Входные параметры: 
//      A - интервальная матрица решаемой системы уравнений  
//      b - интервальный вектор правой части решаемой системы 
// 
//  Реализация процедуры Хансена-Блика-Рона основана на работе 
//      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.
//
////////////////////////////////////////////////////////////////////////////////
  
Infty = %inf;   //  специальный идентификатор для бесконечности  
IterLim = 16;   //  максимальное количество уточняющих итераций Гаусса-Зейделя 
  
//  ввод данных и проверка их корректности 
  
A = interval(A);
b = interval(b);
  
m = size(A,1);
n = size(A,2);
if ( m~=n ) then
    error('Матрица системы должна быть квадратной!')
end
  
mb = size(b,1);
nb = size(b,2);
if ( mb~=m | nb>1 ) then  
    error('Размеры матрицы и правой части должны соответствовать!')
end 
  
////////////////////////////////////////////////////////////////////////////////
//
//  организуем компарант C для матрицы системы A, 
//  а также некоторые другие вспомогательные объекты 
  
for i = 1:n
    for j = 1:n
        if i==j then 
            C(i,j) = mig(A(i,j))
        else 
            C(i,j) = -mag(A(i,j))
        end
    end 
end
  
dA = diag(A);  
E = interval(inv(C)); 
v = abs(E*ones(n,1));
u = C*v; 
  
//  выполняем процедуру Хансена-Блика-Рона  
  
if ~and( min(inf(u)) > 0 ) then
    //  A не является интервальной H-матрицей
    for i = 1:n
        x(i) = interval(-Infty,Infty)
    end 
else 
    //  A является интервальной H-матрицей 
    dC = diag(C); 
    C = C*E - eye(n,n); 
    w = zeros(1,n); 
    for i = 1:n
        w = max(w,sup(-C(i,:)/u(i)))
    end ; 
    E = E + v*w; 
    u = E*mag(b); 
    d = diag(E); 
    for i = 1:n
        alpha = sup(dC(i) + (-1.)/d(i));  
        betta = sup(u(i)/d(i) - mag(b(i)));  
        x(i) = (b(i) + mradius(0,betta))/(dA(i) + mradius(0,alpha));
    end
end  
  
////////////////////////////////////////////////////////////////////////////////
//
//  запускается интервальный итерационый метод Гаусса-Зейделя, 
//  причём брус, полученный процедурой Хансена-Блика-Рона, 
//  берётся начальным внешним приближением для множества решений 
  
y = x;
  
for k = 1:IterLim;
    for i = 1:n;
        s = b(i);
        for j = 1:n;
            if i~=j then
                s = s - A(i,j)*y(j);
            end
        end
        s = intersection( s/A(i,i), y(i) ); 
        y(i) = s;
    end
end
  
endfunction
  
////////////////////////////////////////////////////////////////////////////////