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=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?