wjin

CUDA教程7 -- 双共轭梯度法

 5243  11 分钟

§稀疏矩阵的表示

以下面的矩阵为例：

$$ A=\begin{bmatrix} 1 & 2 & 0 & 0 \\ 0 & 3 & 4 & 5 \\ 0 & 6 & 7 & 0 \\ 0 & 0 & 8 & 9 \end{bmatrix} $$

COO(Coordinate)格式将它存储在三个数组中：

row: 0  0  1  1  1  2  2  3  3
col: 0  1  1  2  3  1  2  2  3
val: 1  2  3  4  5  6  7  8  9

将行号进行压缩可以得到CSR(Compressed Sparse Row)格式：

row: 0  2  5  7  9
col: 0  1  1  2  3  1  2  2  3
val: 1  2  3  4  5  6  7  8  9

相邻的两个元素标识了矩阵一行的非零元素在col和val中的下标范围[row[i], row[i+1])。为保证最后一行也不越界，一般在row末尾再添加一个哨兵元素，它的值是矩阵非零元素个数，如上面的9。

§矩阵转置

§稠密矩阵

稠密矩阵的转置比较简单，将上三角的元素与下三角逐一交换即可：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

void Transpose(double *A, int M, int N, int lda = N) {
    // MxN的矩阵，一行需要跳过lda个元素
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < i; ++j) {
            // swap(A[i, j], A[j, i])
            std::swap(A[i * lda + j], A[j * lda + i]);
        }
    }
}

void Transpose(double *A, double *At, int M, int N, int lda = N, int ldat = M) {
    // MxN的矩阵A，一行需要跳过lda个元素
    // 输出到NxM的矩阵At,一行需要跳过ldat个元素
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < i; ++j) {
            At[j * ldat + i] = A[i * lda + j];
        }
    }
}

§COO格式

稀疏矩阵的转置就不太容易了。对于COO格式的矩阵，如果算法对元素存储顺序没有要求，那么可以简单地交换每个元素的行号和列号，如矩阵

$$ A=\begin{bmatrix} 1 & 2 & 0 & 0 \\ 0 & 3 & 4 & 5 \\ 0 & 6 & 7 & 0 \\ 0 & 0 & 8 & 9 \end{bmatrix} $$

的转置为：

$$ A^T=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 2 & 3 & 6 & 0 \\ 0 & 4 & 7 & 8 \\ 0 & 5 & 0 & 9 \end{bmatrix} $$
它们的COO格式分别表示为： 
A:
row: 0  0  1  1  1  2  2  3  3
col: 0  1  1  2  3  1  2  2  3
val: 1  2  3  4  5  6  7  8  9

At:
row: 0  1  1  2  3  1  2  2  3
col: 0  0  1  1  1  2  2  3  3
val: 1  2  3  4  5  6  7  8  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

void Transpose(IndexType *rows, IndexType *cols, const IndexType num_entries) {
    for (IndexType i = 0; i < num_entries; ++i)
        std::swap(rows[i], cols[i]);
}

// COO格式矩阵的转置乘向量
void SpMTvCOO(
        Vector *result,
        const IndexType *rows,
        const IndexType *cols,
        const ValueType *vals,
        const IndexType num_rows,
        const IndexType num_cols,
        const IndexType num_entries,
        const Vector &v
) {
    if (num_cols != v.size)
        throw std::runtime_error{"SpMTvCOO: 矩阵的列数和向量长度不一致！"};
    
    // 仅交换指针
    std::swap(rows, cols);
    
    result->resize(num_cols);
    for (IndexType i = 0; i < num_cols; ++i)
        (*result)[i] = 0;
    
    for (IndexType i = 0; i < num_entries; ++i)
        (*result)[rows[i]] += vals[i] * v[cols[i]];
}

§CSR格式

对于CSR格式就不能简单地交换行号和列号了，可以使用按行填充的办法来进行转置。算法的基本想法是：按行遍历矩阵，对每个非零元素，可以得知其行号和列号，列号为它在转置矩阵中的行号，行号为它在转置矩阵中的列号，把非零元素逐个填充到转置矩阵中对应行的位置。此处有两个问题：1、转置矩阵中某行的起始下标是多少；2、如何保证填充过去的元素有序。

第一，转置矩阵中各行的起始下标在数值上等于它前面行的非零元素个数之和，因此可以这样计算：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

CsrMatrix CsrMatrix::transpose() const {
    CsrMatrix At(num_cols, num_rows, num_entries);
    std::fill(At.row_offsets.begin(), At.row_offsets.end(), 0);

    // 计算A^T每行有多少个非零元素
    for (IndexType i = 0; i < num_entries; ++i)
        ++(At.row_offsets[column_indices[i]]);

    // 计算A^T每行首个非零元素偏移量
    for (IndexType i = 0, sum = 0; i < At.num_rows; ++i) {
        IndexType temp = At.row_offsets[i];
        At.row_offsets[i] = sum;
        sum += temp;
    }
    At.row_offsets[At.num_rows] = At.num_entries;
    
    // ...
}

第二，转置矩阵的元素顺序需要先按行号排序，行内按列号排序，遍历原矩阵时，从第0行逐行遍历到最后一行，则填充到转置矩阵时列号就是逐渐增大的，因此只需要按行主序遍历原矩阵，得到的转置矩阵就已经满足顺序要求了。完整的算法如下：

 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

namespace Host {
    
CsrMatrix CsrMatrix::transpose() const {
    CsrMatrix At(num_cols, num_rows, num_entries);
    std::fill(At.row_offsets.begin(), At.row_offsets.end(), 0);

    // 计算A^T每行有多少个非零元素
    for (IndexType i = 0; i < num_entries; ++i)
        ++(At.row_offsets[column_indices[i]]);

    // 计算A^T每行首个非零元素偏移量
    for (IndexType i = 0, sum = 0; i < At.num_rows; ++i) {
        IndexType temp = At.row_offsets[i];
        At.row_offsets[i] = sum;
        sum += temp;
    }
    At.row_offsets[At.num_rows] = At.num_entries;

    // 使用curr记录各行已填充到哪个位置
    std::vector<IndexType> curr = At.row_offsets;

    // 将非零元素复制到A^T中
    for (IndexType i = 0; i < num_rows; ++i) {
        for (IndexType j = row_offsets[i]; j < row_offsets[i + 1]; ++j) {
            IndexType col = column_indices[j];
            At.values[curr[col]] = values[j];
            At.column_indices[curr[col]] = i;
            ++(curr[col]);
        }
    }

    return At;
}
    
} // namespace Host

仍然以矩阵

$$ A=\begin{bmatrix} 1 & 2 & 0 & 0 \\ 0 & 3 & 4 & 5 \\ 0 & 6 & 7 & 0 \\ 0 & 0 & 8 & 9 \end{bmatrix} $$

为例，其CSR格式存储为：

row: 0  2  5  7  9
col: 0  1  1  2  3  1  2  2  3
val: 1  2  3  4  5  6  7  8  9

根据列号计算得到转置的每行元素个数为：

rowT: 1 3 3 2

对它累和得到转置矩阵的row:

rowT: 0 1 4 7 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

ind: 0  1  2  3  4  5  6  7  8
col: _  _        _        _
val: _  _        _        _
    
col: 0  _        _        _
val: 1  _        _        _
    
col: 0  0        _        _
val: 1  2        _        _
    
col: 0  0  1     _        _
val: 1  2  3     _        _
    
col: 0  0  1     1        _
val: 1  2  3     4        _
    
col: 0  0  1     1        1
val: 1  2  3     4        5
    
col: 0  0  1  2  1        1
val: 1  2  3  6  4        5
    
col: 0  0  1  2  1  2     1
val: 1  2  3  6  4  7     5
    
col: 0  0  1  2  1  2  3  1
val: 1  2  3  6  4  7  8  5
    
col: 0  0  1  2  1  2  3  1  3
val: 1  2  3  6  4  7  8  5  9

其中下划线是各行第一个非零元素所在位置，其它未填充的位置留白。

