function x=lin_solv(lufcn,A,b,t);

% __________________________________________________________

% lin_solv.m  solves  Ax=b  using a LU method defined by 

% 'lufcn', and arithmetic of t decimal digits

% ----------------------------------------------------------

% USAGE:  x = lin_solv('lufcn', A, b, t);

%

% ==OR==  x = lin_solv('lufcn', A, b);

%

%where for t > 0, use t decimal digits,

%      for t = 0, use full 16 decimal digits, and

%      `lufcn' must be one of lu2, lu3, lu4.

%For usage 2, t=8 is taken as a default.

%

% See lu2, lu3, lu4, chop

 

if nargin < 3, help lin_solv, return, end

if nargin == 3, t = 8; end 

if t == 0, t = 16; end  



[L,U,P]=feval(lufcn,A,t); 



%% y=forward(L,(P*b),t);

n = length(b); y=P*b; 

for k = 1:n-1;     

    y(k+1:n) = chop(y(k+1:n) - chop(y(k)*L(k+1:n,k),t), t); 

end



%% x=backward(U,y,t);

x=y;

for k = n:-1:1;

    x(k)=chop(x(k)/U(k,k),t);

    x(1:k-1) = chop(x(1:k-1) - chop(x(k)*U(1:k-1,k),t), t);

end

