DirectSolver yields incorrect result for random and close to random matrices

example01_cupy.py から始め、n と対角値を変化させて報告されたケースを再現します。DirectSolver の結果を CuPy と NumPy の出力と比較し、issue に示されているノルムと残差も含めます。非特異行列および特異行列で結果が正しくならない原因を特定して修正できれば完了です。

I was considering swapping cupy sparse solvers for nvmath Python/cuDSS DirectSolver, but I am unable to obtain correct results even for small problem instances.

Minimal example

In example01_cupy.py replace

n = 8
...
a += sp.diags([2.0] * n, format="csr", dtype="float64")

with

n = 100
...
# a += sp.diags([2.0] * n, format="csr", dtype="float64")

The same remains true, if I used small, but not tiny values for the diagonal like 0.01.

Example observation 1

||A||: 	40.65133430135082
||b||: 	14.142135623730951

det(A): 7.262125991606542e+28

02-10 17:35:07 userlogger   INFO     = SPECIFICATION PHASE =
02-10 17:35:07 userlogger   INFO     The LHS package is cupyx.
02-10 17:35:07 userlogger   INFO     The RHS package is cupy.
02-10 17:35:07 userlogger   INFO     The device_id=0, dtype = float64, index type = int32.
02-10 17:35:07 userlogger   INFO     The number of equations = 100.
02-10 17:35:07 userlogger   INFO     The operands' memory space is cuda, and the execution space is on device 0.
02-10 17:35:07 userlogger   INFO     The specified stream for the DirectSolver ctor is <cuda.core.experimental._stream.Stream object at 0x14b6709b87c0>.
02-10 17:35:07 userlogger   INFO     The library handle has been created: 94793135484384.
02-10 17:35:07 userlogger   INFO     The sparse direct solver operation has been created.
02-10 17:35:07 userlogger   INFO     Starting solver phase ANALYSIS...
02-10 17:35:07 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:35:07 userlogger   INFO     Starting solver phase FACTORIZATION...
02-10 17:35:07 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:35:07 userlogger   INFO     Starting solver phase SOLVE...
02-10 17:35:07 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:35:07 userlogger   INFO     The DirectSolver object's resources have been released.

 CuDSS:
||x||: 		29391823863.993862
||Ax - b||: 	92575181911.9934

CuPy:
||x||: 		10.647411227999223
||Ax - b||: 	4.289313285599112e-14

NumPy:
||x||: 		10.647411227999223
||Ax - b||: 	4.29904389535207e-14

Note that I computed the determinant to check whether A is invertible.

Example observation 2

Vice versa I made the observation that DirectSolver sometimes returns not nan even though the determinant is 0. Note that here b = cp.ones((n, 1), order="F") was only a vector

||A||: 	5.6922584694712395
||b||: 	10.0

det(A): 0.0

02-10 17:30:48 userlogger   INFO     = SPECIFICATION PHASE =
02-10 17:30:48 userlogger   INFO     The LHS package is cupyx.
02-10 17:30:48 userlogger   INFO     The RHS package is cupy.
02-10 17:30:48 userlogger   INFO     The device_id=0, dtype = float64, index type = int32.
02-10 17:30:48 userlogger   INFO     The number of equations = 100.
02-10 17:30:48 userlogger   INFO     The operands' memory space is cuda, and the execution space is on device 0.
02-10 17:30:48 userlogger   INFO     The specified stream for the DirectSolver ctor is <cuda.core.experimental._stream.Stream object at 0x14c6b6504580>.
02-10 17:30:48 userlogger   INFO     The library handle has been created: 94064601263200.
02-10 17:30:48 userlogger   INFO     The sparse direct solver operation has been created.
02-10 17:30:48 userlogger   INFO     Starting solver phase ANALYSIS...
02-10 17:30:48 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:30:48 userlogger   INFO     Starting solver phase FACTORIZATION...
02-10 17:30:48 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:30:48 userlogger   INFO     Starting solver phase SOLVE...
02-10 17:30:48 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:30:48 userlogger   INFO     The DirectSolver object's resources have been released.

 CuDSS:
||x||: 		1.5379474527090538e+102
||Ax - b||: 	1.5379474527090537e+89

CuPy:
||x||: 		nan
||Ax - b||: 	nan

NumPy:
||x||: 		nan
||Ax - b||: 	nan