solvernxn.m

fastsolving linear equation

function x=SolverNxN(Amxn,b)
% SolverNxN is designed to solve many small linear system in a vectorized
%  way in order to improve time perfomaces. It is in practice a faster
%  alternative to the matlab commands:
% 
% for i=1:number_of_sytems
%     x=A\b;
% end
% 
% The Matlab looping speed can visibly slow down this process so most of 
%  the tme is spent for the for loop and not to solve the system. This is
%  particulary true when A is a small matrix. In this case SolverNxN
%  can be a usefull tool.
% The speed improvement is more visibile the greater number are the systems
%  and the smaller is their dimension.
% The algorithm works assembling all the small elementary systems into a
%  sparse matrix so the \ command is called just one time.
%VERY IMPORTANT:
% Because of algorithm structure is very important that no system is
%  singular. One singular system can ruin the whole solution. So use this
%  routine only if you are sure all the systems are solvable.
%
%  
% 
% Let's say we want to solve Nsyst m x m linear systems. 
% 
% The function call is:
% 
%      x=SolverNxN(Amxn,b);
% 
% INPUT:
% 
% Amxn: is a matrix of dimensions [m ,(m x Nsyst)], so it is just an
%       horizontal concatenation of the small elementary systems.
% 
% b:  is a [(m x Nsyst),1] vector formed by the vertical
%       concatenation of the elematry system vectors.
% 
% 
% OUTPUT:
% 
% x:  is a [(m x Nsyst),1] vector formed by the vertical
%        concatenation of the elematry system solutions.
%   
% 
% Here is offered an example:
% 
% 
% m=3;%we want to solve mxm sytems
% Nsyst=100000;%number of systems  
% 
% Amxn=rand(m,Nsyst*m);%elements of nsyst mxm systems
% b=rand(Nsyst*m,1);%system vector
% 
% %launch the function 
% tic
% x=SolverNxN(Amxn,b);
% toc
% 
% 
% 
% For question,suggestion,bugs report: [e-mail]
% It is brand new so please report problems to my e-mail
% 
% Author: Giaccari Luigi
% Created: 12/02/2009
% Last Update: 12/02/2009


%errors check
if nargin~=2
    error('Two inputs required')
end

%check the matrix dimension

[m,n]=size(Amxn);%We have to solve n/m  mxm systems



if n~=floor(n/m)*m
    error('The number of columns is not multiple of the number of rows, check the input matrix')
end

%check the vector dimension
[m2,n2]=size(b);
if m2~=n
    error('The system vector dimension dismatch the matrix dimension')
end
if n2~=1
    error('Second input must be a vector')
end
clear m2 n2

x=zeros(n,1);


%Solve the system in small part each times, this avoids help memory errors
%Briefly the global system is divded is a series of smaller packages
numberofcut=floor((n*m)/(100001))+1;%100001 is the cutsize
indexsize=floor(n/numberofcut)*m;%we solve indexsize each loop time

%solve from i1 to i2
i1=1;
i2=indexsize;

%call the solver for each package
while i2<n

    x(i1:i2)=Solver(Amxn(:,i1:i2),b(i1:i2));
    
    %increase the package indexes
    i1=i1+indexsize;
    i2=i2+indexsize; 
end

%the last package
x(i1:end)=Solver(Amxn(:,i1:end),b(i1:end));


end

 
function x=Solver(Amxn,b)

 
%get system size
[m,n]=size(Amxn);


%Get the coordinates of the elements and place them in the index arrays
%this is usefull to vectorize the operation

index=zeros(n*m,2);%can not be converted to integer type to save memory since
                       %the sparse command requires doubles
                                              
m2=m*m;
k=1;
%at each loop we assemble the element(k,ceil(i/m)) for all of the
%elementary systems
 for i=1:m2
     
index(i:m2:end,1)=k:m:n;
index(i:m2:end,2)=ceil(i/m):m:n;

if k>=m
    k=1;
else
    k=k+1;
end

 end

%Assembly the global system
A=sparse(index(:,1),index(:,2),Amxn(:),n,n);
clear index


% to check if well assembled
% figure(1)
% hold on
% title('Assembled system','fontsize',14);
% spy(A)

x=A\b;%solve the system

end