function [X,N,info] = lineareq(A,B,method,tol)
% PURPOSE: solves linear equation AX + B = 0
%
% if the exists any solution to this equation, then
% X = Y + N PSI
% for any matrix PSI with appropriate size.
%
% ---------------------------------------------------
% USAGE: [X,N,info] = lineareq(A,B)
% where:
% A,B matrices, A need not be square
% method a method used (optional)
% 1 - qr decomposition with pivoting
% 2 - svd decomposition
% tol uses the tolerance tol when calculating
% null subspaces (optional), on default
% tol = max(size(A)') * norm(A,p) * eps,
% where p=1 in case of method 0,1 and p=2
% in case of method 2
%
% X a matrix representing a basic solution
% N a matrix determining a space of all
% solutions
% info return 1 is solution exists and zero
% otherwise. If number of output
% variables is lower than tree, then
% an error is thrown if there are no
% solutions.
%
% COMMENTS:
%
% Copyright (c) Pawel Kowal (2006)
% All rights reserved
% LREM_SOLVE toolbox is available free for noncommercial academic use only.
% pkowal3@sgh.waw.pl
if nargout==3
throw_error = false;
else
throw_error = true;
end
info = 1;
if size(A,1)==0
if size(B,1)==0
X = zeros(size(A,2),size(B,2));
N = eye(size(A,2));
return;
else
error('Inner matrix dimensions must agree.');
end
end
if size(A,2)==0
if size(B,1)==size(A,1)
if norm(B,1)>0
if throw_error
error('there are no solutions');
else
X=[];
N=[];
info = 0;
return;
end
end
X = zeros(0,size(B,2));
N = zeros(0,0);
return;
else
error('Inner matrix dimensions must agree.');
end
end
if nargin<4
tol = 10*max(size(A)') * norm(A,1)* eps;
end
n = size(A,1);
switch method
case 2
[U,S,V] = svd(A);
if size(A,1)>1 && size(A,2)>1
s = diag(S);
else
s = S;
end
tol = 10*max(size(A)') * eps(max(s));
q = sum(s>tol);
B = U'*B;
if rank(B(q+1:end,:),tol)>0
if throw_error
error('there are no solutions');
else
X=[];
N=[];
info = 0;
return;
end
end
B = B(1:q,:);
N = V(:,q+1:end);
V = V(:,1:q);
X = -V*diag(s(1:q).^-1)*B;
case 1
[A,U,k] = putv(A,1,tol);
n = size(A,2);
if k>0
if norm(U(:,k+1:end)'*B,1)>tol
if throw_error
error('there are no solutions');
else
X=[];
N=[];
info = 0;
return;
end
end
B = U(:,1:k)'*B;
end
if size(A,1)==size(A,2)
X = -A\B;
if issparse(A)
N = sparse(n,0);
else
N = zeros(n,0);
end
return;
end;
[A,U,k] = putv(A',1,tol);
X = -U(:,1:k)*(A'\B);
N = U(:,k+1:end);
end