将矩阵 A 转换为稀疏格式 CSR、COO 等

Question

我有一个小问题，我想在 CSR 中转换一个矩阵 10*10或 COO稀疏矩阵/格式。矩阵是:

 1.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
-0.45  0.10 -0.45  0.00  0.00  0.00  0.00  0.00  0.00  0.00
 0.00 -0.45  0.10 -0.45  0.00  0.00  0.00  0.00  0.00  0.00
 0.00  0.00 -0.45  0.10 -0.45  0.00  0.00  0.00  0.00  0.00
 0.00  0.00  0.00 -0.45  0.10 -0.45  0.00  0.00  0.00  0.00
 0.00  0.00  0.00  0.00 -0.45  0.10 -0.45  0.00  0.00  0.00
 0.00  0.00  0.00  0.00  0.00 -0.45  0.10 -0.45  0.00  0.00
 0.00  0.00  0.00  0.00  0.00  0.00 -0.45  0.10 -0.45  0.00
 0.00  0.00  0.00  0.00  0.00  0.00  0.00 -0.45  0.10 -0.45
 0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  1.00

我正在使用 "CUSP"功能，但它没有工作，一旦矩阵 A 我只想转换为其他格式。 你能帮帮我吗？

好吧，我也想用这个矩阵来求解系统Ax=b，使用bicgstab :

b=
   0.00000
   0.34202
   0.64279
   0.86603
   0.98481
   0.98481
   0.86603
   0.64279
   0.34202
   0.00000

我的代码是:

int n = 10, r;

cusp::coo_matrix<int,float,cusp::device_memory> A(n, n, 3*n - 4);
cusp::array1d<float, cusp::device_memory> x(A.num_rows, 0);
cusp::array1d<float, cusp::device_memory> b(A.num_rows, 1);

   b[0]=0.00000;
   b[1]=0.34202;
   b[2]=0.64279;
   b[3]=0.86603;
   b[4]=0.98481;
   b[5]=0.98481;
   b[6]=0.86603;
   b[7]=0.64279;
   b[8]=0.34202;
   b[9]=0.00000;

i=0;
// row 0
A.row_indices[i] = 0.0;
A.column_indices[i] = 0.0;
A.values[i] = 1.00;

++i;
// rows 1 through n - 2
for (r = 1; r != n - 1; ++r) {
  A.row_indices[i] = r;
  A.column_indices[i] = r - 1;
  A.values[i] = -0.45;
  ++i;
  A.row_indices[i] = r;
  A.column_indices[i] = r;
  A.values[i] = 0.10;
  ++i;
  A.row_indices[i] = r;
  A.column_indices[i] = r + 1;
  A.values[i] = -0.45;
  ++i;
}
// row n - 1
A.row_indices[i] = n - 1;
A.column_indices[i] = n - 1;
A.values[i] = 1.00;
++i;

// set stopping criteria:
//  iteration_limit    = 100
//  relative_tolerance = 1e-3
    cusp::verbose_monitor<ValueType> monitor(b, 100, 1e-3);

// set preconditioner (identity)
cusp::identity_operator<ValueType, MemorySpace> M(A.num_rows, A.num_rows);

// solve the linear system A x = b
cusp::krylov::bicgstab(A, x, b, monitor, M);

cusp::print(x);

使用 Octave 的结果应该类似于:

   0.00000
   0.32441
   0.60970
   0.82144
   0.93411
   0.93411
   0.82144
   0.60970
   0.32441
   0.00000

但也有负数，所以错误。

Answer 1

对于 COO，您必须为每个条目设置三个数组元素:行和列索引以及值。您可以使用这样的代码为 COO 创建一个类似于您描述的矩阵:

int n = 10, i = 0, r;
cusp::csr_matrix<int,float,cusp::host_memory> A(n, n, 3*n - 4);
// row 0
A.row_indices[i] = 0;
A.column_indices[i] = 0;
A.values[i] = 1.00;
++i;
// rows 1 through n - 2
for (r = 1; r != n - 1; ++r) {
  A.row_indices[i] = r;
  A.column_indices[i] = r - 1;
  A.values[i] = -0.45;
  ++i;
  A.row_indices[i] = r;
  A.column_indices[i] = r;
  A.values[i] = 0.10;
  ++i;
  A.row_indices[i] = r;
  A.column_indices[i] = r + 1;
  A.values[i] = -0.45;
  ++i;
}
// row n - 1
A.row_indices[i] = n - 1;
A.column_indices[i] = n - 1;
A.values[i] = 1.00;
++i;

对于 CSR，您必须为每个条目指定一列和一个值，以及每一行的第一个条目的索引，包括最后一个行的最后一个索引。 CSR 的一段类似代码:

int n = 10, i = 0, r = 0;
cusp::csr_matrix<int,float,cusp::host_memory> A(n, n, 3*n - 4);
// row 0
A.row_offsets[r] = i;
A.column_indices[i] = 0;
A.values[i] = 1.00;
++i;
// rows 1 through n - 2
for (++r; r != n - 1; ++r) {
  A.row_offsets[r] = i;
  A.column_indices[i] = r - 1;
  A.values[i] = -0.45;
  ++i;
  A.column_indices[i] = r;
  A.values[i] = 0.10;
  ++i;
  A.column_indices[i] = r + 1;
  A.values[i] = -0.45;
  ++i;
}
// row n - 1
A.row_offsets[r] = i;
A.column_indices[i] = r;
A.values[i] = 1.00;
++i;
++r;
A.row_offsets[r] = i;

要从其他格式“转换”矩阵，您必须让我们知道您的原始数据以何种形式存储。从 cusp::array2d 转换而来应该通过简单地将数组传递给构造函数来工作。一般来说，像上面的代码一样首先以稀疏格式创建矩阵将提供更好的可扩展性。

另请注意，您的示例矩阵排列在对角带中，因此 cusp::dia_matrix会更适合，无论是在易于编码还是在更好的性能方面。要创建这样的三对角矩阵，您可以使用以下代码:

int n = 10, r = 0;
cusp::dia_matrix<int,float,cusp::host_memory> A(n, n, 3*n - 4, 3);
A.diagonal_offsets[0] = -1;
A.diagonal_offsets[1] = 0;
A.diagonal_offsets[2] = 1;
// row 0
A.values(r,0) = A.values(r,2) = 0.00;
A.values(r,1) = 1.00;
// rows 1 through n - 2
for (++r; r != n - 1; ++r) {
  A.values(r,0) = A.values(r,2) = -0.45;
  A.values(r,1) = 0.10;
}
// row n - 1
A.values(r,0) = A.values(r,2) = 0.00;
A.values(r,1) = 1.00;

关于您尝试求解的这个线性方程:Octave 是否在与您粘贴到问题中的矩阵不同的矩阵上运行？因为有 sage我在结果中也得到负数:

n = 10
d = dict()
d[(0,0)] = d[(n-1, n-1)] = 1
for r in range(1, n-1):
    d[(r, r-1)] = d[(r, r+1)] = -45/100
    d[(r,r)] = 1/10
A = matrix(RDF, n, n, d)
b = vector(RDF, [
   0.00000,
   0.34202,
   0.64279,
   0.86603,
   0.98481,
   0.98481,
   0.86603,
   0.64279,
   0.34202,
   0.00000,
   ])
for i in A.solve_right(b):
    print('{:+.5f}'.format(float(i)))

给出以下 vector x:

+0.00000
-0.45865
-0.86197
-1.16132
-1.32062
-1.32062
-1.16132
-0.86197
-0.45865
+0.00000