Set of linear equations

Question

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?

Answer 1

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