[コンパイラ] 17. 並列性検出とループ変換

// データ並列：各反復が独立
for (i = 0; i < n; i++) {
    c[i] = a[i] + b[i];
}

// 4つのプロセッサで分割実行：
// P0: i = 0..n/4-1
// P1: i = n/4..n/2-1
// P2: i = n/2..3n/4-1
// P3: i = 3n/4..n-1

// タスク並列：3つの関数が独立
result1 = processA(data)    // タスク1
result2 = processB(data)    // タスク2
result3 = processC(data)    // タスク3
// 3つのタスクを同時に実行可能

// 3段パイプライン：
// Stage 1: read    -> P0
// Stage 2: process -> P1
// Stage 3: write   -> P2

// 時間：  1      2      3      4      5
// P0:   read1  read2  read3  read4  read5
// P1:          proc1  proc2  proc3  proc4
// P2:                 wrt1   wrt2   wrt3

// フロー依存性（真の依存性）
for (i = 1; i < n; i++) {
    a[i] = a[i-1] + 1;  // 反復iが反復i-1の結果に依存
}
// 依存距離（distance） = 1、並列化不可

// 逆依存性
for (i = 0; i < n-1; i++) {
    a[i] = a[i+1] + 1;  // 読み取り後書き込み
}
// 順序を逆転すれば解決可能

// 出力依存性
for (i = 0; i < n; i++) {
    a[i % 2] = b[i];    // 同じ位置に繰り返し書き込み
}

// 2重ループでの依存ベクトル
for (i = 1; i < n; i++) {
    for (j = 1; j < m; j++) {
        a[i][j] = a[i-1][j-1] + 1;
    }
}

// a[i][j]はa[i-1][j-1]に依存
// 依存ベクトル：(1, 1)（i方向距離1、j方向距離1）

// a[c1*i + c2]とa[c3*j + c4]が同じ位置を参照するには
// c1*i + c2 = c3*j + c4を満たす整数i、jが存在しなければならない
// -> gcd(c1, c3)が(c4 - c2)を割り切る必要がある

// 例：
// a[2*i + 1]とa[2*j]
// gcd(2, 2) = 2、c4 - c2 = 0 - 1 = -1
// 2は-1を割り切れない -> 依存性なし！（並列化可能）

// ループ範囲を考慮して依存性を検査
// a[i + 1]とa[j]（0 <= i,j < n）
// i + 1 = jを満たす(i,j)ペアが範囲内に存在するか？
// 0 <= i < nかつ0 <= i+1 < nなら存在
// -> 依存性あり

// 正確な依存性分析が難しい場合：

// 1. 間接配列アクセス
a[b[i]] = ...  // b[i]の値が不明

// 2. ポインタを通じたアクセス
*p = ...       // pが指す位置が不明

// 3. 関数呼び出し
foo(a, i)      // fooがaを修正するか不明

// こうした場合保守的に依存性ありと仮定（安全だが並列化不可）

// 元：行優先アクセス（Cでは既に最適）
for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
        a[i][j] = 0;

// 交換後：列優先アクセス
for (j = 0; j < m; j++)
    for (i = 0; i < n; i++)
        a[i][j] = 0;

// 依存ベクトル(1, -1)の場合：
// 交換後(-1, 1) -> 第1成分が負 -> 不正！
//
// 依存ベクトル(1, 0)の場合：
// 交換後(0, 1) -> 辞書式正 -> 合法！

// 元：行列乗算
for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
        for (k = 0; k < n; k++)
            c[i][j] += a[i][k] * b[k][j];

// タイリング後（タイルサイズT）：
for (ii = 0; ii < n; ii += T)
    for (jj = 0; jj < n; jj += T)
        for (kk = 0; kk < n; kk += T)
            for (i = ii; i < min(ii+T, n); i++)
                for (j = jj; j < min(jj+T, n); j++)
                    for (k = kk; k < min(kk+T, n); k++)
                        c[i][j] += a[i][k] * b[k][j];

// 元：2つの別々のループ
for (i = 0; i < n; i++)
    a[i] = b[i] + 1;

for (i = 0; i < n; i++)
    c[i] = a[i] * 2;

// 融合後：一つのループ
for (i = 0; i < n; i++) {
    a[i] = b[i] + 1;
    c[i] = a[i] * 2;
}

// 元
for (i = 0; i < n; i++) {
    a[i] = b[i] + 1;        // 独立部分1
    c[i] = c[i-1] + a[i];   // 依存部分2
}

// 分裂後
for (i = 0; i < n; i++)
    a[i] = b[i] + 1;        // 並列化可能！

for (i = 0; i < n; i++)
    c[i] = c[i-1] + a[i];   // 逐次実行必要

// 元の反復空間：(i, j)
for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
        S: a[i][j] = ...

// アフィン変換：新しい反復変数 (p, q) = T * (i, j)
// 例：(p, q) = (i+j, j) -> 傾いた反復空間

// 元：依存ベクトル(0,1)と(1,0) -> 両次元とも逐次的
for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
        a[i][j] = (a[i-1][j] + a[i][j-1]) / 2;

// スキュー：p = i + j, q = j
// 変換後同じp値を持つ反復は並列実行可能

// 変換後の依存ベクトル：(1,1), (1,0) -> p（外部ループ）は逐次、q（内部ループ）は並列可能

// 多面体表現：
// 文Sの反復ドメイン：0 <= i < n, 0 <= j < m
// これは2次元多面体（長方形）

// スケジュール関数：theta(i,j) = (i, j)    （元の順序）
//                 または theta(i,j) = (i+j, j) （スキュー）
//                 または theta(i,j) = (j, i)   （交換）

// 時間的局所性が良いコード
for (i = 0; i < n; i++) {
    sum += a[i];  // sumが毎反復再利用される
}

// 時間的局所性が悪いコード（タイリングで改善可能）
for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
        b[j] += a[i][j];  // b[j]が外部ループでのみ再利用

// 空間的局所性が良いコード（行優先アクセス、C配列レイアウトと一致）
for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
        a[i][j] = 0;      // a[i][0], a[i][1], ... の順でアクセス

// 空間的局所性が悪いコード（列優先アクセス）
for (j = 0; j < n; j++)
    for (i = 0; i < n; i++)
        a[i][j] = 0;      // a[0][j], a[1][j], ... の順でアクセス
                           // 毎アクセスごとにキャッシュラインが異なる

// 最適化なし：O(n^3) キャッシュミス
// ループ交換：空間的局所性改善
// タイリング（B x B）：O(n^3 / B) キャッシュミス
// タイリング + 交換：最適に近い性能

// 実際の性能差（n = 1000、L1キャッシュ32KB基準）：
// 最適化なし：    ~10秒
// ループ交換：    ~5秒
// 32x32タイリング：~1秒
// 性能差は10倍に達し得る！

// DOALL：完全に並列化可能
for (i = 0; i < n; i++)
    a[i] = b[i] + c[i];

// OpenMPへの自動変換：
// #pragma omp parallel for
// for (i = 0; i < n; i++)
//     a[i] = b[i] + c[i];

// DOACROSS：依存距離分の並列性が存在
for (i = 1; i < n; i++)
    a[i] = a[i-1] + b[i];  // 依存距離1

// パイプライン並列化：
// P0: a[0], a[4], a[8], ...
// P1:（待機）a[1], a[5], a[9], ...
// P2:（待機）（待機）a[2], a[6], ...
// P3:（待機）（待機）（待機）a[3], a[7], ...