这种算法是串行算法中工作得最好的，因为它只需要遍历矩阵一次，时间复杂度为O(M+NNZ)，其中M为矩阵行数，NNZ为非零元素个数。

另一个简单的思路是将CSR格式中压缩的行还原成COO格式，交换行号和列号，再按行号排序、行内按列号排序，再将它压缩成稀疏矩阵，这个算法的时间复杂度为O(NNZ log NNZ)，明显高于前一个算法，它的好处在于可以更方便地实现为并行算法。此处给出使用CUDA设计的算法：

  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

namespace SingleGPU {
    
__global__ void Csr2CooKernel(
        IndexType *__restrict__ row_indices,
        const IndexType *__restrict__ row_offsets,
        const IndexType num_rows
) {
    const int thread_id = threadIdx.x + blockIdx.x * blockDim.x;
    const int lane = thread_id & 15;
    const int num_rows_per_block = (blockDim.x * gridDim.x) / 16;

    for (IndexType row = thread_id / 16; row < num_rows; row += num_rows_per_block) {
        for (IndexType j = row_offsets[row] + lane; j < row_offsets[row + 1]; j += 16)
            row_indices[j] = row;
    }
}
    
__global__ void Coo2CsrKernel(
        IndexType *__restrict__ row_offsets,
        const IndexType *__restrict__ sorted_row_indices,
        const IndexType num_rows,
        const IndexType num_entries
) {
    IndexType i = threadIdx.x + blockIdx.x * blockDim.x;

    if (i == 0) {
        row_offsets[0] = 0;

        if (num_entries > 0) {
            for (IndexType j = 1; j <= sorted_row_indices[0]; ++j)
                row_offsets[j] = 0;
        }
    }

    const int stride = blockDim.x * gridDim.x;
    for (; i < num_entries; i += stride) {
        IndexType start = sorted_row_indices[i] + 1;
        IndexType end = i < num_entries - 1 ? sorted_row_indices[i + 1] : num_rows;
        for (IndexType j = start; j <= end; ++j)
            row_offsets[j] = i + 1;
    }
}

CsrMatrix CsrMatrix::transpose() const {
    // csr转为coo，只需要将row_offsets转为row_indices
    IndexType *row_indices = nullptr;
    CHECK(cudaMalloc(&row_indices, sizeof(IndexType) * num_entries));
    LaunchKernel(Csr2CooKernel, 0, 0, row_indices, row_offsets, num_rows);

    // 交换row_indices和column_indices
    CsrMatrix At(num_cols, num_rows, num_entries);

    CHECK(cudaMemcpy(
            At.values, values, sizeof(ValueType) * num_entries,
            cudaMemcpyDeviceToDevice
    ));

    CHECK(cudaMemcpy(
            At.column_indices, row_indices, sizeof(IndexType) * num_entries,
            cudaMemcpyDeviceToDevice
    ));

    // 按照行排序，At的行即A的列
    IndexType *At_row_indices = nullptr;
    CHECK(cudaMalloc(&At_row_indices, sizeof(IndexType) * num_entries));

    CHECK(cudaMemcpy(
            At_row_indices, column_indices, sizeof(IndexType) * num_entries,
            cudaMemcpyDeviceToDevice
    ));
    
    thrust::stable_sort_by_key(
            thrust::device_ptr<IndexType>(At_row_indices),
            thrust::device_ptr<IndexType>(At_row_indices + num_entries),
            thrust::device_ptr<ValueType>(At.values)
    );

    CHECK(cudaMemcpy(
            At_row_indices, column_indices, sizeof(IndexType) * num_entries,
            cudaMemcpyDeviceToDevice
    ));
    
    thrust::stable_sort_by_key(
            thrust::device_ptr<IndexType>(At_row_indices),
            thrust::device_ptr<IndexType>(At_row_indices + num_entries),
            thrust::device_ptr<IndexType>(At.column_indices)
    );

    // 将At_row_indices转换成At.row_offsets，由于At_row_indices已经排好序了，直接
    // 找各个数字第一次出现的位置即可。如0 0 0 1 1 2 2 2 => 0 3 5 8
    LaunchKernelWithNumThreads(
            Coo2CsrKernel, At.num_entries, 0, 0,
            At.row_offsets, At_row_indices, At.num_rows, At.num_entries
    );

    CHECK(cudaFree(row_indices));
    CHECK(cudaFree(At_row_indices));

    return At;
}
    
} // namespace SingleGPU

