[MITgcm-devel] vectorization of layers

[Martin Losch]

Hi Ryan,

I didn't want to make it sound as if the code is inefficient. The code is probably efficient for the type of current cpu that most people use, but for vector cpus it is absolutely important that the code "vectorizes"; the "scalar" performance of these computers is very poor. I am attaching an example to illustrate how very small changes can speep up the code dramatically (especially in the calc_r_star example).

For vectorization to work, the operations within a loop need to be independent of each other. Often I don't understand why something vectorizes or not, but there are analysis tools that help with that. In the present case the while loops are not vectorizable, because they clearly depend on the previous iterate (the while conditions is a function kgu, that is begin computed within the loop). Usually, what I do in cases like these is to move the problematic loop outside the i,j loops. In this case it is a little more difficult, because the while iteration will be different at each (i,j) point. 
This is how I think this need to be changed in order to vectorize:
- turn TatU/V into 2D arrays and precompute them for each kk (small change)
- split the if-block starting at l143 ("IF (TatU .GE. layers_bounds(Nlayers,iLa)) THEN"). The first 2 (3, but in the third case nothing happens) case are simple and can remain in the i,j-loop. The 4th and 5th case (the ones with the while-loops) could be replace by a second kk-loop and appropriate if statements, like this;
DO kk
DO i,j
IF ( appropriate conditions) THEN kgu(i,j) = something
otherwise do nothing
ENDDO
ENDDO
But this is where I am not sure if it works at all. Maybe it would be good ot set up or use a test-experiment (in the spirit of the verification experiments) that checks the results of the layers package, so that I know when I break the algorithm. 
Do you have something simple like that? That would already help greatly. I would then try to come up with a solution and you can tell me if you like it. Most likely, this solution will involve a "ifdef VECTORIZE then do the new code, #else do the old code".

Martin

On Jan 4, 2013, at 6:51 PM, Ryan Abernathey <ryan.abernathey at gmail.com> wrote:

> I wrote this code. What is going on here is the most important and expensive part of the LAYERS package: the sorting of each point in the water column into the appropriate density class. This is essentially a search problem (it has to hunt for the correct bin), and I based the code on one of the search algorithms from Numerical Recipes. I did some basic testing and optimization of the code in MATLAB before putting it into the model, but I am sure there is lots of room for improvement.
> 
> The "hunting" takes place inside the while loops. kgu is an index for layers_bounds that tracks which target density bin is being tested. If the bin is too cold for the current point (whose density value is TatU), kgu is augmented. The reason for using a while-loop is that it only runs as long as it has to until the search condition is met. This seemed like the fastest thing at the time.
> 
> This search has to run on each point in the refined vertical grid (Nr * FineGridFact), so any improvement to its execution time would really multiply.  
> 
> I have only a vague understanding of what vectorization is, but I would be happy to work with someone who understands such things better in order to make the code more efficient. 
> 
> -Ryan
> 
> 
> On Fri, Jan 4, 2013 at 10:58 AM, Martin Losch <Martin.Losch at awi.de> wrote:
> Hi there,
> I am not sure how much effort one still should invest into vectorizing the code, since our vector machine will only run for another year or so, but can you think of a good way of rewrite layers_fluxcalc.F so that it vectorizes? Currently these lines (and the corresponding for TatV):
> 
>          ELSE IF (TatU .GE. layers_bounds(kgu(i,j),iLa)) THEN
> C        have to hunt for the right bin by getting hotter
>           DO WHILE (TatU .GE. layers_bounds(kgu(i,j)+1,iLa))
>            kgu(i,j) = kgu(i,j) + 1
>           ENDDO
> C         now layers_bounds(kgu(i,j)+1,iLa) < TatU <= layers_bounds(kgu(i,j)+1,iLa)
>          ELSE IF (TatU .LT. layers_bounds(kgu(i,j)+1,iLa)) THEN
> C        have to hunt for the right bin by getting colder
>           DO WHILE (TatU .LT. layers_bounds(kgu(i,j),iLa))
>            kgu(i,j) = kgu(i,j) - 1
>           ENDDO
> C         now layers_bounds(kgu(i,j)+1,iLa) <= TatU < layers_bounds(kgu(i,j)+1,iLa)
>          ELSE
> 
> inhibit vectorization (because the while-loops don't vectorize), and thus making it impossible to use the layers-pkg on a vector computer (This routines has a vector operation ratio of 5% and takes 86% of the computer time). Currently I have no clue what to do here, but I don't quite understand, what's going on in thie package/routine.
> 
> Martin
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
> 
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel