第8章：脉动阵列验证方法

测试矩阵维度分类：

1. 最小矩阵：1×1×1，验证退化情况
   - 测试单个PE的功能
   - 验证控制逻辑的最简路径
   - 确保特殊情况不会导致死锁

2. 小于阵列：M,K,N < P，验证未充分利用情况
   - PE利用率低的场景
   - 验证idle PE的正确处理
   - 测试部分阵列激活逻辑

3. 等于阵列：M=K=N=P，验证完美匹配
   - 理想情况，100% PE利用率
   - 验证满负载下的功能正确性
   - 性能基准测试点

4. 轻微超出：M,K,N = P+1，验证最小分块
   - 触发分块逻辑
   - 验证块间切换开销
   - 测试padding处理

5. 2的幂次：N = {1,2,4,8,16,32,64,128,256}
   - 地址计算简单，便于优化
   - 常见的实际应用维度
   - 验证位操作优化的正确性

6. 质数维度：N = {13,17,23,31,37}，最难对齐
   - 最大化padding开销
   - 挑战地址生成逻辑
   - 验证通用性

7. 实际层维度：来自ResNet、BERT的真实层参数
   - ResNet-50: conv1(64×64×3×7×7), fc(1000×2048)
   - BERT-Base: attention(768×768), ffn(3072×768)
   - 验证实际应用场景

A[i][j] = i * 1000 + j  // 行列编码，便于追踪
B[i][j] = (i == j) ? 1 : 0  // 单位矩阵
期望：C[i][j] = i * 1000 + j  // 易于验证

Pattern A: [1,0,1,0,...]
Pattern B: [0,1,0,1,...]
验证相邻PE之间无数据污染

时刻T=0: 开始权重加载
时刻T=P: 权重就绪，开始数据注入
时刻T=2P-1: 首个输出出现在(0,0)位置
时刻T=2P: 第二个输出出现在(0,1)和(1,0)
时刻T=3P-2: 对角线输出达到稳态
时刻T=M+K+N-2: 最后一个输入进入
时刻T=M+K+N+2P-3: 最后一个输出产生

基础维度约束：
1 ≤ M, K, N ≤ 4096  // 覆盖实际应用范围
权重分布：

- 70%: 常见维度 [32, 64, 128, 256, 512, 1024]
- 20%: 边界情况 [1, P-1, P, P+1, 2P-1, 2P, 2P+1]
- 10%: 随机维度，包括质数

对齐约束（以P=32为例）：
M % 32 的分布：

- 40%: 0 (完美对齐)
- 20%: 1 (最小非对齐)
- 20%: 31 (最大非对齐)
- 10%: 16 (半对齐)
- 10%: 其他随机值

Layer 1: 维度组合 (M, K, N)
Layer 2: 数据分布 (uniform, gaussian, sparse)
Layer 3: 数值范围 (full range, small values, boundary values)
Layer 4: 执行模式 (continuous, interrupted, pipelined)

使用SystemVerilog约束：
constraint matrix_dims {
  // 基础约束
  M inside {[1:4096]};
  K inside {[1:4096]};
  N inside {[1:4096]};

  // 相关性约束
  (M > 1000) -> (K < 100);  // 大M配小K，测试极端长宽比

  // 分布约束
  M dist {
    32 := 10,
    64 := 10,
    128 := 20,
    256 := 20,
    512 := 15,
    1024 := 15,
    [1:31] := 5,
    [33:63] := 5
  };
}

covergroup matrix_dims_cg;
  M_cp: coverpoint M {
    bins small = {[1:31]};
    bins aligned = {32, 64, 96, 128};
    bins large = {[256:4096]};
  }
  K_cp: coverpoint K {
    bins small = {[1:31]};
    bins medium = {[32:255]};
    bins large = {[256:4096]};
  }
  N_cp: coverpoint N {
    bins small = {[1:31]};
    bins medium = {[32:255]};
    bins large = {[256:4096]};
  }
  // 三维交叉
  cross M_cp, K_cp, N_cp;
endgroup

property no_x_propagation;
  @(posedge clk) disable iff (!rst_n)
  !$isunknown(pe_output);
endproperty

assert property(no_x_propagation)
  else $error("X propagation detected");

cover property(no_x_propagation);  // 覆盖断言触发

