循环优化与依赖分析

// 一段普通的矩阵乘：三重循环，最内层执行 n^3 次
for (int i = 0; i < n; i++)
  for (int j = 0; j < n; j++) {
    double acc = 0.0;
    for (int k = 0; k < n; k++)
      acc += A[i][k] * B[k][j];
    C[i][j] = acc;
  }

// 循环交换：把 j 与 k 互换，让最内层访问连续内存
for (int i = 0; i < n; i++)
  for (int k = 0; k < n; k++) {
    double a = A[i][k];
    for (int j = 0; j < n; j++)
      C[i][j] += a * B[k][j];   // B 按行访问，连续
  }

// 循环倾斜：把迭代空间沿对角线切割，让依赖方向变成水平，从而可向量化
for (int i = 0; i < n; i++)
  for (int j = 0; j < n; j++)
    A[i][j] = A[i-1][j+1] + 1;      // 依赖方向 (-1, +1)，斜向

// 倾斜后：i' = i + j，依赖变成 (-1, 0)，可以按 i' 并行
for (int i2 = 0; i2 < 2 * n - 1; i2++)
  for (int j = max(0, i2 - n + 1); j <= min(n - 1, i2); j++)
    A[i2-j][j] = A[i2-j-1][j+1] + 1;

// 循环展开：因子 4，减少 3/4 的循环控制开销
for (int i = 0; i + 3 < n; i += 4) {
  s += a[i]; s += a[i+1]; s += a[i+2]; s += a[i+3];
}
for (; i < n; i++) s += a[i];      // 余数循环

// 循环分块：块大小 B 让 A 的一块与 B 的一块能同时在缓存里
for (int ii = 0; ii < n; ii += B)
  for (int kk = 0; kk < n; kk += B)
    for (int jj = 0; jj < n; jj += B)
      for (int i = ii; i < ii + B; i++)
        for (int k = kk; k < kk + B; k++) {
          double a = A[i][k];
          for (int j = jj; j < jj + B; j++)
            C[i][j] += a * B[k][j];
        }

// 同一个循环里的四类依赖
for (int i = 1; i < n; i++) {
  a[i] = a[i-1] + 1;     // 流依赖（RAW）：先读后写，方向 i-1 -> i
  b[i-1] = b[i] * 2;     // 反依赖（WAR）：先写后读
  c[i] = c[i] + 1;       // 输出依赖（WAW）：两次写
  d[i] = d[i-1];         // 输入依赖（RAR）：两次读，不影响重排
}

for (int i = 0; i < n; i++) {
  int t = a[i];      // 循环无关：只在本轮内
  b[i] = t + 1;
  c[i] = d[i-1] + 1; // 循环携带：跨迭代，阻碍并行
}

// 二维循环嵌套：A[i][j] 写，A[i-1][j+1] 读
for (int i = 1; i < n; i++)
  for (int j = 0; j < n - 1; j++)
    A[i][j] = A[i-1][j+1] + 1;

# 用方向向量判断最内层能否并行
def innermost_parallel(dirs):
    for d in dirs:
        if d[-1] != "=":        # 最内层存在携带依赖
            return False, f"最内层方向 {d[-1]}，阻碍并行"
    return True, "最内层无携带依赖，可并行或向量化"

// 判断 A[2*i] 与 A[2*j+1] 是否可能访问同一位置：解 2i = 2j + 1
// 左边恒为偶数，右边恒为奇数，无解 —— 精确判定为无依赖

# GCD 测试：若 gcd(a, b) 不能整除 (c2 - c1)，则无解
from math import gcd

def gcd_test(a, b, c1, c2):
    """访问 a*i + c1 与 b*j + c2 是否可能相交"""
    if a == 0 and b == 0:
        return c1 == c2            # 都是常量
    g = gcd(abs(a), abs(b))
    return (c2 - c1) % g == 0      # 有解才有依赖

// 源循环
for (int i = 0; i < n; i++)
  for (int j = 0; j < n; j++)
    A[i][j] = A[i][j-1] + A[i-1][j];

// 迭代空间（多面体）：
//   { (i, j) : 0 <= i < n, 0 <= j < n }
// 依赖：
//   (i, j) -> (i, j-1)  距离 (0, 1)
//   (i, j) -> (i-1, j)  距离 (1, 0)

# 检查一个调度对一条依赖是否合法
def legal(schedule, src, dst):
    """src -> dst 是依赖，要求 schedule(src) <lex schedule(dst)"""
    a, b = schedule(*src), schedule(*dst)
    return a < b        # 元组的字典序比较

# 循环交换调度：依赖 (0,1) 是否合法？
#   src=(1,0) -> sched=(0,1); dst=(1,-1)? 实际按原始坐标比较：
#   依赖 src=(i,j) -> dst=(i,j-1)，sched 后 (j,i) -> (j-1,i)，合法
print(legal(lambda i, j: (j, i), (1, 1), (1, 0)))   # True

// 可向量化：无携带依赖，访存连续
for (int i = 0; i < n; i++)
  c[i] = a[i] + b[i];

// 不可向量化：存在距离为 1 的流依赖
for (int i = 1; i < n; i++)
  a[i] = a[i-1] + b[i];

// 用 restrict 与 #pragma omp simd 帮助向量化器
void add(const double *restrict a, const double *restrict b, double *restrict c, int n) {
    #pragma omp simd
    for (int i = 0; i < n; i++) c[i] = a[i] + b[i];
}

# 让编译器报告向量化决策
clang -O3 -Rpass=loop-vectorize -Rpass-missed=loop-vectorize add.c -c
gcc -O3 -fopt-info-vec-all add.c -c

// 陷阱：这个循环看似可向量化，但 a 与 c 可能别名
for (int i = 0; i < n; i++) c[i] = a[i] + 1;
// 若 c == a + 1（重叠），向量化后语义改变
// 解决：加 restrict，或用运行期别名检查

// 运行期别名检查：两个版本，运行时分派
if (c + n <= a || a + n <= c) {
    vectorized_loop(c, a, n);      // 无重叠，走向量版本
} else {
    scalar_loop(c, a, n);          // 可能重叠，走标量版本
}