LINEAR SYSTEMS SOLVER (for small Systems)

Solves many small linear systems in a vectorized way. It avoids to call the \ command inside a loop.

Description

The description has been moved to:

http://www.advancedmcode.org/linear-systems-solver-for-small-systems.html

• MATLAB release: MATLAB 7.6 (R2008a)
• Other requirements: Should work on all platforms
• Other Files: license.txt, ; SolverNxN.m, ; TestSolverNxN.m

[John D'Errico]

Pretty good. I like the help, example, etc.

For those interested, the sparse block diagonal matrix is built in a reasonable way. While it uses a loop the could be replaced with a call to ndgrid (see my own blktridiag for a mildly different way to build such a system) the loop here is a simple one that should not be costly.

I also like that the code is batched, so the matrices are limited in size, even for huge sets of problems to solve.

The only thing I'd have done differently would be to supply the matrix systems as a 3-d array, rather than concatenating them all horizontally. Make it a PxQxN matrix, where each plane of the array is a separate system to solve. This way if P>Q, then the linear solves will be least squares problems, all solved in a batch.

[Luigi Giaccari]

Hi John thanks for the comment,

Honestly I didn't know the ndgrid command. More honestly, now that I know, I don't realize how to use it.

I thought about the 3D arrays as input but I personally prefer the 2D concatenation. Anyway in the next revision I will add this feature, since I thing the assembling indexes stays the same and it wont take a long time to code it.

I am not so confortable of least squares problems, I have just used lsqlin a few times. But I thing that this kind of solvers are slower that the ones for linear system, so the advantages of a vectorization will be minor. The for loop and function call time will be hidden by the solver time. As soon as I have a little free time I will run a few speed test and decide if add this new feature.

Thank you again

[John D'Errico]

I have not tested this for speed compared to the looped solution, but it is elegant. It is what I use in blktridiag.

%% Sparse block diagonal matrix, with no explicit loops
% Consider this matrix:
m = 4;
nblocks = 5;
A = rand(m,m*nblocks);

% We wish to generate a sparse block diagonal array,
% with five 4x4 blocks
[ind1,ind2,ind3] = ndgrid(0:m-1,0:m-1,0:nblocks-1);
rind = 1+ind1(:)+m*ind3(:);
cind = 1+ind2(:)+m*ind3(:);

% build the final array all in one call to sparse
Abd = sparse(rind,cind,A(:),nblocks*m,nblocks*m);

% Use spy to see the pattern is correct
spy(Abd)

%% Now do the same thing, but with 4x3 blocks
% Consider this matrix:
p = 4;
q = 3;
nblocks = 5;
A = rand(p,q,nblocks);

% use ndgrid
[ind1,ind2,ind3] = ndgrid(0:p-1,0:q-1,0:nblocks-1);
rind = 1+ind1(:)+p*ind3(:);
cind = 1+ind2(:)+q*ind3(:);

% build the final array all in one call to sparse
Abd = sparse(rind,cind,A(:),nblocks*p,nblocks*q);

% Use spy to see the pattern is correct
spy(Abd)

[Luigi Giaccari]

Hi john,

Your solution seem to be about 5% faster, but the best is it more flexible about the block size,
probably in the next days I will update code with the modification you suggested.

Thank you very much

[Marika]

Hello
I whant to know how to find an aigen vectors which correspond to known eigenvalue. For excample, how to solve the problem: M = [0 0 .5 .5; 0 0 .5 .5; 1 0 0 0; 0 1 0 0] and M'*P=P, where P is the vector, P=[P1; P2; P3; P4] and sum of P's elements equals 1.
Thanks for help.