Присвоение нескольким подматрицам матрицы одновременно. Возможна оптимизация по векторизованным показателям - PullRequest
Новая
       10

Присвоение нескольким подматрицам матрицы одновременно. Возможна оптимизация по векторизованным показателям

2 голосов
/ 27 марта 2011

Есть ли умный способ векторизации цикла for, который присваивает элементы подматрицам матрицы?
Изначально у меня было два цикла for: 

U=zeros(6*(M-2),M-2);
for k=2:M-3  
    i=(k-1)*6+1; 
    for j=2:M-3
        U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
    end
end

Затем я векторизовал внутренний циклтакой, что код теперь читается как 

U=zeros(6*(M-2),M-2);
j=2:M-2;
for k=2:M-3
    i=(k-1)*6+1;
    U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);
end

Это уменьшило время моего процессора более чем на 90%, поэтому я подумал, смогу ли я сделать то же самое с внешним циклом, но это кажется немного сложным,Я назначаю (6x1) -матрицы в матрице U.Я попытался 

U=zeros(6*(M-2),M-2);
k=2:M-3;
i=(k-1)*6+1;
j=2:M-2;
U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

, но это не удалось, так как i: i + 5 вынимает только первые 6 индексов, которые я хочу.

Я также пытался использовать функцию reshape () для преобразования матрицы в вектор, но все еще кажется трудным назначать несколько блоков элементов одновременно.Всего в коде три таких цикла for, поэтому я думаю, что альтернативная оптимизация заключается в их параллельном распараллеливании.Однако без доступа к параллельному набору инструментов векторизация представляется мне хорошим решением, если это возможно.

Код является частью подпрограммы в численном методе конечных разностей для решения системы из 6 уравнений на сетке,поэтому этот вопрос может быть актуален для всех, кто работает с матричными вычислениями в системах уравнений, в частности, PDEs .Будем весьма благодарны за предложения по оптимизации кода!

Ответы [ 2 ]

0 голосов
/ 27 марта 2011

Чтобы понять, как вы можете написать назначение в одну строку без циклов, это может помочь нарисовать массив temp в виде прямоугольника. Тогда различные слагаемые, которые будут объединены в U, представляют собой не что иное, как под прямоугольники temp (или подсетки, если вы хотите отслеживать отдельные элементы в temp, что приведет к определенному элементу U), которые сдвинуты влево, вправо, сверху, снизу соответственно. 

%# define row, column shifts
rowShift = 6;
colShift = 1;

%# That's how we'd like to shift 
%# U(i:i+5,j)=A*temp(i:i+5,j)+B*temp(i:i+5,j-1)+C*temp(i:i+5,j+1)+
%# D*temp(i-6:i-1,j)+E*temp(i+6:i+11,j);

%# assign U
U = A * temp(rowShift+1 : end-rowShift, colShift+1 : end-colShift) +... 
    B * temp(rowShift+1 : end-rowShift, 1 : end-2*colShift) + ...
    C * temp(rowShift+1 : end-rowShift, 2*colShift+1 : end) + ...
    D * temp(1 : end-2*rowShift, colShift+1 : end-colShift) + ...
    E * temp(2*rowShift+1 : end, colShift+1 : end-colShift);
0 голосов
/ 27 марта 2011

Чтобы выбрать непрямоугольную часть матрицы, вам нужно использовать линейные индексы: в матрице 3x3 A A(3,3)==A(9) и A([1 3 5 7 9]) - это вектор, который не может быть достигнут через строку /метод индексации столбцов.

Функция sub2ind преобразует индексы строк / столбцов в линейные индексы, поэтому вы можете использовать его в виде sub2ind(size(U),i:i+5,j) для получения линейных индексов одного блокаU. Измените ваш цикл, чтобы он выполнял только работу по сбору линейных индексов, и тогда вы можете сказать вне цикла: 

U(ind_U) = A*temp(ind_A) + B*temp(ind_B) ...

Кроме того, в любое время, когда вы имеете дело с FDM или FEM, подумайте,следует использовать разреженные матрицы.

...