┌─────────────────────────────────────────────────────────┐
│                    Test Environment                      │
├─────────────────────────────────────────────────────────┤
│  ┌─────────────┐                    ┌────────────────┐  │
│  │  Test Case  │                    │  Config Object │  │
│  └──────┬──────┘                    └────────────────┘  │
│         │                                               │
│  ┌──────▼────────────────────────────────────────────┐  │
│  │                    ENV (Environment)               │  │
│  ├────────────────────────────────────────────────────┤  │
│  │  ┌────────────┐  ┌────────────┐  ┌─────────────┐ │  │
│  │  │  Agent_In  │  │ Agent_Out  │  │  Scoreboard │ │  │
│  │  ├────────────┤  ├────────────┤  └─────────────┘ │  │
│  │  │ Sequencer  │  │  Monitor   │                   │  │
│  │  │   Driver   │  │            │   ┌─────────────┐ │  │
│  │  │  Monitor   │  └────────────┘   │  Coverage   │ │  │
│  │  └────────────┘                   └─────────────┘ │  │
│  └────────────────────────────────────────────────────┘  │
└─────────────────────────────────────────────────────────┘
                            │
                    ┌───────▼────────┐
                    │      DUT        │
                    │ (Systolic Array)│
                    └────────────────┘

class gemm_sequence extends uvm_sequence#(gemm_transaction);
  rand int unsigned m, k, n;
  rand data_pattern_e pattern;

  constraint dims_c {
    m inside {[1:1024]};
    k inside {[1:1024]};
    n inside {[1:1024]};
    // 权重分布控制
    m dist {32:=10, 64:=20, 128:=30, 256:=20, [1:31]:=10, [257:1024]:=10};
  }

  task body();
    gemm_transaction tr;
    tr = gemm_transaction::type_id::create("tr");

    // 配置事务
    start_item(tr);
    assert(tr.randomize() with {
      tr.m == local::m;
      tr.k == local::k;
      tr.n == local::n;
      tr.pattern == local::pattern;
    });
    finish_item(tr);
  endtask
endclass

class systolic_driver extends uvm_driver#(gemm_transaction);
  virtual systolic_if vif;

  task run_phase(uvm_phase phase);
    forever begin
      seq_item_port.get_next_item(req);
      drive_transaction(req);
      seq_item_port.item_done();
    end
  endtask

  task drive_transaction(gemm_transaction tr);
    // 配置阶段
    vif.config_m <= tr.m;
    vif.config_k <= tr.k;
    vif.config_n <= tr.n;
    vif.config_valid <= 1'b1;
    @(posedge vif.clk);
    vif.config_valid <= 1'b0;

    // 权重加载
    for(int i = 0; i < tr.k; i++) begin
      for(int j = 0; j < tr.n; j++) begin
        vif.weight_data <= tr.weight_matrix[i][j];
        vif.weight_valid <= 1'b1;
        @(posedge vif.clk);
      end
    end
    vif.weight_valid <= 1'b0;

    // 输入数据注入（斜向注入模式）
    fork
      inject_input_data(tr);
      collect_output_data(tr);
    join
  endtask
endclass

class systolic_monitor extends uvm_monitor;
  virtual systolic_if vif;
  uvm_analysis_port#(output_transaction) ap;

  task run_phase(uvm_phase phase);
    forever begin
      output_transaction tr;
      collect_output(tr);
      ap.write(tr);  // 广播给scoreboard
    end
  endtask

  task collect_output(output output_transaction tr);
    @(posedge vif.clk iff vif.output_valid);
    tr.row = vif.output_row;
    tr.col = vif.output_col;
    tr.data = vif.output_data;
    tr.timestamp = $time;
  endtask
endclass

class systolic_scoreboard extends uvm_scoreboard;
  uvm_tlm_analysis_fifo#(gemm_transaction) input_fifo;
  uvm_tlm_analysis_fifo#(output_transaction) output_fifo;

  // 参考模型实例
  systolic_reference_model ref_model;

  task run_phase(uvm_phase phase);
    forever begin
      gemm_transaction in_tr;
      output_transaction out_tr;

      input_fifo.get(in_tr);

      // 调用参考模型
      ref_model.compute(in_tr);

      // 收集所有输出并比对
      repeat(in_tr.m * in_tr.n) begin
        output_fifo.get(out_tr);
        check_result(out_tr, ref_model.get_expected(out_tr.row, out_tr.col));
      end
    end
  endtask

  function void check_result(output_transaction out_tr, real expected);
    real error_margin = 0.001;  // 容错范围

    if(abs(out_tr.data - expected) > error_margin) begin
      `uvm_error("MISMATCH", 
        $sformatf("Output[%0d][%0d]: Expected %f, Got %f", 
          out_tr.row, out_tr.col, expected, out_tr.data))
    end else begin
      `uvm_info("MATCH", 
        $sformatf("Output[%0d][%0d] correct: %f", 
          out_tr.row, out_tr.col, out_tr.data), UVM_HIGH)
    end
  endfunction