首先第一步是将CSR格式中压缩的行再展开成COO格式中的行号，这个非常简单，按行遍历矩阵，将[row[i], row[i+1])范围内的行号都置为i即可，如

row: 0 1 4 7 9

展开成

row: 0 1 1 1 2 2 2 3 3

这一步由核函数Csr2CooKernel来做。

其次是交换行号和列号，即将上一步展开的行号作为转置矩阵的列号，将原矩阵的列号作为行号。

然后是对转置矩阵的非零元素排序，先按行号排序，同行号按列号排序。由于排序之前，行号相同的元素其列号已经是从小到大排列的，只是不相邻而已，因此使用稳定的排序算法排序可以直接得到行内列号有序的结果。也可以使用自定义的比较函数进行排序，这个比较函数的比较方法为先比较行号，如果行号相等再比较列号，大致为：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

struct RowCol {
    IndexType row;
    IndexType col;
};

bool cmp(RowCol a, RowCol b) {
    if (a.row == b.row)
        return a.col < b.col;
    else
        return a.row < b.row;
}

std::vector<RowCol> coordinates;
std::vector<ValueType> values;
sort_by_key(coordinates.begin(), coordinates.end(), values.begin());

此处选用第一种方案，原因很简单：不想多引入一个结构体，在GPU上对结构体的操作不好封装，而且对结构体数组的访问性能一般不如数组结构体。

串行的排序算法我们已经很熟悉了，比如：快速排序、归并排序、堆排序等，时间复杂度都是O(n log n)。当然快速排序最坏情况下是O(n^2)，可以通过一些策略来改善其性能，比如GNU stl(SGI stl)的采用的方案为内省式排序(Introspective Sort)，由Musser在1997年[1]提出，GNU stl的实现细节可以参考[2]。

然后并行的排序算法基本都不能直接由快速排序等经典排序算法并行化而来，必须专门设计。比较常见的并行排序方法有：奇偶排序[3]、双调排序[4]、基数排序[5]等。实现这些算法太过复杂，此处调用thrust库完成排序。thrust库是一个类似于stl的模板库，已经集成在CUDA安装包中，默认会安装。其使用方法与stl类似，可以参考其文档[6]。thrust库的排序方法使用基数排序，拥有非常好的性能，此处直接调用。

排完序后，最后一步是将COO格式的行号压缩成CSR格式。一个简单地办法是统计各行的非零元素个数，然后进行扫描生成各行第一个非零元素的偏移：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

__global__ void CountNNZKernel(
        IndexType *nnz_per_rows,
        const IndexType *row,
        const IndexType num_entries
) {
    IndexType i = threadIdx.x + blockIdx.x * blockDim.x;
    const int stride = blockDim.x * gridDim.x;
    for (; i < num_entries; i += stride) {
        atomicAdd(&nnz_per_rows[row[i]], 1);
    }
}

CHECK(cudaMemset(At.row_offsets, 0, sizeof(IndexType) * (At.num_rows + 1)));

LaunchKernel(
        CountNNZKernel, 0, 0,
        At.row_offsets, column_indices, At.num_entries
);

thrust::exclusive_scan(
        thrust::device_ptr<IndexType>(At.row_offsets),
        thrust::device_ptr<IndexType>(At.row_offsets + At.num_rows + 1),
        thrust::device_ptr<IndexType>(At.row_offsets)
);

或者利用这样的性质：行号变化的位置就是对应行的第一个非零元素所在的位置。如

ind: 0 1 2 3 4 5 6 7 8
row: 0 1 1 1 2 2 2 3 3

行号在下标为1、4、7的位置变化，它们就是对应行的第一个非零元素，因此不难得到row_offsets为：0, 1, 4, 7, 9。这个过程由核函数Coo2CsrKernel完成。需要注意：这个算法虽然简单，但边界条件非常多，要写正确非常困难，比如一个稍微极端的例子：

