Sort2 arrays issue

Hello, dear community

I have an issue trying to sort a matrix.
I was using the subroutine sort2(n, arr, brr) being n the number of data, to sort the first array and then to rearrange the second array, but know I have a table with 3 arrays and I dont know how to modify it to sort the first array in ascending order(arr) and to rearrange the corresponding two (brr and crr) could anyone help me with this?

It is not clear where you are getting the routine SORT2() from or exactly what it does; but
it sounds like you want what is called an indexed sort; which this could be used to generate.
You could also just run SORT2() twice, I suspect. So a better description of SORT2() would help.

But as a guess, if you created a second array that is just the whole numbers to the size of the first array; like

integer,allocatable :: indx(:)
      :
 indx=[(i,i=1,size(arr))]
 call sort2(n,arr,indx)
brr=brr(indx)
crr=crr(indx)
arr=arr(indx)   ! if arr is not actually sorted by SORT2(), although it likely is
deallocate(indx)

but that is just a guess, and could be wrong. No way to tell without knowing more details about SORT2(); and I think it is likely if you have one sort routine you have others which are likely to do this more efficiently. I am assuming SORT2() sorts array ARR and does the same element moves on BRR that it does on ARR when called with “SORT2(ARR,BRR)” ? No way of telling if BRR can be INTEGER or not from the description so far, but that can be gotten around.

There are several packages available from various sources; several of which are on github or other places and freely available for doing sorts like what I think you want. If you use fpm(1) or stdlib(1) there are several appropriate sort routines at your fingertips.

Edit: Added N to SORT per below; but see below.

1 Like

This sounds like the sort2() from “Numerial Recipes” :wink:

The solution proposed by @urbanjost is mostly OK, except that it misses the size in the sort2() call, and that the second array has to be real. This can be fixed, but it’s actually simpler to use the indexx() routine from the same source, which purpose is precisely to build the index array.

integer, allocatable :: indx(:)
allocate( indx(n) )
call indexx(n,arr,indx)
arr = arr(indx)
brr = brr(indx)
crr = brr(indx)

Note that each array indexing construct will generate a temporary array allocation, which could make them inefficient if you have very large arrays.

The alternative (without temporary allocations) would be to modify sort2() and create sort3(n,arr,brr,crr) with an additional array to be sorted. This is actually quite simple, you just jave to spot where brr is used and do the same things with crr.

But at the end, if your arrays are small ones, repeatedly calling sort2() is also perfectly reasonable.

Sorting the elements, or computing an indx(:) array that ranks the elements, is an O(n*log(n)) kind of process. Once the indx(:) array has been computed, then rearranging the elements is just an O(n) process. So modifying a sort routine to do multiple rearrangments is almost always inefficient compared to the two-step process. There are some exceptions, such as when it is known that only one or two elements are out of order, in which case the sort step requires <n effort to begin with. Of course, the final rearrangement step using the indx(:) array usually requires some temp workspace memory allocation, it is tricky to do the general rearrangement in-place, but if the programmer has available a work array anyway, then that is not a problem.

2 Likes

I can recommend a sorting package I’ve put together, sortff (I’m happy to accept contributions as well). There’s a video that explains the ideas behind it as well, Vector Subscripts For Fortran Array Access - YouTube

The code will then look similar to some suggestions above. I.e

use sortff, only: sorted_order

associate(order => sorted_order(arr))
  arr = arr(order)
  brr = brr(order)
  crr = crr(order)
end associate

You could easily put such an operation in a subroutine (e.g. sort3) if desired.

1 Like

You may consider some improvements for small-size arrays. I’ve done some tests using sorting networks (you can try those here), but it turns out that, unless the sorted array size is very small (say <=4), a loop is always better than a sorting network.

gcc/gfortran also very good at optimizing, such that sorting network are never really necessary, see for example here

My understanding is that it indeed can provide some performance benefits to manually handle small array cases, but unless you are sorting lots of small arrays, you probably won’t really see an impact from the lower performance of the general case. Also, I believe it’s likely that it would be somewhat system and hardware specific exactly how to get peak performance. So if it turns out you really do need a highly tuned implementation, there’s no library that’s going be exactly what you need, because how to tune it will almost always be specific to your particular case. But if you’re just trying to get something working, pick up a library that makes that part easy.

From what I see, the sorting network only improves over the worst-case scenario (fully reverted array), because, even though the compiler can optimize the sequence of the sorting net for SIMD, it cannot skip any operations!

While the loops do, if numbers are already partially sorted.