OpenMP Parallelization in SPECTUB

[jchan192]

Hello STIR community,

I have been working on enabling OpenMP parallelization in the SPECT UB module and would like to share my changes, along with my current results. I’d appreciate any feedback on whether my implementation introduces any potential issues.

First, I removed the omp critical in this line to allow parallelization

STIR/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx

Lines 897 to 906 in 42aebde

	#ifdef STIR_OPENMP
	# pragma omp critical(PROJMATRIXBYBINUBONEVIEW)
	#endif
	if (!subset_already_processed[kOS])
	{
	if (!this->keep_all_views_in_cache)
	{
	this->clear_cache();
	subset_already_processed.assign(prj.NOS, false);
	}

Secondly , I declared two new variables , local_wm and local_wmh, at the beginning of
ProjMatrixByBinSPECTUB::compute_one_subset. Then, I replace all the this->wm and this->wmh with local_wm and local_wmh in the ProjMatrixByBinSPECTUB::compute_one_subset. Since the declaration is a bit lengthy, I will show if requested.

I was testing the implementation by running OSMAPOSL, and the computational time does decrease when using 2 and 4 threads. The resulting reconstructed images do show some differences between 1 thread and multi-threads, primarily in regions with higher numerical values, but the differences are relatively small (~1e-5). Here I uploaded a slice of the difference

Do you think the approach to parallelizing SPECTUB seem reasonable?
Thanks in advance,
Jowett

[KrisThielemans]

Thanks for trying this. However, I don't think this can work correctly/efficiently.

Certainly, making wm and wmh local variables in compute_one_subset is a good idea. It's a further step along the path taken in #1169. In fact, it makes me wonder if we need the wm class member at all, which would be really good.

However, by removing the omp critical, the parallelisation will now call compute_one_subset in multiple threads. That means that a huge amount of work is actually repeated in all those threads. It might work, but it's very inefficient. (There might be a problem with some LORs due to a bug I just saw, see #1576).

When @danieldeidda and I discussed multi-threading SPECTUB previously, we thought it'd need either major restruction of the original SPECTUB code (which we're not so eager to do...), or using multi-threading inside the weight calculation. @danieldeidda started #460 for this, but we ran out of time/steam.

Happy to discuss this further!

[jchan192]

thanks for the response!

Since the current approach has each thread handling a different viewgram independently (code below), I don't understand how "a huge amount of work is actually repeated in all those threads" . Even though there are repeated parts like the nested loops in the wm_calculation, i think they are unavoidable because we need to process each viewgram anyway?

STIR/src/recon_buildblock/distributable.cxx

