function z = linpps0(A,b,nu)
% 
%   The simplest PPS-algorithm, implementing parameter partitioning method, 
%   for computing outer component-wise estimates of the solution set to 
%   interval linear algebraic systems of equations. 
% 
%   The code linpps0 is designed for estimating the minimal value of the nu-th 
%   component for the point from the solution set. The maximal component-wise 
%   estimates can be found by negating the minimums for the system with the 
%   right-hand side vector taken with 'minus' sign. 
%  
%   As a basic enclosure method, linpps0 uses interval Gauss method (function 
%   'IGauss'). For simplicity, only the interval matrix is bisected, which 
%   may cause 'non-optimality' of the estimates computed after the natural 
%   termination of the algorithm. 
% 
%   Input argumetns:
%       A  is the interval matrix of the system under solution, 
%       b  is the interval right-hand side vector of the system, 
%       nu is the number of the component along which the solution 
%             set should be estimated. 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  MaxStepNumber = 100;  %   limitation on the number of steps to be executed 
                        %   by the algorithm 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   inputting the data and checking their correctness 
%  
  A = intval(A);
  b = intval(b);
  
  m = size(A,1);
  n = size(A,2);
  if ( m~=n )
      error('The input matrix must be square!')
  end
  
  k = size(b,1);
  l = size(b,2);
  if ( k~=m | l>1 ) 
      error('Right-hand side size does no match to that of the matrix!')
  end 
  
  if ( nu < 1 | nu > n ) 
      error('Estimated component number is invalid!')
  end 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
%   launching the PPS-algorith proper 
%  
x = IGauss(A,b); 
  
NStep = 1; 
WorkList(1) = struct( 'Matrix', A, 'Estim', inf(x(nu)) );
  
while ( NStep < MaxStepNumber ) 
    
    %   searching the leading ILAS in WorkList 
    LL = size(WorkList,2); 
    nuEst = realmax; 
    for k = 1:LL; 
        p = WorkList(k).Estim;
        if p < nuEst
            nuEst = p;
            Leading = k; 
        end
    end
    
    display(NStep), display(WorkList(Leading).Matrix) 
    
    %   bisect the matrix of the leading ILAS to descendants Matr1 è Matr2  
    Matr1 = WorkList(Leading).Matrix;   Matr2 = Matr1; 
    maxrad = -realmax; 
    for i = 1:n; 
        for j= 1:n;
            mm = rad(Matr1(i,j));
            if maxrad < mm
                maxrad = mm;
                mi = i; 
                mj = j;
            end 
        end
    end 
    
    if maxrad == 0
        break
    end
    
    s = Matr1(mi,mj); 
    Matr1(mi,mj) = intval(inf(s)); 
    Matr2(mi,mj) = intval(sup(s)); 
    
    %   compute enclosures of the solution sets for systems-descendants 
    x1 = IGauss(Matr1,b);   x2 = IGauss(Matr2,b);  
    
    %   construct the records-descendants and put them into WorkList 
    Desc1 = struct( 'Matrix', Matr1, 'Estim', inf(x1(nu)) );
    Desc2 = struct( 'Matrix', Matr2, 'Estim', inf(x2(nu)) );
    WorkList = [ WorkList  Desc1 ];
    WorkList = [ WorkList  Desc2 ];
    
    WorkList(Leading) = []; %   deleting the former leading record from WorkList 
    NStep = NStep + 1; 
  
end
  
z = nuEst; 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
    function x = IGauss(A,b)
    %  interval Gauss method with partial pivoting 
  
        C = A;  d = b;  
  
        for j = 1:n-1;
            [maxmig,index] = max(mig(C(j:n,j)));
            if maxmig <= 0 
                error('All possible partial pivots contain zero!')
            end
            k = index + j - 1;
            if k ~= j       %   transpose the rows if necessary 
                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);     %   backward substitution 
    
        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
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  