Assegnazione a più sotto-matrici di una matrice contemporaneamente.Possibile ottimizzazione per indici vettorializzati

Question

C'è un modo intelligente per vettorizzare un ciclo for che assegna elementi alle sottomatrici di una matrice?
Inizialmente, avevo due for-loop:

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

Quindi ho vettorizzato il ciclo interno, in modo tale che il codice ora legga

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

Questo ha ridotto il mio tempo di CPU di oltre il 90%, quindi mi sono chiesto se potessi fare lo stesso con il ciclo esterno, ma sembra un po ' complicato, dal momento che assegno a (6x1)-matrici all'interno della matrice U.Ci ho provato

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

ma questo fallisce, dal momento che i: i + 5 estrae solo i primi 6 indici che voglio.

Ho anche provato a usare la funzione reshape () per convertire la matrice in un vettore, ma sembra ancora difficile assegnare a più blocchi di elementi contemporaneamente.Ci sono in totale tre di questi for-loop nel codice, quindi immagino che un'ottimizzazione alternativa sia quella di parallelizzarli in qualche modo.Tuttavia, senza accesso alla casella degli strumenti parallela, la vettorializzazione mi sembra una buona soluzione se possibile.

Il codice fa parte di una subroutine in un metodo numerico a differenza finita per risolvere un sistema di 6 equazioni su una griglia, quindi questa domanda potrebbe essere rilevante per chiunque lavori con calcoli matriciali su sistemi di equazioni, in particolare PDE.Suggerimenti per ottimizzare il codice sarebbero molto apprezzati!

Answer 1

Per capire come è possibile scrivere l'assegnazione in una riga senza loop, può aiutare a disegnare l'array temp come un rettangolo.Quindi, le diverse sommatorie che si uniranno a U non sono altro che sotto-rettangoli di temp (o sotto-griglie, se si desidera tenere traccia dei singoli elementi in temp ciò si tradurrà in un elemento specifico di U) che vengono spostati a sinistra, destra, in alto, in basso, rispettivamente.

%# 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);

Answer 2

Per selezionare una porzione non rettangolare di una matrice, è necessario utilizzare gli indici lineari: in una matrice 3x3 A, A(3,3)==A(9) e A([1 3 5 7 9]) è un vettore che non può essere ottenuto tramite il metodo di indicizzazione della riga / colonna.

the sub2ind Function converte gli indici di riga / colonna in indici lineari,Quindi puoi usarlo nel modulo sub2ind(size(U),i:i+5,j) per ottenere gli indici lineari di un blocco di U. cambia il tuo loop per fare solo il lavoro di raccolta degli indici lineari, e quindi puoi dire al di fuori del loop:

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

.

Inoltre, ogni volta che hai a che fare con FDM o FEM, considerare se dovresti usare le matrici sparse.