Confused about transposes / strides

[peterukk]

Hi,

After reading "When to transpose data" I'm afraid I'm still a bit confused. It seems the "normal" way to use FTorch is to have your Fortran and PyTorch arrays be similar, and have FTorch get around the column/row major differences of the languages by working out the appropriate strides in torch_tensor_from_array. But wouldn't this imply a performance penalty or a transpose somewhere else? For example, SIMD vectorization is generally much faster with unit strides. I imagine the PyTorch model wants to work with contiguous arrays for various optimizations e.g. vectorization. If the Fortran and PyTorch code are set up with similar arrays, there must be a transpose or non-unit-stride access happening somewhere.

I've previously used the PyTorch-Fortran library and with that I constructed my Fortran input arrays with the transposed layout w.r.t. PyTorch, for instance my_fortran_array(nx, nbatch) where nx is the contiguous dimension in memory and my understanding is that the library takes a pointer and will simply re-interpret the array as pytorch_array(nbatch, nx) where again nx is the contiguous dimension in physical memory, as it should be. How can I replicate this behaviour with FTorch and avoid any under-the-hood transposes / non contiguous memory access?. Do I do what I did before but pass layout=[2,1] to torch_tensor_from_array? A bit brain-dead today, sorry if I missed something obvious!

[jatkinson1000]

Hey @peterukk thanks for asking.

So, you are correct that the default behaviour of FTorch is to have users index their arrays identically between Fortran and PyTorch with the aim of not having to worry about whether they have transposed correctly or not, or which point in the process to transpose.

To achieve this we make use of strided access in Torch - we want to avoid doing data transfer or augmentation as these can also be costly processes so have Fortran and Torch read from the same data in memory using pointers as you state. Under the hood FTorch will calculate the appropriate strides to use based on the dimensions of your array. See the code here in the versatile torch_tensor_from_blob() routine.

Now, as you say, this could itself incur overheads, particularly in the form of cache misses compared to sequential reading in memory. Whether this will have an impact, and whether this is larger than the impact of transposition operations very much depends on the size of the arrays you are using, and the specification of the machine (cache sizes etc.). There is no definitive ruling, and whilst we have not found it to be a bottleneck in projects so far it is something we have talked about exploring with some benchmarking comparing transpose vs. stride for different size arrays on different machines.
We chose to stick with 'matching indexing' as our default approach since one of the aims of FTorch is to be easily usable by a wide non-technical audience. However...

If this is a concern we do provide users who are comfortable with the ability to index the arrays as they like. This is done with the layout argument you used previously (made optional in #348). To use this take the transpose in Fortran and then call torch_tensor_from_array() with the following argument list:

torch_tensor_from_array(tensor, data_in, layout, device_type, device_index, requires_grad)

(with device_index and requires_grad being optional).

From the sound of your previous code this is perhaps not quite what you are after as if you have already performed a transpose in Fortran and call torch_tensor_from_array() with layout=[2, 1] (or more generally layout=[n, n-1, ... 2, 1]) to get the expected values but with data now being read sequentially in memory from the Torch side you will still have an array with the same shape as the Fortran transpose in Torch.
Of course, we always encourage testing to verify numbers appear as expected!

For more control you could look at torch_tensor_from_blob which gives users full control over shape and strides.

Hopefully this helps, but please do ask any follow-up questions, and if you have any suggestions on how we can clarify this - e.g. changes/additions to the documentation or the introduction of an example showing this we would be happy to incorporate them. The source for the transposing docs can be found here and contributions are welcome!

[peterukk]

Thanks Jack that was very helpful! I hope to benchmark the performance difference between the two approaches. It could be that the non-sequential reading of the input arrays, and non-sequential writing of output at the end, is very rarely a bottleneck in practice as modern ML libraries such as PyTorch are liberal with making new memory allocations, assuming cost is dominated by matrix-matrix multiplies (a good assumption for anything but truly small neural networks). So I guess there would be only one non-sequential read at the beginning of the model, but intermediate arrays allocated by PyTorch would be contiguous in memory. Still, this performance aspect could be worth testing for someone like myself working with fairly small RNNs, which may or may not directly operate on the input array.

[peterukk]

I'm trying to understand the stride calculation in torch_tensor_from_blob. Let's say nx=16, nbatch=100 using the previous example where nx should be memory-contiguous dimension, so tensor_shape = [16,100] for the Fortran array.
layout = [1,2] --> strides = 1 16
layout = [2,1] --> strides = 100 1
Why does strides=(100,1) lead to data being read sequentially in memory in Torch for an array whose row-major python shape is (100,16)?

[jatkinson1000]

Hey @peterukk

I wrote a bunch of stuff below, but to be honest whilst self-checking what I said I came across this article that I don't feel I can improve upon to be honest: https://martinlwx.github.io/en/how-to-reprensent-a-tensor-or-ndarray/

---

In the docs for the C++ implementation of from_blob strides specifies the strides to be taken to move to the next item in each dimension.

For row-major languages dimension 0 is always the "slowest" dimension, and dimension n is always the fastest.

As is often the case, the PyTorch docs are usually better explained with examples, however:

>>> x = torch.tensor([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10]])
>>> x.stride()
(5, 1)