endclass

class SystolicRefModel {
private:
  int array_size;
  bool weight_stationary;
  float *weight_buffer;

public:
  void configure(int m, int k, int n, int p) {
    this->m = m; this->k = k; this->n = n;
    this->array_size = p;
  }

  void compute_gemm(float* A, float* B, float* C) {
    // 分块计算，模拟硬件行为
    int m_tiles = (m + array_size - 1) / array_size;
    int n_tiles = (n + array_size - 1) / array_size;
    int k_tiles = (k + array_size - 1) / array_size;

    for(int mt = 0; mt < m_tiles; mt++) {
      for(int nt = 0; nt < n_tiles; nt++) {
        for(int kt = 0; kt < k_tiles; kt++) {
          compute_tile(A, B, C, mt, nt, kt);
        }
      }
    }
  }

  void compute_tile(float* A, float* B, float* C, 
                   int mt, int nt, int kt) {
    int m_start = mt * array_size;
    int n_start = nt * array_size;
    int k_start = kt * array_size;

    int m_end = min(m_start + array_size, m);
    int n_end = min(n_start + array_size, n);
    int k_end = min(k_start + array_size, k);

    // 精确模拟脉动阵列计算顺序
    for(int i = m_start; i < m_end; i++) {
      for(int j = n_start; j < n_end; j++) {
        for(int l = k_start; l < k_end; l++) {
          C[i*n + j] += A[i*k + l] * B[l*n + j];
        }
      }
    }
  }
};

import numpy as np

class SystolicReference:
    def __init__(self, array_size=32, dtype='float16'):
        self.array_size = array_size
        self.dtype = dtype

    def compute(self, A, B, weight_stationary=True):
        """
        计算矩阵乘法 C = A @ B
        A: [M, K], B: [K, N]
        """
        M, K = A.shape
        K2, N = B.shape
        assert K == K2, "矩阵维度不匹配"

        # 转换数据类型模拟硬件
        if self.dtype == 'nvfp4':
            A = self.quantize_nvfp4(A)
            B = self.quantize_nvfp4(B)

        # 分块计算
        C = np.zeros((M, N), dtype=A.dtype)
        P = self.array_size

        for m in range(0, M, P):
            for n in range(0, N, P):
                for k in range(0, K, P):
                    # 提取tile
                    m_end = min(m + P, M)
                    n_end = min(n + P, N)
                    k_end = min(k + P, K)

                    A_tile = A[m:m_end, k:k_end]
                    B_tile = B[k:k_end, n:n_end]

                    # 计算并累加
                    C[m:m_end, n:n_end] += self.systolic_compute(
                        A_tile, B_tile, weight_stationary
                    )

        return C

    def systolic_compute(self, A_tile, B_tile, weight_stationary):
        """模拟单个tile的脉动阵列计算"""
        if weight_stationary:
            # Weight-stationary数据流
            return self.ws_dataflow(A_tile, B_tile)
        else:
            # Output-stationary数据流
            return self.os_dataflow(A_tile, B_tile)

    def ws_dataflow(self, A, B):
        """Weight-stationary精确时序模拟"""
        M, K = A.shape
        K2, N = B.shape
        C = np.zeros((M, N))

        # 模拟逐周期计算
        for cycle in range(M + K + N - 2):
            for i in range(M):
                for j in range(N):
                    # 计算数据到达时间
                    k_idx = cycle - i - j
                    if 0 <= k_idx < K:
                        C[i, j] += A[i, k_idx] * B[k_idx, j]

        return C

    def quantize_nvfp4(self, x, bias=1):
        """nvfp4量化模拟"""
        # E2M1格式：1位符号，2位指数，1位尾数
        sign = np.sign(x)
        abs_x = np.abs(x)

        # 计算指数
        exp = np.floor(np.log2(abs_x)) + bias
        exp = np.clip(exp, 0, 3)  # 2位指数

        # 计算尾数
        mantissa = abs_x / (2 ** (exp - bias)) - 1
        mantissa = np.round(mantissa * 2) / 2  # 1位尾数

        # 重构数值
        return sign * (1 + mantissa) * (2 ** (exp - bias))

