Can we sort a structure of arrays using qsort() in C?

Here is a version using Glibc's qsort_r to sort an array of indices according to a comparison function that uses a pointer to the myfluid as context, and then permutes the arrays inside the structure in-place according to the sorted indices.

If qsort_r is not available, qsort could be used, but the context would need to be passed using global variables. FreeBSD's standard C library also includes the qsort_r function.

Note: There may be portability issues with qsort_r because some versions have the parameters of qsort_r in a different order and the parameters of the comparison function in a different order. See this answer by Schwern for the sordid history of this function and its relatives.

The permutation part of the code after the qsort_r call in the tracer_sort function is stolen from this answer by AnT stands with Russia.

#define _GNU_SOURCE // for glibc's qsort_r

#include<stdio.h>
#include<stddef.h>
#include<stdlib.h>
#include<complex.h>

// 2D fluid (SoA)
typedef struct { 
  ptrdiff_t n;
  ptrdiff_t size;
  double *u0;
  double *u1;
  ptrdiff_t numTracers;
  ptrdiff_t *t_ID;
  double *t_r0;
  double *t_r1;
} myfluid;

void allocate_memory (myfluid *f){
  f->u0 = malloc (f->size * sizeof(*f->u0));
  f->u1 = malloc (f->size * sizeof(*f->u1));
  f->t_ID = malloc (f->numTracers * sizeof(*f->t_ID));
  f->t_r0 = malloc (f->numTracers * sizeof(*f->t_r0));
  f->t_r1 = malloc (f->numTracers * sizeof(*f->t_r1));
}

void initialize (myfluid *f) {
  // random velocity component in the range (-1, 1)
  srand48 (1);
  for (ptrdiff_t i = 0; i < f->size; i ++){
    f->u0[i] = -1. + 2 * drand48 ();
    f->u1[i] = -1. + 2 * drand48 ();
  }
  // random positions inside the unit square (0, 1)
  for (ptrdiff_t i = 0; i < f->numTracers; i ++){
    f->t_ID[i] = i;
    f->t_r0[i] = drand48 ();
    f->t_r1[i] = drand48 ();
  }
}

// a and b are pointers into an array of indices of type ptrdiff_t
// context is a pointer to myfluid
int compare (const void *a, const void *b, void *context){
  const myfluid *f = context;
  double r0_a = f->t_r0[*(const ptrdiff_t *)a];
  double r0_b = f->t_r0[*(const ptrdiff_t *)b];

  return (r0_a > r0_b) - (r0_a < r0_b);
}

void printout (myfluid *f) {
  printf("\n2D fluid\n");
  for (ptrdiff_t i = 0; i < f->n; i ++){
    for (ptrdiff_t j = 0; j < f->n; j ++){
      ptrdiff_t idx = j + i * f->n;
      printf ("(%.2lf, %.2lf)\t", f->u0[idx], f->u1[idx]);
    }
    printf ("\n");
  }

  printf("\nTracers\n");
  for (ptrdiff_t i = 0; i < f->numTracers; i ++){
    printf ("%td\t %.2lf\t %.2lf\n", f->t_ID[i], f->t_r0[i], f->t_r1[i]);
  }
}

void tracer_sort(myfluid *fluid,
    int (*cmp)(const void *, const void *, void *)) {
  ptrdiff_t n = fluid->numTracers;
  ptrdiff_t index[n];

  // initialize index
  for (ptrdiff_t i = 0; i < n; i++) {
    index[i] = i;
  }
  // sort index
  qsort_r(index, n, sizeof(index[0]), cmp, fluid);
  // permute the tracer elements according to the index
  // this will clobber the index
  for (ptrdiff_t i_dst_first = 0; i_dst_first < n; i_dst_first++) {
    ptrdiff_t i_src = index[i_dst_first];

    if (i_src != i_dst_first) {
      // follow the permutation cycle, clobbering the index as we go
      ptrdiff_t i_dst = i_dst_first;
      ptrdiff_t t_ID = fluid->t_ID[i_dst];
      double t_r0 = fluid->t_r0[i_dst];
      double t_r1 = fluid->t_r1[i_dst];

      do {
        fluid->t_ID[i_dst] = fluid->t_ID[i_src];
        fluid->t_r0[i_dst] = fluid->t_r0[i_src];
        fluid->t_r1[i_dst] = fluid->t_r1[i_src];
        index[i_dst] = i_dst;

        i_dst = i_src;
        i_src = index[i_src];
      } while (i_src != i_dst_first);

      fluid->t_ID[i_dst] = t_ID;
      fluid->t_r0[i_dst] = t_r0;
      fluid->t_r1[i_dst] = t_r1;
      index[i_dst] = i_dst;
    }
  }
}
 
void deallocate_memory (myfluid *fluid) {
  free (fluid->u0);
  free (fluid->u1);
  free (fluid->t_ID);
  free (fluid->t_r0);
  free (fluid->t_r1);
}

int main (int argc, char **argv){
  myfluid fluid;

  fluid.n = 8;
  fluid.numTracers = 8;

  fluid.size = fluid.n * fluid.n;

  allocate_memory (&fluid);

  initialize (&fluid);

  // sort the tracers along r0 direction
  tracer_sort(&fluid, compare);

  printout (&fluid);

  deallocate_memory(&fluid);

  return 0;
}

Note: if variable length arrays (VLAs) are not supported, replace ptrdiff_t index[n]; in tracer_sort with dynamically allocated storage, which must be freed before returning from the function.

Note about complexity: The space complexity is O(n) due to the use of an auxiliary array of indices. The sorting of the indices has time complexity O(n log n) and the subsequent permutation of the various arrays according to the sorted indices has time complexity O(n), so the overall time complexity is O(n log n).

This is not directly possible with qsort, but you could map your structure of arrays to an array of structures, which you could sort with qsort, then map those values back into a structure of arrays. It seems roundabout because it is (and begs the question of why you're not simply storing an array of structures in the first place), but it is possible.

E.g.

#include <stdlib.h>
#include <stdio.h>

struct foo { 
    int a[10];
    int b[10];
};

struct pair { 
    int a;
    int b;
};

int cmp_pair(const void *a, const void *b) {
    const struct pair *c = a;
    const struct pair *d = b;

    if (c->a < d->a) return -1;
    if (c->a > d->a) return 1;
    if (c->b < d->b) return -1;
    if (c->b > d->b) return 1;

    return 0;
}

int main(void) {
    struct foo x = {
        {1, 5, 8, 2, 4, 5, 6, 7, 8, 0},
        {4, 8, 7, 6, 1, 2, 6, 0, 9, 2}
    };

    struct pair y[10];

    for (size_t i = 0; i < 10; i++) {
        y[i] = (struct pair){ .a = x.a[i], .b = x.b[i] };
    }

    qsort(y, 10, sizeof(struct pair), cmp_pair);

    for (size_t i = 0; i < 10; i++) {
        x.a[i] = y[i].a;
        x.b[i] = y[i].b;
    }

    for (size_t i = 0; i < 10; i++) {
        printf("%d, %d\n", x.a[i], x.b[i]);
    }
}

Output:

0, 2
1, 4
2, 6
4, 1
5, 2
5, 8
6, 6
7, 0
8, 7
8, 9