The strides indicates the number of items to move in memory to get to the next item in the corresponding dimension.
Since PyTorch is row major then to move to the next item in a column (dimension 0) would require a step of n_cols (100 in your case) whilst moving to the next item in a row (dimension 1) requires a step of 1.

Click to see more elaborate demo of dimension and stride order expected in Torch.

>>> import torch
>>>
>>> my_tensor = torch.triu(torch.ones(5, 3))
>>> my_tensor
tensor([[1., 1., 1.],
        [0., 1., 1.],
        [0., 0., 1.],
        [0., 0., 0.],
        [0., 0., 0.]])
>>>
>>> my_tensor.dim()
2
>>>
>>> my_tensor.size()
torch.Size([5, 3])
>>>
>>> my_tensor.dim_order()
(0, 1)
>>>
>>> my_tensor.stride()
(3, 1)

So for a torch tensor of shape (Size) [5, 3] where dimension 0 is along columns ($n$ = row number) and dimension 1 is along rows ($m$ = column number) we expect strides of (3, 1).

For contiguous row-major data then moving 1 'item' in memory will take us to the next element in a row, i.e. stepping in dimension 1, so the strides to move in dimension 1 is 1, whilst the stride to move in dimension 0 (down a column) is 3.

In contrast, if the data was saved from column-major Fortran using the same layout then the stride to move along dimension 0 would be 1, whilst the stride to move along dimension 1 (rows) would be equal to n_rows (16 in your example). This is what is generated by the default layout = [1, 2] for an 'FTorch tensor'.

---

Note that Torch sometimes chooses to leave tensors alone in memory rather than reshaping for some operations, and uses strided access to represent this. e.g. taking the transpose:

Click to see strided access in PyTorch for transpose.

>>> import torch
>>>
>>> my_tensor = torch.triu(torch.ones(5, 3))
>>> my_tensor
tensor([[1., 1., 1.],
        [0., 1., 1.],
        [0., 0., 1.],
        [0., 0., 0.],
        [0., 0., 0.]])
>>>
>>> my_tensor.dim()
2
>>>
>>> my_tensor.size()
torch.Size([5, 3])
>>>
>>> my_tensor.dim_order()
(0, 1)
>>>
>>> my_tensor.stride()
(3, 1)
>>>
>>> my_transposed_tensor = my_tensor.transpose(0, 1)
>>>
>>> my_transposed_tensor
tensor([[1., 0., 0., 0., 0.],
        [1., 1., 0., 0., 0.],
        [1., 1., 1., 0., 0.]])
>>>
>>> my_transposed_tensor.dim()
2
>>>
>>> my_transposed_tensor.size()
torch.Size([3, 5])
>>>
>>> my_transposed_tensor.dim_order()
(1, 0)
>>>
>>> my_transposed_tensor.stride()
(1, 3)
>>>

So it may well be that with knowledge of the stride pattern the internals of Torch are set up to do sensible things, though I have no basis for that claim.

---

This article from PyTorch may also clarify some ideas (as well as showing things can get more complicated internally for certain models/operations like CNNs).

[peterukk]

@jatkinson1000 Thanks so much. Lol to be honest I'm still a bit confused - I guess I thought strides should be (16,1) because