ind: 0 1 2 3 4 5 6 7 8
row: 0 0 0 0 2 2 2 4 4

压缩后应为：

row: 0 4 4 7 7 9

而不是：

row: 0 4 7 9

这个矩阵有5行，压缩后的行号应当有6个元素，第1、3两行没有非零元素。此外还要考虑第1个0以及最后一个哨兵元素的生成，还有注意下标不要越界！

§双共轭梯度法

双共轭梯度法可以求解非对称、非正定的方程组，其思想在于同时求解$Ax=b$和$A^Tx=b$，详细过程如下：

$$ \begin{aligned} & r_0=b-Ax^{(0)} \\ & choose\ r_0^*:(r_0, r_0^*) \ne 0 \\ & p_1=r_0 \\ & p_1^*=r_0^* \\ & for\ k=1,2,\dots,n \\ &\qquad \xi_k=(r_{k-1},r_{k-1}^*)/(Ap_k,p_k^*) \\ &\qquad x^{(k)}=x^{(k-1)}+\xi_k p_k \\ &\qquad r_k=r_{k-1}-\xi_k Ap_k \\ &\qquad if\ \|r_k\| < error \\ &\qquad\qquad break \\ &\qquad end \\ &\qquad r_k^*=r_{k-1}^*-\xi_k A^Tp_k^* \\ &\qquad \theta_k=(r_k,r_k^*)/(r_{k-1},r_{k-1}^*) \\ &\qquad p_{k+1}=r_{k} + \theta_k p_k \\ &\qquad p_{k+1}^*=r_k^*+\theta_k p_k^* \\ & end \end{aligned} $$
其中除了涉及到共轭梯度法需要矩阵和向量基本操作外，还涉及到矩阵的转置乘向量。为实现矩阵的转置乘向量，一是可以直接用原矩阵进行计算，方法为用向量的各个元素作为系数对矩阵的各行作线性组合： 
 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

#if !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600
#else
static __inline__ __device__ double atomicAdd(
        double *address,
        double val
) {
    unsigned long long *address_as_ull = (unsigned long long *) address;
    unsigned long long old = *address_as_ull;
    unsigned long long expected;

    do {
        expected = old;
        old = atomicCAS(address_as_ull, expected,
                        __double_as_longlong(val + __longlong_as_double(expected)));
    } while (expected != old);

    return __longlong_as_double(old);
}
#endif

__global__ void SpMTvKernel(
        ValueType *__restrict__ result,
        const ValueType *__restrict__ values,
        const IndexType *__restrict__ row_offsets,
        const IndexType *__restrict__ column_indices,
        const IndexType num_rows,
        const ValueType *__restrict__ vec
) {
    const int thread_id = blockIdx.x * blockDim.x + threadIdx.x;
    const int lane = thread_id & 15;           // 线程束内偏移量
    const int num_warps = gridDim.x * blockDim.x / 16;

    for (int row = thread_id / 16; row < num_rows; row += num_warps) {
        for (int j = row_offsets[row] + lane; j < row_offsets[row + 1]; j += 16)
            atomicAdd(&result[column_indices[j]], vec[row] * values[j]);
    }
}

void SpMTv(Vector *result, const CsrMatrix &A, const Vector &v) {
    if (A.num_rows != v.size)
        throw std::runtime_error{
            "SingleGPU::SpMTv: 参与运算的矩阵行数与向量维度不一致！"
        };

    result->resize(A.num_cols);
    result->fill(0.0);
    LaunchKernel(
            SpMTvKernel, 0, 0,
            result->data(),
            A.values, A.row_offsets, A.column_indices, A.num_rows,
            v.data()
    );
}

由于计算能力6.0以下的设备无双精度浮点类型的原子操作，所以需要自己实现一个，可以参考CUDA文档Programming Guide的B.14节Atomic Functions。

这种方案的缺点是需要使用原子操作，在数万甚至数十万个线程并发时非常影响性能，而且这个操作在循环中，需要反复调用，很容易成为性能瓶颈。因此此处使用第二种方案：先将矩阵转置，直接使用转置矩阵乘向量：

 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

