CUDSS reordering discrepency

benjamin.donnot · July 7, 2026, 3:23pm

Hello,

When developping an application that leverages cudss to solve linear system on the GPU I found an issue.

The returned solution, with default reordering algorithm gives a garbage solution. By that I mean solving J.dx = F and then computing the max residual ||J.dx - F|| gives really huge number in some cases (1e36 for nested dissection for example).

I manage to scale down the problem to a tiny (302x302) matrix with 1842 nnz and a corresponding RHS. that I am allowed to share if needed.

The checking script is (I think) rather “simple” (though it’s cuda c++ + pybind11 for python manipulation) :

load the data
initialize a “custom cdss setting” (transfer data to gpu, perform analysis etc.)
solve the problem J.dx = F with scipy, no error
for all permutation algorithm I could fine: solve the problem with GPU (cudss), compute the max residual from the solution

(I can also share the script, it is part of package that is in the process of being properly open sourced)

And I obtained the following results:

reordering_alg      cuDSS max|dx|       residual        verdict
default             6.56959e+36         3.98087e+37     GARBAGE
btf_colamd          0.433558            1.07721e-16     OK
colamd              0.433558            1.07721e-16     OK
amd                 2.4111e+09          5.50653e+10     GARBAGE
nested_dissection   1.6424e+36          9.95217e+36     GARBAGE
none                0.433558            2.49677e-15     OK

I did not validated with other things, like CUDSS_CONFIG_MATCHING_ALG for example.

My real usecase involve uniform batch size, which is not supported by COLAMD or BTF_COLAMD at the moment. So my only fallback would be “None”.

I am using cudss 0.8.0 with cuda 12.4 and the drivers are Driver Version: 550.144.03.

Let me know if you need further information.

Best

Benjamin

kvoronin · July 7, 2026, 4:26pm

Hi @benjamin.donnot

Could you share the matrix and how you setup cudss (incl all if any no-default options and how you create matrices)? Ideally we would want to see a matrix in the mtx format + a short snippet of C/C++ code but we can also handle other data as well.

Also you might want to check for any errros in the log after you enable CUDSS_LOG_LEVEL=5 in the env (or use other logging APIs from cuDSS Logging Features — NVIDIA cuDSS ).

We have a guide for simple gotchas cuDSS Tips and tricks — NVIDIA cuDSS and accuracy issues cuDSS Tips and tricks — NVIDIA cuDSS .

Without any other data, I would think that maybe you used the wrong indexing (but BTF_COLAMD giving the right residual suggests otherwise). Also, suspicious is that default != nested_dissection as they currently should give exactly the same result.

If the input matrix has a valid structure and indexing, then likely it is the number of perturbed pivots which skyrockets and the default pivoting epsilon likely does not fit the scaling of your matrix very well. Then, as the guide above suggests, you might need to change the pivoting epsilon, or turn on matching with scaling.

Thanks,
Kirill

benjamin.donnot · July 8, 2026, 7:22am

bug_cudss.zip (26.4 KB)
Hello,

Thanks for your quick response. I attached a zip folder with inside:

RHS.npy the RHS vector saved using python / numpy
J.mtx (the sparse matrix)
cudss_raw_repro.py the script to compute the error

The code is now part of a larger code base that I can trim down for a minimal reproducible example if needed. The codebase is here: GitHub - Grid2op/gpusim2grid · GitHub but only a tiny portion of the code is required to reproduce the bug. I made a “solve_cudss_raw” that uses cudss to solve only one system.

Also, suspicious is that default != nested_dissection

I think the default is correct, but the output of cudss does not seem to be “reproducible” from one experiment to another, for example if I start it once:

> python cudss_raw_repro.py J.mtx RHS.npy
loaded J (302, 302) nnz=1842, F (302,)

scipy  spsolve(J, F): max|dx| = 0.433558  (reference; residual = 1.11022e-16)

