# Set of linear equations

I have $N^2$ matrices. Each one is a $3 \times 3$ matrix. I want to concatenation them to a $3N \times 3N$ matrix. I know that the big matrix is symmetric. Evaluation of whole $3 \times 3$ matrices are time consuming so I want to speed up my program.

Do you have any suggestion to speed it up?

clc
load a                          %#  Nx3 matrix [x1 y1 z1;x2 y2 z2;... ].
%#  vectors construct a sphere
L=(301:500)*1e-9;  K=2*pi./L;   %# 1x200 array
%# some constants =========================================================
I=eye(3);
e0=1;
[npoints,ndims]=size(rv);
N=npoints;
d0=(4*pi/(3*N))^(1/3)*5e-9;
V=N*d0^3; aeq=(3*V/(4*pi))^(1/3);
E0y=ones(N,1);
E0z=E0y;
Cext=zeros(1,length(L));
Qext=zeros(1,length(L));
A=zeros(3,3,N^2);
% =================================
r=sqrt(sum(rv,2).^2);           %# [Nx1] lengrh of each rv vector
for i=1:N
for j=1:N
dx(i,j)=rv(i,1)-rv(j,1); %# The x component of distances between each 2 point [NxN]
dy(i,j)=rv(i,2)-rv(j,2); %# The y component of distances between each 2 point [NxN]
dz(i,j)=rv(i,3)-rv(j,3); %# The z component of distances between each 2 point [NxN]
end
end
dv=cat(3,dx,dy,dz);  %# d is the distance between each 2 point (a 3D matrix) [NxNx3]
d=sqrt(dx.^2+dy.^2+dz.^2);                     %# length of each dv vector
nx=dv(:,:,1)./d; ny=dv(:,:,2)./d; nz=dv(:,:,3)./d;
n=cat(3,nx,ny,nz);                              %# Unit vectors for points that construct my sphere

for s=1:length(L)
E0x=exp(1i*K(s)*rv(:,1))'; % ' #closing the quote for syntax highlighting
% 1x200 array  in direction of x(in Cartesian coordinate system)
% Main Loop    =================================================
p=1;
for i=1:N
for j=1:N
if i==j
A(:,:,p)=a(s)*eye(3);           %# 3x3 , a is a 1x200 constant array
p=p+1;                          %# p is only a counter
else
A(:,:,p)=-exp(1i*K(s)*d(i,j))/d(i,j)*(-K(s)^2*([nx(i,j);ny(i,j);nz(i,j)]...
*[nx(i,j) ny(i,j) nz(i,j)]-I)+(1/d(i,j)^2-1i*K(s)/d(i,j))...
*(3*[nx(i,j);ny(i,j);nz(i,j)]*[nx(i,j) ny(i,j) nz(i,j)]-I));
p=p+1;
end
end
end

%===============================================================
B = reshape(permute(reshape(A,3,3*N,[]),[2 1 3]),3*N,[]).';    %# From :gnovice
%# concatenation of N^2 3x3 matrixes into a 3Nx3N matrix

for i=1:N
E00(:,i)=[E0x(i) E0y(i) E0z(i)]';
end
b=reshape(E00,3*N,1);
P=inv(B)*b;
Cext(s)=(4*pi*K(s))*imag(b'*P);
Qext(s)=Cext(s)/(pi*aeq^2);
end

Qmax=max(Qext); Qext=Qext/Qmax;
L=L*1e9;
plot(L,Qext,'--');figure(gcf)

%# The B matrix is symmetric(B_ij=B_ji) so I can reduce the computations
%# I think I should write sth like this
% for i=1:N
%        for j=i:N
%            if i==j
%                A(:,:,p)=a(s)*eye(3);           %# 3x3 , a is a 1x200 constant array
%                p=p+1;                          %# p is only a counter
%            else
%            A(:,:,p)=... [just same as above]
%            p=p+1;
%            end
%        end
%    end
% But how concatenate them like befor in order?


Pls get a.mat : http://dl.dropbox.com/u/21031944/Stack/a.mat

     rv.mat:   http://dl.dropbox.com/u/21031944/Stack/rv.mat


Or sth like this (for symmetry part):

for i=1:N
for j=1:N
if i==j
A(:,:,p)=a(s)*eye(3);           %# 3x3 , a is a 1x200 constant array
p=p+1;                          %# p is only a counter
elseif i>j
A(:,:,p)=... [just same as above]
p=p+1;
else
A(:,:,p)=A(:,:,??);
p=p+1;
end
end
end


I don't know is program clear?

• A bigger rv May 15, 2011 at 4:58
• May 15, 2011 at 4:59
• Could you run the profiler to find out which part exactly requires the speedup? Dec 27, 2012 at 14:52

In my experience with finite element problems, I found that creating the large matrix A easily took as long, if not longer, than calculating the solution x = A/b. Matlab does a good job with that matrix inverse.
In (old) FORTRAN creation of A is done with nested loops. In Matlab this iteration is very slow. One trick is to reverse the order of the nesting, so you iterate on small dimensions (the 3x3), and use matrix operations on the large ones (N). A further step is to do everything with matrix operations (on 4D arrays if needed), but that is harder to conceptualize.
I suspect your A matrix is sparse (lots of eye components). In that case the Matlab sparse matrix operations can be your friend. Basically you create 3 large vectors, I, J, V, where 2 identify the coordinates, and the third the (nonzero) value at each. If a specific pair of coordinates is repeated, the corresponding values are added, which in some matrix definitions is convenient. But that is more than I can explain here.
https://codereview.stackexchange.com/a/30101/27783 takes a simple heat diffusion model, using loops and indexing, and speeds it up in several steps. Simply vectorizing the inner loops helps a lot. The best speed came generating a matrix using eye (or diag), and doing fast matrix multiplication (with np.dot).
• If the goal is just solving linear equations, then indeed using ldivide or rdivide seems to be the way to go. Aug 26, 2013 at 12:53