// DPI-C函数声明
import "DPI-C" function void c_ref_model_init(int array_size);
import "DPI-C" function void c_ref_model_compute(
    input real A[], input real B[], 
    output real C[], 
    input int m, k, n
);

class dpi_reference_model extends uvm_object;
  function new(string name = "dpi_reference_model");
    super.new(name);
    c_ref_model_init(32);  // 初始化32x32阵列
  endfunction

  function void compute(gemm_transaction tr, ref real result[]);
    real A[], B[], C[];

    // 展平矩阵数据
    A = new[tr.m * tr.k];
    B = new[tr.k * tr.n];
    C = new[tr.m * tr.n];

    flatten_matrix(tr.A_matrix, A);
    flatten_matrix(tr.B_matrix, B);

    // 调用C++模型
    c_ref_model_compute(A, B, C, tr.m, tr.k, tr.n);

    // 重组结果
    unflatten_matrix(C, result, tr.m, tr.n);
  endfunction
endclass

┌─────────────────────────────────────────────────┐
│          Cycle-Accurate Simulator               │
├─────────────────────────────────────────────────┤
│  ┌──────────────┐    ┌──────────────┐         │
│  │  Compute     │    │   Memory     │         │
│  │  Model       │    │   Model      │         │
│  │  - PE Array  │    │  - SRAM      │         │
│  │  - Pipeline  │    │  - DRAM      │         │
│  │  - Dataflow  │    │  - NoC       │         │
│  └──────┬───────┘    └──────┬───────┘         │
│         │                    │                  │
│  ┌──────▼────────────────────▼──────┐         │
│  │       Timing Model                │         │
│  │  - Clock domains                  │         │
│  │  - Synchronization                │         │
│  │  - Pipeline stages                │         │
│  └───────────────┬───────────────────┘         │
│                  │                              │
│  ┌───────────────▼───────────────────┐         │
│  │     Performance Statistics        │         │
│  │  - Cycle counts                   │         │
│  │  - Utilization                    │         │
│  │  - Bottleneck analysis            │         │
│  └───────────────────────────────────┘         │
└─────────────────────────────────────────────────┘

class PEArrayModel {
private:
    int array_size;
    int pipeline_depth;
    vector<vector<PEState>> pe_states;

public:
    struct PEState {
        bool is_active;
        int current_cycle;
        float accumulator;
        float weight_reg;
        float input_reg;
        float output_reg;
    };

    void cycle_update() {
        // 更新每个PE的状态
        for(int i = 0; i < array_size; i++) {
            for(int j = 0; j < array_size; j++) {
                PEState& pe = pe_states[i][j];

                if(pe.is_active) {
                    // MAC运算
                    pe.accumulator += pe.weight_reg * pe.input_reg;

                    // 数据传递（向右和向下）
                    if(j < array_size - 1) {
                        pe_states[i][j+1].input_reg = pe.input_reg;
                    }
                    if(i < array_size - 1) {
                        pe_states[i+1][j].output_reg = pe.accumulator;
                    }
                }
            }
        }

        current_cycle++;
    }

    int compute_latency(int M, int K, int N) {
        // 计算总延迟
        int num_m_tiles = (M + array_size - 1) / array_size;
        int num_n_tiles = (N + array_size - 1) / array_size;
        int num_k_tiles = (K + array_size - 1) / array_size;

        int cycles_per_tile = array_size * num_k_tiles;
        int pipeline_fill = 2 * array_size - 1;
        int pipeline_drain = array_size - 1;

        int total_cycles = num_m_tiles * num_n_tiles * 
                          (cycles_per_tile + pipeline_fill + pipeline_drain);

        return total_cycles;
    }
};

流水线阶段分析：
Stage 1: 指令译码 (1 cycle)
Stage 2: 地址计算 (1 cycle)  
Stage 3: 存储读取 (2-3 cycles, 取决于bank冲突)
Stage 4: 数据对齐 (1 cycle)
Stage 5: PE计算 (1 cycle MAC)
Stage 6: 累加更新 (1 cycle)
Stage 7: 结果写回 (2-3 cycles)