reordering_alg      cuDSS max|dx|       residual        verdict
default             1.6424e+36          9.95217e+36     GARBAGE
btf_colamd          0.433558            1.07721e-16     OK
colamd              0.433558            1.07721e-16     OK
amd                 2.20911e+21         5.04523e+22     GARBAGE
nested_dissection   1.6424e+36          9.95217e+36     GARBAGE
none                0.433558            2.59053e-15     OK

Another same run:

reordering_alg      cuDSS max|dx|       residual        verdict
default             3.28479e+36         1.99043e+37     GARBAGE
btf_colamd          0.433558            1.07721e-16     OK
colamd              0.433558            1.07906e-16     OK
amd                 1.12604e+09         2.57173e+10     GARBAGE
nested_dissection   2.05602e+22         1.24585e+23     GARBAGE
none                0.433558            5.96786e-15     OK

As you can see, none of the “garbage” output gives the same result from one run to another.

And if I start the code with the logging feature, I have :

scipy  spsolve(J, F): max|dx| = 0.433558  (reference; residual = 1.11022e-16)

reordering_alg      cuDSS max|dx|       residual        verdict
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssSetStream] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigSet] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateCsr] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 3
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] CUDSS_CONFIG_REORDERING_ALG CUDSS_REORDERING_ALG_NESTED_DISSECTION requires = 97340 bytes (9.734e-05 GB) in host memory
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Using 1 threads on host for the reordering
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 4
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 1008
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDestroy] start
default             4.13562e+24         2.50599e+25     GARBAGE
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssSetStream] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigSet] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateCsr] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 3
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 4
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 1008
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDestroy] start
btf_colamd          0.433558            1.07721e-16     OK
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssSetStream] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigSet] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateCsr] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 3
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 4
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 1008
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDestroy] start
colamd              0.433558            1.07721e-16     OK
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssSetStream] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigSet] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateCsr] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 3
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 4
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 1008
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDestroy] start
amd                 1.05057e+09         2.39932e+10     GARBAGE
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssSetStream] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigSet] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateCsr] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 3
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] CUDSS_CONFIG_REORDERING_ALG CUDSS_REORDERING_ALG_NESTED_DISSECTION requires = 97340 bytes (9.734e-05 GB) in host memory
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Using 1 threads on host for the reordering
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 4
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 1008
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDestroy] start
nested_dissection   4.92719e+36         2.98565e+37     GARBAGE
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssSetStream] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataCreate] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigSet] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateCsr] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixCreateDn] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 3
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 4
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssExecute] start
[2026-07-08 09:05:51][CUDSS][487843][Info][cudssExecute] Run with phase 1008
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssMatrixDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDataDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssConfigDestroy] start
[2026-07-08 09:05:51][CUDSS][487843][Api][cudssDestroy] start
none                0.433558            2.59456e-15     OK

I am not an expert but this looks “rather good” to me.

you setup cudss (incl all if any no-default options and how you create matrices)?

The matrices are part of the Jacobian matrices when solving powerflow (power system analysis) using the “Newton Raphson” algorithm.

I did not run any non default parameters, as far as I know, when setting up cudss.

We have a guide for simple gotchas cuDSS Tips and tricks — NVIDIA cuDSS and accuracy issues cuDSS Tips and tricks — NVIDIA cuDSS .

Thank you very much, I will have a deeper look there.

Without any other data, I would think that maybe you used the wrong indexing

In the codebase, I made lots of tests to make sure everything was correct. There is no CI yet (looking for a CI plateform that would allow me to run cuda “for free” for “open source projects”, but that’s unrelated).

And for small powergrids I don’t have any issues (all tests pass). The issues are when I tried to run it on larger ones.

And when trying to debug, I ended up finding the issues was cudss giving wrong results. I then spent a bit of time to “focus” the issue on a smaller matrices (reducing the original powergrid as much as possible while the problem still remains) that I could share.

