Solver routines are not OpenMP-safe

[conradtchan]

Now that the issue with the OpenMP compiler flags has been resolved (#220), we're ready to test the performance benefits of OpenMP at the Fortran level.

Unfortunately, it's not as straightforward as applying some omp do loops, e.g.

!$omp parallel do private(k, ierr)
do k = 1, nsol
   call newton_tanh(c(k,:), abu(k,:,:), nstar, ierr)
   if (ierr == 1) then
      call psolve_tanh(c(k,:), abu(k,:,:), nstar)
      call newton_tanh(c(k,:), abu(k,:,:), nstar, ierr)
   endif
enddo
!$omp end parallel do

because newton_tanh calls set_abu_data which dynamically allocates a module variable - which isn't OpenMP-safe. The solver routines need to be re-written before OpenMP parallel do loops can be used.

[2sn]

@conradtchan Thanks. What is the usual way around that in OpenMP codes? We stuff data in an object that we allocate and pass a pointer to it - or do all as object methods?

[conradtchan]

The most straightforward way is to allocate the abu variable outside of the omp loop, declare it as private, and then pass it through as an argument.

Independent of the issues with OpenMP, the current implementation seems inefficient because abu is repeatedly allocated and deallocated as it loops over the solutions, even though nstar and nel don't change.

But perhaps I don't fully understand the design. What was the reason for making abu a module variable?

[2sn]

I think at the time I did not want to pass complex data.

[2sn]

It may have been that is it is used as auxiliary data by the Powell algorithm.

[conradtchan]

Would you like me to refactor this routine to pass the abundance as an argument? Or at least make the changes needed to demonstrate whether we get a performance benefit from using an omp parallel loop.

[2sn]

Maybe one can create some kind to type / class to contain all the data. The module data was convenient to keep the code structure simple at first, but not good for, e.g., functional programming (as it should ideally be) with pure functions and parallelisation.

Maybe keep changes in a branch for demo?

In the long run, I think there will be need for the GA to be in Fortran.

StarFit is now being increasingly picked up by the community ...

[conradtchan]

Ah, I see why it's a module variable now. The root finding subroutine doesn't take auxiliary variables, so the abundances had to be passed in that way.

[2sn]

I assume it can be updated, I already did a few modifications to it to make the code more modern.

As mentioned earlier, probably best to create some object/data type in Fortran to all all the data for a star, and maybe one for an abundance set and a container object with a bunch of them, and then maybe pass these. Could be derived from some generic template type that the finder knows, and passes on to the optimisation routine. Would be some major refactor, though. The fitness routines would then take that as input rather than taking it from a common block. I think, though, there is a lot more data that is computed for each set when loaded, some inverse matrices, etc. They all would have to be part of the data as well.

On the other hand. I think the data that is computed for the star and stored in the module data is actually not written by any of the optimisations, so it could be fixed read-only. Could that be done? But model data, I agree, would need to be passed in case t was stored in a common block, that design would have to change.