function [y,c] = linear(A,sigma_x,x,sigma_b,b)

// Ouput variables initialisation (not found in input variables)
y=[];
c=[];

// Display mode
mode(0);

// Display warning for floating point exception
ieee(1);

// [y,c]=linear(A,sigma_x,x,sigma_b,b)
// Linear equations solver (generalization)
// Given the system of linear equations Ax=b
// where A is an nxn matrix of the coefficients
// x is a n-components column vector of the unknowns
// b is a n-components column vector of the coefficients
// Suppose that x_sigma1, x_sigma2,...,x_sigmad
// and b_sigma(d+1),...b_sigman are GIVEN coefficients, 
// then:  
// [y,c]=linear(A,sigma_x,x,sigma_b,b)    where
// sigma_x is a d-components column vector of the indices of the given x''s
// sigma_b is a (n-d)-components column vector of the indices of the given b''s
// x is a d-components column vector of the actual value of the given x''s
// b is a (n-d)-components column vector of the actual value of the given b''s
// Then [y,c] are the wanted solutions (the values of x,b which were unknown). 
// PS: the usual inversion of a linear system is recovered by
//     [y,c]=linear(A,[],[],[1:n]'',[b1;b2....;bn])
// 
// R.Coisson (Parma University) and G. Vernizzi (Oxford University) 
// January 2002. 

n = max(size(A));
lengthx = max(size(x));lengthb = max(size(b));
if n~=(lengthx+lengthb) then
  error("wrong TOTAL number of given unknowns")
end;

//% let us create the ''given'' column vectors.
X = zeros(n,1);B = zeros(n,1);
X(sigma_x) = x; B(sigma_b) = b;
RHS = B-A*X;

//the unknowns are x1,x2,x3...,xn,b1,b2,...bn
//DELETING the ones which appear in the GIVEN LIST.
X_unkn = 1:max(size(X));X_unkn(sigma_x) = [];
B_unkn = 1:max(size(B));B_unkn(sigma_b) = [];
M = A(:,X_unkn);
for k = 1:max(size(B_unkn))
  M(B_unkn(k),lengthb+k) = -1;
end;

//The system we have to solve is M*[y;c]=RHS
solution = inv(M)*RHS;
y = solution(1:lengthb);
c = solution(lengthb+1:n);

endfunction
