CUDA教程9 -- 代数多重网格

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108

void AMG::Selector(int level) {
    const CsrMatrix &A = level == 0 ? _A : _coarseA[level - 1];
    
    // 节点个数
    const IndexType N = A.num_rows;

    // 初始化为0，即未标记
    CHECK(cudaMalloc(&_nodes_prop[level], N * sizeof(int)));
    CHECK(cudaMemset(_nodes_prop[level], 0, N * sizeof(int)));

    // 为每个节点生成一个与它影响的节点个数相关的权重
    float *weight = nullptr;
    CHECK(cudaMalloc(&weight, N * sizeof(float)));

    curandGenerator_t gen;
    CHECK_CURAND(curandCreateGenerator(
            &gen, CURAND_RNG_PSEUDO_DEFAULT
    ));

    // 选择一个随机数种子然后生成随机数，由于不需要真正的伪随机数，只是为了
    // 给度相同的节点赋予不同的权重，所以种子不需要随机性
    CHECK_CURAND(curandSetPseudoRandomGeneratorSeed(gen, 1234ULL));
    CHECK_CURAND(curandGenerateUniform(gen, weight, N));

    // 计算每个节点的强影响出度，加到权重上。即#Si^T = #{j | Sji = 1}
    // Sij = 1, iff i!=j && -Aij > alpha * max_{k!=i} { -Aik }。
    // Sij=1表示i强依赖于j或者j强影响i，Si^T就是i强影响的点集。
    // 遍历矩阵，如果-Aij满足上述要求，则将weight[j]加1
    LaunchKernelWithNumThreads(
            PMISCalWeightKernel, N, 0, 0,
            weight, A.values, A.row_offsets, A.column_indices,
            _thresholds[level].data(), N
    );

    // 使用归约来计算有多少节点被标记了
    IndexType num_marked = 0;

    int max_grid = 0;
    int block = 0;
    LaunchConfig(&max_grid, &block);

    int grid = (N + block - 1) / block;
    if (grid > max_grid)
        grid = max_grid;

    IndexType *counts = nullptr;
    CHECK(cudaMalloc(&counts, grid * sizeof(IndexType)));

    std::vector<IndexType> h_counts(grid);

    // 首先将#Si^T为0的的节点标记为细网格节点
    LaunchKernelWithNumThreads(
            PMISMarkFineIfSiTEqualZeroKernel<>, N, 0, 0,
            _nodes_prop[level], counts, weight, N
    );

    CHECK(cudaMemcpy(h_counts.data(), counts,
                     grid * sizeof(IndexType), cudaMemcpyDeviceToHost));

    for (IndexType n : h_counts)
        num_marked += n;


    // 将拥有局部最大权重的节点标记为粗网格节点，将它的邻居标记为细网格节点。
    // 在每一轮，为消除不同线程间标记的影响，使每轮的标记逐次增大，奇数为粗网格，
    // 偶数为细网格，全部标记后再将它们转换成1(粗)和2(细)
    CsrMatrix At = A.transpose();

    for (IndexType mark = 1; mark <= N; ++mark) {
        LaunchKernelWithNumThreads(
                PMISSelectorMarkCoarseKernel, N, 0, 0,
                _nodes_prop[level],
                A.values, A.row_offsets, A.column_indices,
                _thresholds[level].data(), weight, N, mark
        );

        LaunchKernelWithNumThreads(
                PMISCountMarkedKernel<>, N, 0, 0,
                counts, _nodes_prop[level], N
        );

        CHECK(cudaMemcpy(h_counts.data(), counts,
                         grid * sizeof(IndexType), cudaMemcpyDeviceToHost));

        num_marked = 0;
        for (IndexType n : h_counts)
            num_marked += n;

        if (num_marked == N)
            break;
    }

    // 计算有多少个粗网格和细网格，并将标记转换成1(粗)和2(细)
    LaunchKernelWithNumThreads(
            PMISCountCoarseKernel<>, N, 0, 0,
            counts, _nodes_prop[level], N
    );

    CHECK(cudaMemcpy(h_counts.data(), counts,
                     grid * sizeof(IndexType), cudaMemcpyDeviceToHost));

    _num_coarse[level] = 0;
    for (IndexType n : h_counts)
        _num_coarse[level] += n;

    CHECK(cudaFree(weight));
    CHECK(cudaFree(counts));
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

// 一个线程处理一行
__global__ void PMISSelectorMarkCoarseKernel(
        IndexType *nodes_prop,
        const ValueType *__restrict__ A_values,
        const IndexType *__restrict__ A_row_offsets,
        const IndexType *__restrict__ A_column_indices,
        const ValueType *__restrict__ threshold,
        const float *__restrict__ weight,
        const IndexType num_rows,
        const IndexType mark
) {
    const IndexType odd = 2 * mark - 1;
    const IndexType even = 2 * mark;

    IndexType row = threadIdx.x + blockIdx.x * blockDim.x;
    const unsigned int stride = blockDim.x * gridDim.x;
    for (; row < num_rows; row += stride) {
        // 已经标记的点跳过
        if (nodes_prop[row] != 0)
            continue;

        // 计算row这个节点是否为局部最大权重的点
        bool is_local_max = true;
        for (IndexType j = A_row_offsets[row]; j < A_row_offsets[row + 1]; ++j) {
            IndexType col = A_column_indices[j];

            // 自身，已经标记的(不能是本轮标记的)，非强影响i且i不强影响它的，都不考虑
            if (row == col
                || (nodes_prop[col] != 0           // 读
                    && nodes_prop[col] != even     // 读，旧值新值条件判断结果一样
                    && nodes_prop[col] != odd)     // 读，弱一致性即可
                || (-A_values[j] < threshold[row]
                    && -A_values[j] < threshold[col]))
                continue;

            // 只要有一个邻居的权重比它大，它就不是局部最大节点
            if (weight[row] < weight[col]) {
                is_local_max = false;
                break;
            }
        }

        // 如果它是局部最大节点，将它标记为粗网格，将它的未标记的邻居标记为细网格
        if (is_local_max) {
            nodes_prop[row] = odd;    // 写

            for (IndexType j = A_row_offsets[row]; j < A_row_offsets[row + 1]; ++j) {
                IndexType col = A_column_indices[j];

                if (row != col && nodes_prop[col] == 0
                    && (-A_values[j] >= threshold[row]
                        || -A_values[j] >= threshold[col]))
                    nodes_prop[col] = even;   // 写
            }
        }
    }
}