# Most optimal parallel method for calculating the integral of a 2D function

I posted already this question to SO but got no answer so I try it now here:

In some crunching number program, I have a function which can be just 1 or 0 in three dimensions. I do not know in advance the function, but I need to know the total "surface" of the function which is equal to zero. In a similar problem I could draw a rectangle over the 2D representation of the map of United Kingdom. The function is equal to 0 at sea, and 1 at the earth. I need to know the total water surface. I wonder what is the best parallel algorithm or method for doing this.

I thought first about the following approach; a) divide 2D map area into a rectangular grid. For each point that belongs to the center of each cell, check whether it is earth of water. This can be done in parallel. At the end of the procedure I will have a matrix with ones and zeroes. I will get the area with some precision. Now I want to increase this precision, so b) choose the cells that are in the border regions between zeroes and ones (what is the best criterion for doing this?) and in those cells, divide them again into successive cells and repeat the process until one gets the desired accuracy. I guess that in this process, the critical parameters are the grid size for each new stage, and how to store and check the cells that belong to the border area. Finally the most optimal method, from the computational point of view, is the one that performs the minimal number of checks in order to get the value of the total surface with the desired accuracy.

• Are you looking for a heuristic? If no, what do you know about the boundary between the zero set and the one set?
– user6973
Commented Oct 23, 2011 at 21:33
• yes, some heuristic would be great Commented Oct 24, 2011 at 7:26
• what do you exactly mean about the boundary? the function to integrate or the country to calculate the area to is different on each different run of the program, if that is what you mean Commented Oct 24, 2011 at 7:27

You can partition your domain for subsequent parallel processing as follows. There are a few possible options. You can try a static geometric decomposition, including recursive bisection, quad or oct-trees and a space-filling curve.

Recursive bisection techniques are used to partition a domain (e.g., a finite element grid) into subdomains of approximately equal computational cost while attempting to minimize communication costs, that is, the number of channels crossing task boundaries A divide and conquer approach is taken. The domain is first cut in one dimension to yield two subdomains. Cuts are then made recursively in the new subdomains until we have as many subdomains as we require tasks. Notice that this recursive strategy allows the partitioning algorithm itself to be executed in parallel. Advanced variants include recursive graph bisection and recursive spectral bisection.

Here is an example. In the following picture, assume that magenta points requires twice as much work as cyan ones.

Then, you continue applying recursive bisection by cycling through the axes until # pieces = # processors, obtaining something like the following:

The idea behind quad-trees is the following one. Partition all of space into quadrants (octants in 3-space). If any piece contains too much work, partition it into quadrants, and so on. See the next picture for an example.

The best general-purpose geometric load-balancing comes from space-ﬁlling curves. The order in which points are visited in the space-ﬁlling curve determines how the geometric objects are grouped together to be assigned to the processors. The following image shows an Hilbert space-filling curve.

You problem may also be amenable to a solution using what is known as a scattered decomposition. Used when there is a structured domain space (e.g., an image) and the processing requirements are clustered, such as modeling a crash or processing an image with only a few items of interest. Suppose there are P processors. Cover the problem domain with non-overlapping copies of a grid of size P and assign each processor a cell in each of the grids. If you use more pieces:

• Less load imbalance, i.e., computation time
• More overhead and/or communication time

Finding a good tradeoff may require experimenting running the application and measuring how parallel time is affected by different scattering strategies.

Here is an example. The first picture show an example image that must be processed. Say it contains sea water and an island, and you may want to process the vegetation on the island. The image is partitioned using a very simple domain decomposition technique. The second picture shows an example of scattered decomposition.

Finally, you can use an advanced technique known as Adaptive Mesh refinement (AMR): grids are broken into blocks of ﬁxed extents, when needed blocks are reﬁned into children with same extents. Block size are chosen to optimize balance of cache efﬁciency (optimization possible since loops operate on ﬁxed size array), communication, and load balance. Often balancing blocks per processor sufﬁces; if communication is excessive then use a space-ﬁlling curve. Rebalancing requires only simple collective communication operations to decide where blocks go. The technique can be used for arbitrary dimensions. If block size is 1, adaptive blocks are quad or oct-trees. See the next picture for an example.

Hope this helps.

• thanks a lot for your really impressive and exhaustive answer Commented Oct 25, 2011 at 8:11
• you gave me a lot of new ideas I did not know the exact names for which, I will have a look at them in detail. Last thing, would you recommend any GPU implementation (if any) for AMR and quad tree decomposition? Commented Oct 25, 2011 at 8:12
• Unfortunately, I am not aware of public domain, open source GPU implementations, unless recently released. I know for sure that someone implemented AMR and Quad-trees on GPUs, since I remember reading a few papers about. Commented Oct 25, 2011 at 11:32
• ok thanks a lot; you deserve 1000 points for the answer! last thing, which tools did you use for the figures? Commented Oct 25, 2011 at 13:36
• The figures are not mine; I use them in my course on Parallel Algorithms in order to explain possible options for domain decompositions. Commented Oct 25, 2011 at 13:46