Lines 431 to 443 in 42aebde

	# if _OPENMP < 201107
	# pragma omp for schedule(dynamic)
	# else
	// OpenMP loop over both vs_nums_to_process and tof_pos_num
	# pragma omp for schedule(dynamic) collapse(2)
	# endif
	#endif
	for (int timing_pos_num = min_timing_pos_num; timing_pos_num <= max_timing_pos_num; ++timing_pos_num)
	{
	// note: older versions of openmp need an int as loop
	for (int i = 0; i < static_cast<int>(vs_nums_to_process.size()); ++i)
	{

Thanks for your suggestion and I will now try implementing multi-threading the weight calculation approach and will share any new findings.

[KrisThielemans]

Since the current approach has each thread handling a different viewgram independently (code below), I don't understand how "a huge amount of work is actually repeated in all those threads" .

ok. that's interesting. I think you're right. However, if we'd do list-mode reconstruction, you would no longer be sure that there is 1 view per thread (in fact, there's going to be a whole load of events per thread), so we couldn't rely on this anymore.

Even though there are repeated parts like the nested loops in the wm_calculation, i think they are unavoidable because we need to process each viewgram anyway?

Correct, but the question is in which order. An alternative would be to fill the cache upfront, but that would be impractical if you don't have enough memory.

I suppose an alternative could be to break down the compute_one_subset to compute_one_lor. Then the rest of STIR will parallise everything just fine. That would be ideal, but doesn't look trivial.

[jchan192]

Hi @KrisThielemans
I think I have made the parallelization to work in wm_calculation.
However, the problem now is that OSMAPOSL is divided into two main process. First, it executes the RPC_process_related_viewgrams_sensitivity_computation which calls wm_calculation.
Then it will execute RPC_process_related_viewgrams_gradient<true> which performs the iterative update.

I need to disable line 418 in order to test the parallelized wm_calculation in the first process. Otherwise, it will call RPC_process_related_viewgrams_sensitivity_computation in multiple threads, instead of doing multi-threading within wm_calculation. This is the code that I currently have, and it works fine, and I notice the decreased computational time.

STIR/src/recon_buildblock/distributable.cxx

Lines 415 to 419 in 19bfec5

	#ifdef STIR_OPENMP
	std::vector<double> local_log_likelihoods;
	std::vector<int> local_counts, local_count2s;
	# pragma omp parallel shared(local_log_likelihoods, local_counts, local_count2s)
	#endif

However, we only parallelized the wm_calculation. The second part of OSMAPOSL with iterative update is not parallelized because we disabled it. I tried even enabling it by modifying it to below. It simply means we allow omp parallel if we are using the gradient, and disabled it if we are not.

#ifdef STIR_OPENMP
  std::vector<double> local_log_likelihoods;
  std::vector<int> local_counts, local_count2s;
bool use_parallel = (RPC_process_related_viewgrams == rpc_gradient_true_ptr);
#  pragma omp parallel shared(local_log_likelihoods, local_counts, local_count2s) if (use_parallel)
if (use_parallel)
  std::cout << " now doing reconstruction in parallel\n";
#endif

Unfortunately, I still do not observe any speed up for the second part (iterative update) of the computation.
I ma out of idea on how to parallelized the part that performs the iterative update. Any suggestion will be greatly appreciated. Thanks.

Note that my previous method does parallelize both parts of OSMAPOSL, so im not sure what i am missing for this approach.

[jchan192]

Update. The above problem is now resolved. It was caused by an older version of OPENMP which does not support collapse(2).
The differences seem to be even smaller than my previous approach, which is encouraging.

I will do more tests tmr to confirm.

[KrisThielemans]

That sounds great! I think when we were discussing #460, we planned to get around the "double parallelisation" problem by 2 variables STIR_OPENMP_PROJECTIONS (i.e. in distributable.cxx) and STIR_OPENMP_MATRIX (used in the SPECTUB stuff). I suppose our plan was to disable the former, but enable the latter, although I'd have to look more into this.
Of course, a solution that doesn't need this would be preferable by far.

My current feeling is that you want to keep the critical sections in ProjMatrixByBinSPECTUB to ensure that there's only one wm_calculation running at once. If that is multi-threaded in itself, it doesn't matter.

Of course, the set_num_threads(1) will have to go.

It might be good to discuss this further in person.

[danieldeidda]

Hello STIR community,

I have been working on enabling OpenMP parallelization in the SPECT UB module and would like to share my changes, along with my current results. I’d appreciate any feedback on whether my implementation introduces any potential issues.

First, I removed the omp critical in this line to allow parallelization

STIR/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx

Lines 897 to 906 in 42aebde

	#ifdef STIR_OPENMP
	# pragma omp critical(PROJMATRIXBYBINUBONEVIEW)
	#endif
	if (!subset_already_processed[kOS])
	{
	if (!this->keep_all_views_in_cache)
	{
	this->clear_cache();
	subset_already_processed.assign(prj.NOS, false);
	}

Secondly , I declared two new variables , local_wm and local_wmh, at the beginning of ProjMatrixByBinSPECTUB::compute_one_subset. Then, I replace all the this->wm and this->wmh with local_wm and local_wmh in the ProjMatrixByBinSPECTUB::compute_one_subset. Since the declaration is a bit lengthy, I will show if requested.

I was testing the implementation by running OSMAPOSL, and the computational time does decrease when using 2 and 4 threads. The resulting reconstructed images do show some differences between 1 thread and multi-threads, primarily in regions with higher numerical values, but the differences are relatively small (~1e-5). Here I uploaded a slice of the difference

Do you think the approach to parallelizing SPECTUB seem reasonable? Thanks in advance, Jowett

I am a bit worried that there is any change in voxel values to be honest. These do look small but what if they are data dependent or if they increase with iterations/ Also what are the whites? can't understand from the color scale

[jchan192]

You are right, I have not tested with other set up yet, such as using other data set/ higher iteration number. I will try that after implementing the approach that parallelizes the wm_calculation

The white spots correspond to negative difference which is undefined in log scale, which is indeed misleading.Here is a new plot that shows the absolute value of the difference instead.

[jchan192]

Thanks Daniel for his availability, Here is the document that I showed today, which records the changes and performance tests that we have done prior to today's meeting.
Parallelization Development and Testing.pdf

And here is a summary for today's discussion.

• we examined the implementation of the two parallelization approaches : one that parallelizes over viewgrams, and another that parallelizes the wm_calculation.
• Depending on the number of subset, vs_nums_to_process can fluctuate, and so is the performance. More tests are needed with different num of subset to check the performance
• we should perform tests with higher subitereations. For example, 160 iterations. The more the iterations, the better the results, and the multi-threaded results should behave the same .
• The artifact appearing at the boundary of the reconstructed image is expected if attentuation mask and 3D psf are used. A suggestion is to clip off the boundary. We should perform tests with No Mask while using 3D PSF to confirm
• Things remain to be discussed: the first approach still doesnt work in list-mode reconstruction, convergence tests using cost / objective function

[danieldeidda]

Thanks! Note That I meant 160 subiterations which in your case were 40 iterations.

[jchan192]

Hi all, sorry for the delay in pushing the update to STIR. we are currently facing a minor issue that could use some help with .

When running the parallelized OSMAPOSL in STIR , the reconstruction completes and produces the reconstructed image as expected but the program does not terminate immediately. Instead, it remains active while slowly releasing cached memory. For now, we simply replace EXIT_SUCCESS with _exit(0); to end the code, which is acceptable for our internal use because the output is already generated anyway. If we plan to contribute this parallelization feature back to the STIR branch, i think its better to resolve the slow memory release issue.

After investigation, we notice that problem could be due to releasing memory of cache_collection

The slow memory release can be reproduced by running memory intensive OSMAPOSL with parallelization on, if you put the following code in the destructor function of ProjMatrixByBinSPECTUB

ProjMatrixByBinSPECTUB::~ProjMatrixByBinSPECTUB()
{
  std::cout << "start clearing cache....\n";
  clear_cache();  // <---new code here
  delete_UB_SPECT_arrays();
}

I have tried many ways to release memory of cache_collection but none seems to do the job.
Does anyone have suggestions on how to quickly release the memory of cache_collection so the program can terminate without delay?
thanks in advance