function z = linpps1(A,b,nu)
% 
%   Простейший PPS-алгоритм - метод дробления параметров - для вычисления 
%   внешних покоординатных оценок объединённого множества решений интервальных 
%   систем линейных алгебраических уравнений (ИСЛАУ). 
% 
%   Программа linpps1 предназначена для вычисления минимального значения 
%   nu-ой компоненты точек из множества решений. Максимальная покомпонентная 
%   оценка может быть найдена как минус минимум той же величины для системы 
%   с противоположной правой частью. 
%  
%   В качестве базового метода внешнего оценивания в программе используется
%   интервальный метод Гаусса (функция IGauss). 
% 
%   Входные параметры программы: 
%       A   - квадратная интервальная матрица решаемой ИСЛАУ,
%       b   - интервальный вектор правой части ИСЛАУ, 
%       nu  - номер компоненты, вдоль которой оценивается множество решений. 
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  MaxStepNumber = 100;  %   ограничение на количество итераций PPS-алгоритма 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   вводим данные и проверяем их корректность: квадратная ли матрица, 
%   соответствуют ли друг другу размерности матрицы и правой части? 
%  
  A = infsup(A);
  b = infsup(b);
  
  m = size(A,1);
  n = size(A,2);
  if ( m~=n )
      error('Матрица решаемой системы должны быть квадратной!')
  end
  
  k = size(b,1);
  l = size(b,2);
  if ( k~=m | l>1 ) 
      error('Размеры правой части не соответствуют размерам матрицы!')
  end 
  
  if ( nu < 1 | nu > n ) 
      error('Номер оцениваемой компоненты некорректен!')
  end 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
%   запускаем собственно PPS-алгоритм,  
%   и сначала инициализируем рабочий список WorkList 
  
x = IGauss(A,b);   % находим брус начальной вншней оценки 
  
NStep = 1; 
WorkList(1) = struct( 'Matrix', A, 'RHS', b, 'Estim', inf(x(nu)) );
  
%   запускаем основной цикл уточнения оценки множества решений 
  
while ( NStep < MaxStepNumber ) 
    
    %   находим в рабочем списке ведущую ИСЛАУ,  
    %   она получает номер k 
    LL = size(WorkList,2); 
    nuEst = realmax; 
    for k = 1:LL; 
        p = WorkList(k).Estim;
        if p < nuEst
            nuEst = p;
            Lead = k; 
        end
    end
    
    display(NStep), 
    display(WorkList(Lead).Matrix), 
    display(WorkList(Lead).RHS) 
    
    A1 = WorkList(Lead).Matrix;     A2 = A1; 
    b1 = WorkList(Lead).RHS;        b2 = b1; 
    
    %   определяем в матрице и правой части самый широкий интервальный элемент 
    Amaxrad = -realmax; 
    for i = 1:n; 
        for j= 1:n;
            mm = rad(A1(i,j));
            if Amaxrad < mm
                Amaxrad = mm;
                mAi = i; 
                mAj = j;
            end 
        end
    end 
    
    bmaxrad = -realmax; 
    for i = 1:n; 
        mm = rad(b1(i));
        if bmaxrad < mm
            bmaxrad = mm;
            mbi = i; 
        end 
    end 
    
    %   если все элементы имеют нулевую ширину, 
    %   то останавливаем выполнение алгоритма 
    if Amaxrad == 0 && bmaxrad == 0 
        break
    end
    
    %   дробим матрицу ведущей ИСЛАУ или её правую часть на потомки 
    %   по самому широкому интервальному элементу 
    if Amaxrad >= bmaxrad 
        %   дробим матрицу, а правую часть оставляем неизменной 
        s = A1(mAi,mAj); 
        A1(mAi,mAj) = infsup(inf(s)); 
        A2(mAi,mAj) = infsup(sup(s)); 
    else 
        %   дробим правую часть, матрицу оставляем неизменной 
        s = b1(mbi); 
        b1(mbi) = infsup(inf(s)); 
        b2(mbi) = infsup(sup(s)); 
    end 
    
    %   находим внешние интервальные оценки множеств решений 
    %   для потомков ведущей ИСЛАУ 
    x1 = IGauss(A1,b1);   x2 = IGauss(A2,b2);  
    
    %   формируем записи-потомки и помещаем их в рабочий список 
    Desc1 = struct( 'Matrix', A1, 'RHS', b1, 'Estim', inf(x1(nu)) );
    Desc2 = struct( 'Matrix', A2, 'RHS', b2, 'Estim', inf(x2(nu)) );
    WorkList = [ WorkList  Desc1 ];
    WorkList = [ WorkList  Desc2 ];
    
    %   удаляем бывшую ведущую ИСЛАУ из рабочего списка 
    WorkList(Lead) = [];
      
    NStep = NStep + 1; 
   
end
  
%   в качестве результата работы алгоритма 
%   берём последнюю ведущую оценку   
z = nuEst; 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
    function x = IGauss(A,b)
    %  интервальный метод Гаусса с частичным выбором главного элемента
  
        C = A;  d = b;  
  
        for j = 1:n-1;
            %   находим ведущий элемент - имеющий наибольшую мигнитуду - 
            %   в поддиагональной части текущего столбца 
            [maxmig,index] = max(mig(C(j:n,j)));
            if maxmig <= 0 
                error('Все поддиагональные элементы столбца содержат нули!')
            end
            k = index + j - 1;
            if k ~= j  then     %   при необходимости переставляем строки  
                C([j,k],j:n) = C([k,j],j:n);
                d([j,k]) = d([k,j]);
            end
            %    выполняем преобразования приямого хода метода Гаусса 
            for i = j+1:n;
                r = C(i,j)/C(j,j);
                C(i,j+1:n) = C(i,j+1:n) - r*C(j,j+1:n);
                d(i) = d(i) - r*d(j);
            end
        end
    
        %   обратная подстановка 
        %   для нахождения интервальных оценок компонент 
        x(n) = d(n)/C(n,n); 
    
        for i = (n-1):(-1):1;
            r = d(i);
            for j = i+1:n;
                r = r - C(i,j)*x(j);
            end
            x(i) = r/C(i,i);
        end
  
    end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
