Набор линейных уравнений


У меня \$ н^2 \$ матриц. Каждая из них является \ 3 $ \3 раза \$ матрица. Я хочу конкатенации их \$ 3Н \раз 3Н \$ матрица. Я знаю, что большая матрица симметрична. Оценка всей \ 3 $ \3 раза \$ матрицы отнимает много времени, поэтому я хочу ускорить мою программу.

У вас есть какие-либо пожелания, чтобы ускорить его?

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?

Пожалуйста получите.коврик : http://dl.dropbox.com/u/21031944/Stack/a.mat

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

Или sth как это (для симметрии часть):

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    

Я не знаю четкой программы?



399
4
задан 15 мая 2011 в 04:05 Источник Поделиться
Комментарии
1 ответ

По моему опыту с методом конечных элементов проблемы, я обнаружил, что создают большие матрицы в легко ушло столько же времени, если не больше, чем вычисления решения х = А/B. Матлаб делает хорошую работу с тем, что обратная матрица.

В (старых) Фортран создание это делается с вложенными циклами. В MATLAB эту итерации происходит очень медленно. Один фокус в том, чтобы изменить порядок вложенности, поэтому вы переходите на небольшие размеры (в 3х3), и использование матричных операций на крупных (Н). Следующим шагом является, чтобы сделать все с матричных операций (на 4Д массивов при необходимости), но это труднее осмыслить.

Я подозреваю, что ваш а матрица разрежена (много глаз компоненты). В этом случае разреженной матрицы Матлаб операции может быть ваш друг. В принципе можно создать 3 больших векторов, я, Джей, в, котором 2 определить координаты, и третий (ненулевое) значение в каждой. Если повторяются определенной пары координат, соответствующие значения будут добавлены, что в некоторых определениях матрицы удобно. Но это больше, чем я могу здесь объяснить.

https://codereview.stackexchange.com/a/30101/27783
берет простая модель диффузии тепла, используя петли и индексация, и ускоряет его в несколько шагов. Просто векторизации внутренние петли помогает. Лучшей скорости пошел на создание матрицы с помощью глаз (или рис), и делает быстрое умножение матриц (с НП.точка).

3
ответ дан 22 августа 2013 в 01:08 Источник Поделиться