DirectSolver yields incorrect result for random and close to random matrices

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

Start with example01_cupy.py and reproduce the reported cases by varying n and the diagonal value. Compare DirectSolver results with the CuPy and NumPy outputs, including the norms and residuals shown in the issue. Done means identifying and correcting the cause of the incorrect results for nonsingular and singular matrices.