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 rv
load a %# Nx3 matrix [x1 y1 z1;x2 y2 z2;... ].
%# vectors construct a sphere
L=(301:500)1e-9; K=2pi./L; %# 1x200 array
%# some constants =========================================================
I=eye(3);
e0=1;
[npoints,ndims]=size(rv);
N=npoints;
d0=(4pi/(3N))^(1/3)*5e-9;
V=Nd0^3; aeq=(3V/(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(1iK(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

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).

StackExchange Code Review Q#2413, answer score: 3