Use parallel computing to apply a function to array

Hello everyone,

I am currently building a radiative transfer code to calculate synthetic spectra. The way this works is by having a big table containing the parameter of each line, for example line center nu0, intensity S and lineshape width Width

LineParams =

nu0 S Width ...

--------- ------- -----------

2356 5.3e-20 3.2

2359 6.7e-19 7.2

...

I have a Voigt profile function VoigtProfile which takes as input the parameters of the lineshape from LineParams and an array of frequencies as the x axis nu, returning the spectral line.

For now, this is achieved with a for loop where a finite range of the x axis is taken to compute the lineshape

% Creating the parameter table. In the code comes from scientific databases
% and goes through a series of checks, cleanups and parameter additions
nu0 = [2356 2359];
S = [0.53, 6.7] * 10^(-19);
Width = [3.2 7.2];

LineParams = table(nu0, S, Width);

% Initializing the x-axis and spectrum array
nu = 2200:0.1:2400;
Spectrum = zeros(length(nu), 1)

% Computing the lineshape
for i = 1:length(LineParams.S)
    TempWidth = LineParams.Width(i)
    Tempnu0 = LineParams.nu0(i)
    TempS = LineParams.S(i)
    nuRangeLogic = (nu >= Tempnu0 - 3*TempWidth & nu <= Tempnu0 + 3*TempWidth) % Selecting a range 3 times the width of the lineshape
    Spectrum(nuRangeLogic) = Spectrum(nuRangeLogic) + VoigtProfile(TempS, Tempnu0, Tempwidth, nu(nuRangeLogic)); % Build the spectrum along each section of the nu array
end

function Out = VoigtProfile(S, nu0, Width, X) % For the sake of simplicity, this example actually creates a gaussian profile, but the idea remains the same
    Out = S .* normpdf(X, nu0, Width);
end

The great amount of lines and the high point density required for the spectrum make this implementation very slow. One idea was to implement Parallel Computing to split the load into different workers, hence reducing the amount of time required.

However, is not very clear to me how to proceed. I followed some suggestions from the Matlab editor and read some documentation, resulting in the following implementation

parfor i = 1:length(LineParams.S)
    TempWidth = LineParams.Width(i)
    Tempnu0 = LineParams.nu0(i)
    TempS = LineParams.S(i)
    Spectrum = Spectrum + VoigtProfile(TempS, Tempnu0, Tempwidth, nu);
end

Because of how parfor works, I had to ditch the logical indexing, so the lineshape must be computed on the whole array, which may be even more time-consuming (~100.000 points)!

Is there any better and more orthodox implementation for this case?

Thank you all for any help.

Answer 1

[Harald]

Hi,

you can introduce a temporary variable SpectrumAdd and add this to Spectrum:

Spectrum = zeros(size(nu)); % modified to have same orientation as nu
parfor i = 1:length(LineParams.S)
    SpectrumAdd = zeros(1, length(nu)); % initialize temporary variable with zeros
    TempWidth = LineParams.Width(i);
    Tempnu0 = LineParams.nu0(i);
    TempS = LineParams.S(i);
    nuRangeLogic = (nu >= Tempnu0 - 3*TempWidth & nu <= Tempnu0 + 3*TempWidth); 
    SpectrumAdd(nuRangeLogic) = VoigtProfile(TempS, Tempnu0, TempWidth, nu(nuRangeLogic)); % Use logical indexing to modify part of the temporary variable
    Spectrum = Spectrum + SpectrumAdd; % add the full variable
end

Best wishes,

Harald

[Matt J]

Possibly faster?

[Spectrum,SpectrumAdd] = deal(zeros(size(nu)));

parfor i = 1:length(LineParams.S)
    SpectrumAdd(:)=0;             % <----- use pre-allocated broadcast variable
    TempWidth = LineParams.Width(i);
    Tempnu0 = LineParams.nu0(i);
    TempS = LineParams.S(i);
    nuRangeLogic = abs(nu-Tempnu0) <= 3*TempWidth; % <----- fewer ops
    
    SpectrumAdd(nuRangeLogic) = VoigtProfile(TempS, Tempnu0, TempWidth, nu(nuRangeLogic)); 
    Spectrum = Spectrum + SpectrumAdd;
end

[Francesco]

Thank you all for your help!

This implementation works nicely!

Answer 2

[Walter Roberson]

You are writing to overlapping windows. parfor can only write to non-overlapping windows.

The best you can do is to use parfor to record the start and end of each nuRangeLogic along with the calculated VoigtProfile for that range (as a cell array entry). Then in a serial post processing, accumulate the appropriate portions.

It is not at all clear that this will be any faster. The code you are using is all vectorized and is all eligible for automatic splitting between cores.

The array sizes are potentially too small for automatic splitting, in which case only a single core would be used, so potentially it would be worth using parfor.