bool BiCG(
        Vector *x,
        const CsrMatrix &A,
        const CsrMatrix &At,
        const Vector &b,
        ValueType error,
        IndexType max_steps,
        const ValueType int_bound = 1e-15
) {
    const IndexType N = b.size;

    // r_0 = b - Ax_0
    Vector r(N);
    Residual(&r, A, b, *x);

    ValueType norm_r = Norm2(r);

    if (max_steps <= 0)
        max_steps = std::numeric_limits<IndexType>::max();

    // 初始解精度已经足够高
    if (norm_r < error)
        return false;

    // 取r*_0 = r_0，以满足r*_0与r_0不正交
    Vector r_star = r;

    // p_1 = r_0, p*_1 = r*_0
    Vector p = r;
    Vector p_star = r_star;

    // r_k^T * r*_k and r_{k-1}^T * r*_{k-1}
    ValueType rho0 = Dot(r, r_star);
    ValueType rho1 = 0;

    // A * p_k and A^T * p_k
    Vector q(N);
    Vector q_star(N);

    for (IndexType k = 1; k <= max_steps; ++k) {
        // q_k = A * p_k
        SpMv(&q, A, p);

        // xi_k = (r_{k-1}^T * r*_{k-1}) / (p*_k^T * q_k)
        ValueType xi = rho0 / Dot(p_star, q);

        // xi = 0时算法中断
        if (xi > -int_bound && xi < int_bound)
            return true;

        // x_k = x_{k-1} + xi_k * p_k
        Add(x, *x, p, xi);

        // r_k = r{k-1} - xi_k * q_k
        Add(&r, r, q, -xi);

        norm_r = Norm2(r);

        if (norm_r < error)
            break;

        // q*_k = A^T * p*_k
        SpMv(&q_star, At, p_star);  // A的转置乘向量

        // r*_k = r*_{k-1} - xi_k * A^T * p*_k
        Add(&r_star, r_star, q_star, -xi);

        rho1 = rho0;
        rho0 = Dot(r, r_star);

        // theta_k = (r_k, r*_k) / (r_{k-1}, r*_{k-1})
        ValueType theta = rho0 / rho1;

        // p_{k+1} = r_k + theta_k * p_k
        Add(&p, r, p, theta);

        // p*_{k+1} = r*_k + theta_k * p*_k
        Add(&p_star, r_star, p_star, theta);
    }

    return false;
}

A和转置通过参数传入，因此调用方法为：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

Host::CsrMatrix A;     // generate A
Host::Vector b;        // generate b

// 数据传输
SingleGPU::CsrMatrix d_A(A);
SingleGPU::Vector d_b(b);

// 计算转置
SingleGPU::CsrMatrix d_At = d_A.transpose();

SingleGPU::Vector d_x(N, 1.0);   // init x with N ones, and solve Ax=b

SingleGPU::BiCG(&d_x, d_A, d_At, d_b, 1e-6, 100);  // 残量小于1e-6或者迭代100次后退出

需要注意：计算出来的浮点数在判断等于时不能用==，只有判断它们的差是否足够接近于0。这是因为浮点数的计算有舍入误差，除非是直接赋值，计算出来的浮点数结果和真实结果一般都有误差，所以不能使用==来逐二进制位比较。具体细节可以参考IEEE 754浮点数标准[7,8]。

§参考文献

[1] Musser, David. (1997). Introspective Sorting and Selection Algorithms. Software Practice and Experience.

[2] http://feihu.me/blog/2014/sgi-std-sort/.

[3] https://zh.wikipedia.org/wiki/奇偶排序.

[4] https://en.wikipedia.org/wiki/Bitonic_sorter.

[5] https://en.wikipedia.org/wiki/Radix_sort.

[6] https://thrust.github.io/doc/modules.html.

[7] IEEE Standard for Binary Floating-Point Arithmetic.

[8] https://zh.wikipedia.org/wiki/IEEE_754.

§思考问题

尝试使用CUDA在CSR存储格式下实现稳定双共轭梯度法(BiCGStab)。

相关文章：

加载评论