2012-12-31 64 views
1

我有一個小問題,我想將一個矩陣10 * 10轉換爲CSRCOO稀疏矩陣/格式。該矩陣是:以稀疏格式轉換矩陣A 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"的功能,但它沒有工作,一旦THA矩陣一個我想只是爲了在其他格式的轉換。 你能幫我嗎?

嗯,我也想用這個矩陣來解決系統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); 

使用結果八度應該類似於:

0.00000 
    0.32441 
    0.60970 
    0.82144 
    0.93411 
    0.93411 
    0.82144 
    0.60970 
    0.32441 
    0.00000 

但是也帶有負數,所以錯誤

回答

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; 

關於你試圖解決這個線性方程:會不會是八度上比你粘貼到你的問題的一個不同的矩陣操作?因爲與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 
+0

這個我今天將測試,謝謝。 – Ivan

+0

完美地工作。感謝您的快速幫助。我非常喜歡這個頁面。 – Ivan

+0

我有另外一個問題,但我不想打開一個新的話題。如果使用矩陣「A」,我想要求解方程系統Ax = b,使用** cusp :: krylov :: bicgstab(A,x,b,monitor,M); **與 * * b = [0.00000 0.34202 0.64279 0.86603 0.98481 0.98481 0.86603 0.64279 0.34202 0.00000 ] ** – Ivan