3 Parallelizing the Models

 3.1 High-pass Filter
 3.2 Commom-mode Model
 3.3 Astronomical Sky Model
 3.4 Other Models

The parallelization of the iterative map-maker is accomplished by the individual model. Each model is responsible for ensuring that it is fit to all applicable data and that the best-fit model is distributed to all processes. There are three classes of models to consider; those that require:

(1)
all time samples for each individual detector (e.g. high-pass filter)
(2)
all detector samples at each time slice (e.g. common mode)
(3)
some more-complicated division of the data set (e.g. astronomical signal, which could be solved by dividing the samples by region on the sky)

We can choose the division of the data set between nodes to make one of the model types easy to calculate, but not all of them. Because the high-pass filter model performs FFTs on the detector time streams, which we would hard to perform if an individual detector was spread out across nodes, we decide to split the data between nodes so that each process acts on the complete time streams for a number of detectors.

The parallelization of the example models listed above are now described in some detail to illustrate how the parallelization of the different model-types is accomplished.

3.1 High-pass Filter

Figure 3 illustrates how the existing multithreaded map-maker performs the convolutions necessary for the high-pass filter model. Since each detector time-stream can be processed independently, the detectors are divided into Nthread chunks and each thread performs the convolution of each detector with the filter and stores the result directly in a shared-memory array.


pdfpict
Figure 3: Illustration of the calculation of the high-pass filter model on four detectors with eight time-samples each, using two threads. The red outlines indicate how the data are chunked for processing by each thread.


With full detector time streams on each node, parallelizing the high-pass filter model is trivial, as shown in Figure 4. Each node is responsible for Ndet/Nnode detectors which are then further split up for each thread. No communication between processes is required.


pdfpict
Figure 4: Parallel version of the high-pass filter model. Here we have two processes (represented by black outlines), each with two threads. The model is trivially parallelized, as each detector can be processed independently.


3.2 Commom-mode Model

The calculation of the common-mode model is illustrated in Figure 5. The data are now chunked along the detector dimension so that each thread calculates the average detector value for each of a range of time slices. This is shown in the figure as a two-step process: the signal and number of samples are accumulated into arrays, and then the accumulated sums are divided by the number of samples. It is shown this way so that the generalization to the parallel version is more clear.


pdfpict
Figure 5: Illustration of the calculation of the common-mode model on a single node using two threads. The data are chunked along the detector axis so that each thread can calculate the common-mode signal (average value) at each time sample. The calculation of the average is shown in two steps: first the signal and number of detectors are accumulated at each time slice, then the average at each step is calculated. Note that it is simple to extend this to a weighted average by accumulating wdtxdt and wdt instead of xdt and 1.


In the parallel version, a single process is not able to fully calculate the model. Instead, each process accumulates the sums for its set of detectors, then MPI communications are used to accumulate the arrays across all processes. Each process can then perform the final division of the accumulated arrays. The parallel version is shown in Figure 6.


pdfpict
Figure 6: Parallel version of the common-mode model. An extra commumincation step (represented by green arrows) is inserted, where the arrays are accumulated across all nodes. Multithreading is still used for parallelized calculation of the model.


The code to perform this operation might look something like this:

  for (idet=0; idet<ndet; idet++) {
    for (isamp=0; isamp<nsamp; isamp++) {
      signal[isamp] += data[idet*ndet+isamp];
      count[isamp] += 1;
    }
  }
  
  #if USE_MPI
  MPI_Allreduce(MPI_IN_PLACE, signal, nsamp, MPI_DOUBLE, MPI_SUM,
                MPI_COMM_WORLD);
  MPI_Allreduce(MPI_IN_PLACE, count, nsamp, MPI_DOUBLE, MPI_SUM,
                MPI_COMM_WORLD);
  #endif
  
  for (isamp=0; isamp<nsamp; isamp++) {
    signal[isamp] /= count[isamp];
  }

This program is run by each process, where ndet is the number of detectors “owned” by the process, and is not restricted to be the same in each process. The communications are performed by the MPI command MPI_ALLreduce, which performs the exact operation we need: an array (signal and count in the examples above) is combined across processes with a specified operation (here, we specify “sum”) and the result is returned in place to each process. Following the communication step, each process calculates the average for itself—this is quicker than having a single thread calculate the average and then communicate the result. Also note that, without the “#if USE_MPI/#endif” clause, the code is identical to what is calculated by the serial version.

3.3 Astronomical Sky Model

The astronomical sky model rebins the data time streams, after subtraction of noise models, into a map of the sky, the end product of the map-maker. The projection for detector/time sample to map pixel is a complicated function depending on the focal plane layout and telescope scan pattern, and thus the task cannot be easily subdivided into a number of independent tasks in the same way that is done for the filter and common-mode models. Instead, each thread in the single-node version of the map-maker has its own copy of the map arrays and accumulates a chunk of the data samples; the data could be chunked either by detector or time sample. After each thread has finished accumulating its samples, the accumulator arrays are combined and the final division is performed.

The paralellization of this calculation is analogous to the common-mode model; after the threaded accumulator arrays are combined, the arrays are accumulated across nodes, again using MPI_Allreduce, and the divisions are performed by each process.

3.4 Other Models

Models other than the ones discussed above should fall within one of the categories discussed above. Communication commands other than MPI_Allreduce or operators other than MPI_SUM may be needed for more complicated model-fitting.