>>>  torch.zeros((100,16)).stride()
(16, 1)

The (100,16) sized torch array would be equivalent to (16,100) in Fortran, which was our array in this example, so what's the difference? (a short answer will suffice, you've spent enough time on this!)

Edit: I can see torch_from_blob_c is given as third argument the Fortran shape, not C++ / Python shape - does two wrongs somehow make a right?

[jatkinson1000]

Hi @peterukk I think I now understand what you are after.

You can do this in FTorch, but it requires you to use torch_tensor_from_blob() (the underlying routine for tensor transformation that gives 'power users' full control as described in the docs) and declare the shape and stride pattern yourself.

I can see that it might be useful to have a more generic interface to achieve this, however.

As such, I have implemented it on the C-layout branch.

Please see the description of how to use the default layout, the reverse layout, and the new permute_shape option with reverse layout (which I believe is what you want here) in my comment on the pull request.

From your install of FTorch please run:

git fetch
git checkout C-layout

and then rebuild FTorch as normal.

After this please see if you can get the results you'd expect by:

• Transposing your array in Fortran
• Generating a tensor using
  torch_tensor_from_array(tensor, transposed_array, layout=[2, 1], permute_shape=.true., ...)

If so that's great! (If not, let us know so we can get it working as needed).
I'll then work to complete the PR ready to merge into main by adding the necessary documentation etc.

Any feedback on this process (including my question on the PR about whether it would be preferable for the transpose to be done by FTorch instead of requiring you to do it first) would be really appreciated.

[peterukk]

@jatkinson1000 I have tested the new feature and it works! Thanks a lot for the amazing support!

Speed difference to my old inference code based on PyTorch-Fortran is less than 10%, pretty sure it's working as intended. As expected, doing layout = [2,1] without permute_shape=.true. led to a crash.

I hope to still time this against "default FTorch" with layout=[1,2], permute_shape=.false. but I'm pretty sure that would be slower..what the new feature allows for is contiguous memory access within PyTorch and minimal, user-controlled "transposes". I put transposes in quotation marks because there need not be any transposes with layout=[2,1], permute_shape=.true. It depends on whether the host Fortran model uses a similar array layout as PyTorch (and ML in general) with columns or "batches" being outermost in physical memory. In E3SM, and probably most GCMs out there, that's not the case - columns are innermost, yet as I said before, even then the PyTorch input arrays don't necessarily correspond to any GCM arrays one-to-one so there is still no need for explicit transposes, just implicit ones (non-contiguous read/write operations when "packing the inputs" and "unpacking the outputs". Example with my code:

do i=1,ncol
  do k=1,nlev
    input_torch_rnn_lev(1,k,i) =  state%t(i,k) ! Temperature
    input_torch_rnn_lev(2,k,i) =  state%t(i,k) ! Water vapor
    ...
  
call torch_tensor_from_array(in_tensors(1), input_torch_rnn_lev, layout=[3, 2, 1], device_type=torch_kCPU, permute_shape=.true.)
call torch_tensor_from_array(out_tensors(1), output_torch_rnn_lev, layout=[3, 2, 1], device_type=torch_kCPU, permute_shape=.true.)

call torch_model_forward(model, in_tensors, out_tensors)

do k = 1, nlev
  ptend%s(1:ncol,k) = output_torch_rnn_lev(1,k,1:ncol) ! Temperature tendency
  ...  

There is no reason to assume the Fortran model will always have columns innermost in memory, if they are outermost then we could avoid these non-contiguous memory accesses too.

To answer your question in the pull request, "should we require users to transpose the data themselves before calling this routine, as is the current implementation, or should we instead make the Fortran transpose part of the internals?" - I don't fully understand this because the premise is a Fortran transpose is always needed, which is not the case? But no, I would strongly favor users having control of any transposes - in the best case there are not even implicit transposes (non-contiguous reads/writes when preparing and unpacking the inputs and outputs) needed. The performance-oriented user that wants to do this needs to of course be aware of Fortran being column major and PyTorch being row major, so that the input arrays on the Fortran side look different (even if they're actually the same) - that is a very reasonable expectation especially if this information is in the docs