search

exanauts / CUDSS.jl

MIT License
18 stars 1 forks source link

Irreproducible results for the same linear system #22

Closed yuwenchen95 closed 5 months ago

yuwenchen95 commented 5 months ago

I found cudss can't reproduce the same result:

using CUDA, CUDA.CUSPARSE
using CUDSS 
using SparseArrays, LinearAlgebra
using Random
Random.seed!(666)

n = 1000
T = Float64
A_cpu = sprand(T, n, n, 0.05) + 0.01*I
x_cpu = zeros(T, n)
b_cpu = rand(T, n)

A_gpu = CuSparseMatrixCSR(A_cpu)
x_gpu = CuVector(x_cpu)
b_gpu = CuVector(b_cpu)

solver = CudssSolver(A_gpu, "G", 'F')

cudss("analysis", solver, x_gpu, b_gpu)
cudss("factorization", solver, x_gpu, b_gpu)
cudss("solve", solver, x_gpu, b_gpu)

x1 = deepcopy(x_gpu)
ldiv!(x_gpu,solver,b_gpu)
x2 = deepcopy(x_gpu)

println("error is ", norm(x1-x2))

It just want to solve the same linear system twice but the values of x_gpu are different, i.e. norm(x1-x2) is nonzero.

Is it a bug or the internal property of the method underlying cudss?

amontoison commented 5 months ago

@yuwenchen95 cuDSS is not deterministic. It's a very difficult property to ensure on GPU. Some libraries in CUDA, like CUSOLVER has an option for deterministic behaviour.