If the input matrix has a valid structure and indexing,

I think the structure is valid for different reasons:

it’s the same codebase that pass regular tests (same path in the code as working examples)
scipy is able to solve the problem well without any issues
cudss is able to solve the problem (if using BTF and COLAMD or COLAMD alone), but for the “big picture” I cannot use COLAMD (with or without BTF) when using “uniform batch mode” (present in the codebase but not required to reproduce the bug)

then likely it is the number of perturbed pivots which skyrockets and the default pivoting epsilon likely does not fit the scaling of your matrix very well. Then, as the guide above suggests, you might need to change the pivoting epsilon, or turn on matching with scaling.

Again thank you for the tips and for the links, I will try them and see how well they work.

Thank you and let me know if you want a minimal repo with only python + pybind11 + minimal cuda file repo to reproduce the bug.

Best regards

Benjamin

benjamin.donnot · July 8, 2026, 7:56am

bug_cudss_simple.zip (32.3 KB)
Hello again,

I made an even smaller reproducible example above. No need to pull anything from github. To run it you need: python + pybind11 + numpy + pip and a working c++17 compiler (and of course cuda and cudss :-) )

Everything is self contained and inside this repository (see readme if it can’t find python or pybind11 headers) :

mkdir build && cd build
cmake .. && make -j
PYTHONPATH=. python ../cudss_raw_repro.py ../J.mtx ../RHS.npy

Withtout any code from anywhere else.

Benjamin

kvoronin · July 8, 2026, 5:04pm

Hi @benjamin.donnot!

Thanks for the reproducer, I will run experiments locally and report the findings.

Best regards,
Kirill

kvoronin · July 9, 2026, 12:02am

Hi @benjamin.donnot!

Update: checked your matrix.

Findings:

You need to enable matching (for the considered matrix, any choice would work, but I recommend the default setting CUDSS_MATCHING_ALG_AUTO). Note: we might enable matching by default in one of the next releases.
Without matching, you get perturbed pivots due to having 0s on the diagonal (and the fact that reordering is symmetric and can’t move them away) which degrade the accuracy quite badly so that iterative refinement cannot recover it.
Using CUDSS_CONFIG_PIVOT_EPSILON_ALG = CUDSS_PIVOT_EPSILON_ALG_SCALED does not help either.

Best regards,
Kirill

benjamin.donnot · July 10, 2026, 7:01am

Hello,

Thanks a lot for your time and help.

Indeed matching seems to work for this matrix.

Best

Benjamin

Topic		Replies	Views
cuDSS is sometimes wrong where cuSparse and Umfpack succeed GPU-Accelerated Libraries cuda , cusparse , cudss	3	270	August 28, 2025
Cudss generates different solutions GPU-Accelerated Libraries cudss	3	281	July 17, 2024
Cusolver solve sparse Ax=b wrong GPU-Accelerated Libraries cusolver	7	313	August 19, 2024
cuSolverSp_LinearSolver produces different results with different order methods CUDA Programming and Performance	6	1594	September 21, 2017
Sparse matrix solve fails with cuda-sparse with cuDSS enabled GPU-Accelerated Libraries cudss	1	52	April 23, 2026
How to obtain only the re-ordered matrix from cuDSS? GPU-Accelerated Libraries cuda , cudss	1	100	August 12, 2025
What is the strategy for reorder in the cuDss? GPU-Accelerated Libraries cudss	2	485	July 3, 2024
Example using cusparse and cusolverSpDcsrlsvchol GPU-Accelerated Libraries cusolver , cusparse	22	639	May 13, 2025
Wrong result with cuDSS GPU-Accelerated Libraries cuda , cudss	3	195	June 10, 2025
cuDSS in a Fortran MPI program leading to segfault nvc, nvc++ and nvfortran cusparse , cudss	7	305	July 18, 2025

CUDSS reordering discrepency

Related topics