我有一个小问题,我想转换矩阵10*10在CSR或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,使用bicstab:
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
But也有负数,所以WRONG
对于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;
关于你试图解决的这个线性方程:可能是这个八度是在一个不同的矩阵上运行,而不是你粘贴到你的问题中?因为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)))
给出以下向量x:
+0.00000
-0.45865
-0.86197
-1.16132
-1.32062
-1.32062
-1.16132
-0.86197
-0.45865
+0.00000