0

我正在尝试使用 CSparse 库（它是https://github.com/DrTimothyAldenDavis/SuiteSparse套件的一部分）对稀疏矩阵执行 LU 分解以求解稀疏线性矩阵方程。

为了验证功能，我在 MATLAB 中使用任意行、列、值三元组值定义了一个A大小为 的稀疏矩阵，条件是确保系统可以被分解/求解，以及一个大小为 的密集矩阵。(10, 10)cond(full(A)) < 100B(10)

使用下面代码中定义的值，A矩阵可以简单地在 MATLAB 中使用[L, U, P] = lu(A).

#include "cs.h"

int A_nonzero_elements = 19;
int A_indices_x_accessor[] = {1, 6, 5, 4, 9, 0, 2, 6, 8, 1,
                              2, 7, 4, 2, 3, 5, 7, 6, 7};
int A_indices_y_accessor[] = {0, 0, 1, 2, 2, 3, 4, 5, 5, 6,
                              6, 6, 7, 8, 8, 8, 8, 9, 9};
double A_values_accessor[] = {-0.563833742338016, -0.591908940985609, 0.559849215154234,
                              -0.253963929329736, -0.746974224826767, 0.148071275980455,
                              -0.730804912317641, 0.711443743939293,  -1.04374045545080,
                              -0.965352009768567, -0.565568775086541, 0.860497537960044,
                              1.03140136113829,   -0.941106875100149, -0.164924307450810,
                              0.205878365083051,  0.233891294833039,  -1.57903145117562,
                              0.385992304174899};
int m = 10;
int n = 10;
// Construct a sparse matrix A
cs *A = cs_spalloc(m, n, A_nonzero_elements, 1, 1);
for (int i = 0; i < A_nonzero_elements; i++) {
     assert(cs_entry(A, (csi)A_indices_x_accessor[i], (csi)A_indices_y_accessor[i], A_values_accessor[i]));
}
cs *A_compressed = cs_compress(A);
cs_spfree(A);
// Construct a dense matrix B
double *B = (double *)malloc(n * sizeof(double));
for (int i = 0; i < n; i++) {
     B[i] = 1.0;
}
// Solve AX=B using LU factorization
auto ok = cs_lusol(2, A_compressed, B, std::numeric_limits<double>::epsilon());
printf("%d\n", ok);

但是，此例程失败，因为cs_lu top = cs_spsolve (L, A, col, xi, x, pinv, 1)返回值10，因此循环for (p = top ; p < n ; p++)不会执行/迭代。

我尝试使用不同的order和tol值（cs_lusol分别是 的第一个和最后一个参数），但表现出相同的行为。

MATLAB 返回以下 L、U 和 P 矩阵：

L= (1,1)       1.0000
   (7,1)       0.9526
   (2,2)       1.0000
   (3,3)       1.0000
   (8,3)       0.3400
   (4,4)       1.0000
   (5,5)       1.0000
   (6,6)       1.0000
   (7,6)       0.6493
   (7,7)       1.0000
   (9,7)      -0.8914
   (8,8)       1.0000
   (9,9)       1.0000
  (10,9)      -0.7051
  (10,10)      1.0000



U = (1,1)      -0.5919
    (2,2)       0.5598
    (3,3)      -0.7470
    (4,4)       0.1481
    (5,5)      -0.7308
    (1,6)       0.7114
    (6,6)      -1.0437
    (5,7)      -0.5656
    (7,7)      -0.9654
    (8,8)       1.0314
    (2,9)       0.2059
    (5,9)      -0.9411
    (9,9)       0.2339
    (1,10)     -1.5790
    (7,10)      1.5041
    (9,10)      1.7268
    (10,10)      1.2176

P = (4,1)        1
    (7,2)        1
    (5,3)        1
    (10,4)       1
    (8,5)        1
    (2,6)        1
    (1,7)        1
    (9,8)        1
    (6,9)        1
    (3,10)       1

，并且能够毫无问题地解决系统。

是否有人能够识别上述代码的任何语义错误或其他问题？任何帮助将不胜感激。

4

0 回答 0