用SparseSolver求解双精度方程组时出现问题

综合讨论区
#1

大家好，我在使用SparseSolver求解双精度方程时出现问题（返回solution均为0且solver.info()==false）。

在下面这段示例代码中，A为可逆的对角矩阵。b所有元素为1。无法求解x。

import taichi as ti
import numpy as np



arch = ti.cpu # or ti.cuda
ti.init(arch=arch, default_fp=ti.f64)

n = 6



K = ti.linalg.SparseMatrixBuilder(n, n, dtype=ti.f64, max_num_triplets=100 )
b = ti.ndarray(ti.f64, shape=n)


@ti.kernel
def fill(A: ti.types.sparse_matrix_builder(), b: ti.types.ndarray(), interval: ti.i32):
    for i in range(n):
        A[i, i] += 2.0

        if i % interval == 0:
            b[i] += 1.0


fill(K, b, 1)

A = K.build()
print(">>>> Matrix A:")
print(A)
print(">>>> Vector b:")
print(b.to_numpy())
# outputs:
# >>>> Matrix A:
# [2, 0, 0, 0]
# [0, 2, 0, 0]
# [0, 0, 2, 0]
# [0, 0, 0, 2]
# >>>> Vector b:
# [1. 1. 1. 1.]
solver = ti.linalg.SparseSolver(solver_type="LLT")
solver.analyze_pattern(A)
solver.factorize(A)
x = solver.solve(b)
success = solver.info()
print(">>>> Solve sparse linear systems Ax = b with the solution x:")
print(x.to_numpy())
print(f">>>> Computation succeed: {success}")
# outputs:
# >>>> Solve sparse linear systems Ax = b with the solution x:
# [0.5 0.5 0.5 0.5]
# >>>> Computation was successful?: True
#2

在 @houkensjtu 的帮助下解决了，需要对于ti.linalg.SparseSolver也指定数据类型ti.f64

可以得到正确解的代码：

import taichi as ti
import numpy as np



arch = ti.cpu # or ti.cuda
ti.init(arch=arch, default_fp=ti.f64)

n = 6



K = ti.linalg.SparseMatrixBuilder(n, n, dtype=ti.f64, max_num_triplets=100 )
b = ti.ndarray(ti.f64, shape=n)



@ti.kernel
def fill(A: ti.types.sparse_matrix_builder(), b: ti.types.ndarray(), interval: ti.i32):
    for i in range(n):
        A[i, i] += 2.0

        if i % interval == 0:
            b[i] += 1.0





fill(K, b, 1)

A = K.build()
print(">>>> Matrix A:")
print(A)
print(">>>> Vector b:")
print(b.to_numpy())
# outputs:
# >>>> Matrix A:
# [2, 0, 0, 0]
# [0, 2, 0, 0]
# [0, 0, 2, 0]
# [0, 0, 0, 2]
# >>>> Vector b:
# [1. 1. 1. 1.]
solver = ti.linalg.SparseSolver(solver_type="LLT", dtype=ti.f64)
solver.analyze_pattern(A)
solver.factorize(A)
x = solver.solve(b)
success = solver.info()
print(">>>> Solve sparse linear systems Ax = b with the solution x:")
print(x.to_numpy())
print(f">>>> Computation succeed: {success}")
# outputs:
# >>>> Solve sparse linear systems Ax = b with the solution x:
# [0.5 0.5 0.5 0.5]
# >>>> Computation was successful?: True