总流水线深度: 9-11 cycles

class MemoryHierarchy {
private:
    struct CacheLevel {
        int size_kb;
        int latency_cycles;
        int bandwidth_gbps;
        float hit_rate;
    };

    vector<CacheLevel> levels = {
        {32,    1,   2048, 0.95},  // L0: PE本地寄存器
        {256,   3,   1024, 0.85},  // L1: Tile SRAM
        {4096,  10,  512,  0.75},  // L2: Global SRAM
        {32768, 100, 256,  1.0}    // L3: HBM/DDR
    };

public:
    int access_latency(int address, int size) {
        int total_latency = 0;

        for(auto& level : levels) {
            if(random() < level.hit_rate) {
                // 命中当前级别
                int transfer_cycles = (size * 8) / 
                    (level.bandwidth_gbps * 1e9 / CLOCK_FREQ);
                total_latency = level.latency_cycles + transfer_cycles;
                break;
            }
            // 未命中，继续下一级
            total_latency += level.latency_cycles;
        }

        return total_latency;
    }

    void model_prefetch(int stride, int count) {
        // 预取建模
        for(int i = 0; i < count; i++) {
            int addr = base_addr + i * stride;
            prefetch_queue.push(addr);
        }
    }
};

class SRAMBankModel {
    static const int NUM_BANKS = 16;
    static const int BANK_WIDTH = 256;  // bits

    struct BankState {
        bool is_busy;
        int busy_until_cycle;
        queue<AccessRequest> pending_requests;
    };

    BankState banks[NUM_BANKS];

    int schedule_access(int address, int cycle) {
        int bank_id = (address / BANK_WIDTH) % NUM_BANKS;
        BankState& bank = banks[bank_id];

        if(bank.busy_until_cycle <= cycle) {
            // 无冲突
            bank.busy_until_cycle = cycle + 1;
            return cycle;
        } else {
            // Bank冲突，需要等待
            int actual_cycle = bank.busy_until_cycle;
            bank.busy_until_cycle = actual_cycle + 1;
            conflict_count++;
            return actual_cycle;
        }
    }
};

优点：权重复用最大化，减少权重加载开销
缺点：输入/输出数据流动开销大

性能模型：
T_ws = T_weight_load + T_compute + T_output_collect
其中：

- T_weight_load = P × (K/P) = K cycles
- T_compute = M × N / P² × K cycles  
- T_output_collect = M × N cycles

优点：部分和保持在PE本地，减少累加开销
缺点：权重和输入都需要流动

性能模型：
T_os = T_init + T_compute + T_writeback
其中：

- T_init = M × N / P² cycles
- T_compute = K × max(M/P, N/P) cycles
- T_writeback = M × N / P² cycles

优点：平衡各种数据复用
缺点：控制复杂度高

性能模型需要考虑行级数据驻留：
T_rs = Σ(T_row_compute[i]) for i in [0, M/P)

PE利用率：
η_PE = 实际活跃PE周期数 / (P² × 总周期数)

对于非对齐矩阵：

- M=100, K=100, N=100, P=32
- 需要4×4×4=64个tiles
- 每个tile利用率：(100%32)²/32² = 4²/32² = 1.56%
- 整体利用率很低！

Event ID | Event Type           | Counter Value
---------|---------------------|---------------
0x01     | GEMM_START         | Timestamp
0x02     | WEIGHT_LOAD_BEGIN  | Cycle count
0x03     | WEIGHT_LOAD_END    | Cycle count
0x04     | COMPUTE_BEGIN      | Cycle count
0x05     | FIRST_OUTPUT       | Cycle count
0x06     | COMPUTE_END        | Cycle count
0x07     | MEMORY_STALL       | Stall cycles
0x08     | BANK_CONFLICT      | Conflict count

符号位(S) | 指数(E1E0) | 尾数(M0)
    1     |     2      |    1

最大累加次数（nvfp4）：
N_max = 2^24 / max_value = 2^24 / 6 ≈ 2.8M

测试矩阵：

- M=17, K=33, N=65：全部需要padding
- M=16, K=31, N=16：仅K维需要padding
- M=1, K=1, N=1：最小矩阵
- M=15, K=15, N=15：接近但不等于阵列大小