Performing inline 3D matrix calculations
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
All:
I have the following code:
% f, tau, jonesx, jonesy are (1, numParams) vectors
% sxyt = 1e6; numFrames = 10; numParams = 18;
A = single(zeros(sxyt, numFrames));
B = single(zeros(sxyt, numFrames));
osc = single(zeros(sxyt, numFrames, numParams));
mag = single(zeros(sxyt, numFrames, numParams));
DxDyOut = single(zeros(sxyt, numFrames, 2));
for k = 1:numParams
osc(:,:,k) = exp((2.*pi.*f(k).*t(:,:) + randomPhase(k)).*1i);
mag(:,:,k) = a(k) .* (k1(k) + k2(k).*exp(-t(:,:)/(tau(k) + 0.0001)));
A = A + real( jonesx(k) .* osc(:,:,k) .* mag(:,:,k) );
B = B + real( jonesy(k) .* osc(:,:,k) .* mag(:,:,k) );
end
DxDyOut(:,:,1) = A;
DxDyOut(:,:,2) = B;
I would like to perform these operations inline rather than in the loop. This code gets called many times and for each time it takes ~7.5s for first equation and 2.5s for 2nd thru 4th equations (for sxyt = 1e6; numFrames = 10; numParams = 18), so roughly 400 ms and 140 ms per execution of a line.
Is there any way one can eliminate the loop? I figure it will be about 10x faster inline.
I looked at mtimes and mtimesx (in user submissions) but can't quite figure out how either may be used for this. I thought I could use mtimesx submission as follows:
- Convert the (1, numParams) vectors to (numFrames,numFrames,numParams) size - each value in the (:,:,i) would be the same as the (1,i) element.
- Convert the (sxyt, numFrames) vectors to (sxyt,numFrames, numParams) size
Any help or suggestion on how to use mtimesx or any other method is much appreciated.
2 commentaires
Walter Roberson
le 22 Mai 2015
I would still like to know the speed difference between doing this as single precision compared to doing it in double precision.
Réponses (1)
Walter Roberson
le 22 Mai 2015
Pull values out of the loop when practical. For example 2.*pi.*f(k) can be calculated ahead of the loop as can k2(k)/(tau(k)+0.0001). You can also factor out a(k)
f_2pi = single(2.*pi.*f);
k2_tau = single(f2 ./ (tau + 0.0001));
a_k1 = a .* k1;
for k = 1:numParams
osc(:,:,k) = exp(f2pi(k).*t(:,:) + randomPhase(k)).*1i);
mag(:,:,k) = a_k1(k) + a(k) .* k2_tau(k).*exp(-t(:,:));
A = A + real( jonesx(k) .* osc(:,:,k) .* mag(:,:,k) );
B = B + real( jonesy(k) .* osc(:,:,k) .* mag(:,:,k) );
end
and then clearly a(k) can be factored into k2_tau(k):
f_2pi = single(2.*pi.*f);
a_k2_tau = single(a .* f2 ./ (tau + 0.0001));
a_k1 = single(a .* k1);
for k = 1:numParams
osc(:,:,k) = exp(f2pi(k).*t(:,:) + randomPhase(k)).*1i);
mag(:,:,k) = a_k1(k) + a_k2_tau(k).*exp(-t(:,:));
A = A + real( jonesx(k) .* osc(:,:,k) .* mag(:,:,k) );
B = B + real( jonesy(k) .* osc(:,:,k) .* mag(:,:,k) );
end
Also experiment,
for k = 1:numParams
osck = single(exp(f2pi(k).*t + randomPhase(k)).*1i));
magk = single(a_k1(k) + a_k2_tau(k).*exp(-t));
osc_magk = osck .* magk;
A = A + real( jonesx(k) .* osc_magk );
B = B + real( jonesy(k) .* osc_magk );
osc(:,:,k) = osck;
magk(:,:,k) = magk;
end
The next question would be whether jonesx and jonesy consist entirely of real values. If they do, then you can move the real():
for k = 1:numParams
osck = single(exp(f2pi(k).*t + randomPhase(k)).*1i));
magk = single(a_k1(k) + a_k2_tau(k).*exp(-t));
osc_magk = single(real(osck .* magk));
A = A + jonesx(k) .* osc_magk;
B = B + jonesy(k) .* osc_magk;
osc(:,:,k) = osck;
magk(:,:,k) = magk;
end
Possibly some of the single() calls can be omitted.
1 commentaire
Voir également
Catégories
En savoir plus sur Function Creation dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!