能否在C语言中用qsort()对SoA形式的被动示踪剂排序？

我正在开发一款湍流模拟代码，需要利用流体速度分量平流输运一组被动示踪剂。目前流体采用**数组结构（SoA）定义，方便后续调用FFTW做速度分量变换；被动示踪剂则用结构数组（AoS）**定义，这样能按任意方向快速排序。

我写了如下代码：

#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;
} myfluid;

// Passive tracers (AoS)
typedef struct { 
  ptrdiff_t ID;
  double r0;
  double r1;
} mytracer;

void allocate_memory (myfluid *f, mytracer **t){
  f->u0 = malloc (f->size * sizeof(double));
  f->u1 = malloc (f->size * sizeof(double));

  *t  = malloc (f->numTracers * sizeof(mytracer));
}

void initialize (myfluid *f, mytracer *t) {
  // 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 ++){
    t[i].ID = i;
    t[i].r0 = drand48 ();
    t[i].r1 = drand48 ();
  }
}

int compare (const void *a, const void *b){
  mytracer *A = (mytracer *)a;
  mytracer *B = (mytracer *)b;
  return (A->r0 > B->r0) - (A->r0 < B->r0);
}

void printout (myfluid *f, mytracer *t) {
  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", t[i].ID, t[i].r0, t[i].r1);
  }
}
 
void deallocate_memory (myfluid *fluid, mytracer **t) {
  free (fluid->u0);
  free (fluid->u1);

  free (*t);
}

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

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

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

  allocate_memory (&fluid, &tracer);

  initialize (&fluid, tracer);
  
  // sort the tracers along r0 direction
  qsort(tracer, fluid.numTracers, sizeof(mytracer), compare);

  printout (&fluid, tracer);

  deallocate_memory(&fluid, &tracer);

  return 0;
}

请问是否可以将示踪剂也改为SoA形式，同时不丢失当前的快速排序能力？

回答： #

完全可以将示踪剂改为SoA结构，同时保留按任意维度排序的能力，核心思路是维护一个索引数组，通过排序索引而非直接排序数据数组来实现需求，具体方案如下：

1. 定义SoA形式的示踪剂结构 #

把原来的AoS结构体拆分成独立的数组，保证内存连续，适配后续的向量化或FFTW操作：

// 被动示踪剂（SoA）
typedef struct {
  ptrdiff_t numTracers;
  ptrdiff_t *ID;
  double *r0;
  double *r1;
} mytracer_soa;

2. 用索引数组实现排序逻辑 #

不需要直接对SoA的r0/r1数组进行排序（这会破坏不同数组间的数据关联），而是创建一个ptrdiff_t类型的索引数组indices，其中indices[i]表示排序后第i个位置对应的原始示踪剂下标。排序时仅操作这个索引数组，通过比较索引指向的r0/r1值来实现排序。

比如针对r0方向升序排序的比较函数（使用POSIX标准的qsort_r来传递SoA上下文，保证线程安全）：

// 按r0升序排序的比较函数（针对索引数组）
int compare_r0(const void *a, const void *b, void *arg) {
    mytracer_soa *t = (mytracer_soa*)arg;
    ptrdiff_t idx_a = *(const ptrdiff_t*)a;
    ptrdiff_t idx_b = *(const ptrdiff_t*)b;
    return (t->r0[idx_a] > t->r0[idx_b]) - (t->r0[idx_a] < t->r0[idx_b]);
}

如果编译器不支持qsort_r，也可以将SoA结构设为全局变量配合标准qsort使用，但前者是更优的方案。

3. 完整改造后的核心代码示例 #

内存分配函数

void allocate_memory_soa(myfluid *f, mytracer_soa *t) {
    f->u0 = malloc(f->size * sizeof(double));
    f->u1 = malloc(f->size * sizeof(double));

    t->ID = malloc(t->numTracers * sizeof(ptrdiff_t));
    t->r0 = malloc(t->numTracers * sizeof(double));
    t->r1 = malloc(t->numTracers * sizeof(double));
}

初始化函数

void initialize_soa(myfluid *f, mytracer_soa *t) {
    srand48(1);
    // 初始化流体速度
    for (ptrdiff_t i = 0; i < f->size; i++) {
        f->u0[i] = -1. + 2 * drand48();
        f->u1[i] = -1. + 2 * drand48();
    }
    // 初始化示踪剂SoA数据
    for (ptrdiff_t i = 0; i < t->numTracers; i++) {
        t->ID[i] = i;
        t->r0[i] = drand48();
        t->r1[i] = drand48();
    }
}

主函数中的排序与访问

int main(int argc, char **argv) {
    myfluid fluid;
    mytracer_soa tracer;
    ptrdiff_t *indices;

    fluid.n = 8;
    tracer.numTracers = 8;
    fluid.size = fluid.n * fluid.n;

    allocate_memory_soa(&fluid, &tracer);
    // 分配并初始化索引数组
    indices = malloc(tracer.numTracers * sizeof(ptrdiff_t));
    for (ptrdiff_t i = 0; i < tracer.numTracers; i++) {
        indices[i] = i;
    }

    initialize_soa(&fluid, &tracer);
    
    // 按r0方向排序索引数组
    qsort_r(indices, tracer.numTracers, sizeof(ptrdiff_t), compare_r0, &tracer);

    // 按排序后的顺序访问示踪剂
    printf("\n排序后的示踪剂（按r0升序）\n");
    for (ptrdiff_t i = 0; i < tracer.numTracers; i++) {
        ptrdiff_t idx = indices[i];
        printf("%td\t %.2lf\t %.2lf\n", tracer.ID[idx], tracer.r0[idx], tracer.r1[idx]);
    }

    // 释放内存
    free(indices);
    free(tracer.ID);
    free(tracer.r0);
    free(tracer.r1);
    free(fluid.u0);
    free(fluid.u1);

    return 0;
}

4. 方案优势 #

• 保留SoA优势：内存连续布局，适合后续向量化操作（如SIMD），也与流体的SoA结构统一，便于代码维护。
• 不丢失排序能力：只需修改比较函数，就能实现按r0、r1甚至ID的任意维度排序。
• 性能更优：排序时仅操作小尺寸的索引数组（ptrdiff_t通常为8字节），比排序AoS结构体（每个至少24字节）更快，缓存命中率更高。

内容的提问来源于stack exchange，提问作者ashwin