Skip to content

Pure-Julia Sparse Cholesky #721

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Aug 12, 2025
Merged

Pure-Julia Sparse Cholesky #721

merged 3 commits into from
Aug 12, 2025

Conversation

samuelsonric
Copy link
Contributor

@samuelsonric samuelsonric commented Aug 12, 2025

This PR adds the Cholesky factorization algorithm implemented in CliqueTrees.jl. It is written in pure-Julia and accepts arbitrary numeric types. Somewhat slower than CHOLMOD. See this discourse thread.

Example

using BenchmarkTools, CliqueTrees, LinearSolve, MatrixMarket, SuiteSparseMatrixCollection
ssmc = ssmc_db(); name = "nasasrb";
A = mmread(joinpath(fetch_ssmc(ssmc[ssmc.name .== name, :]; format="MM")[1], "$(name).mtx"))
b = rand(size(A, 2))
problem = LinearProblem(A, b)

println("CholeskyFactorization:")
@btime solve(problem, CholeskyFactorization())

println("CliqueTreesFactorization:")
@btime solve(problem, CliqueTreesFactorization());

output:

CholeskyFactorization:
  227.126 ms (137 allocations: 327.23 MiB)
CliqueTreesFactorization:
  255.163 ms (103 allocations: 251.56 MiB)

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

Comment on lines 48 to 59
A = [
3 1 0 0 0 0 0 0
1 3 1 0 0 2 0 0
0 1 3 1 0 1 2 1
0 0 1 3 0 0 0 0
0 0 0 0 3 1 1 0
0 2 1 0 1 3 0 0
0 0 2 0 1 0 3 1
0 0 1 0 0 0 1 3
];

b = rand(8)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this not supposed to be a sparse matrix?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sparse, now. Converted internally in any case: I was trying to avoid using SparseArrays from some reason.


b = rand(8)
prob = LinearProblem(A, b)
sol = solve(prob) # in case cliquetrees is used as default
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
sol = solve(prob) # in case cliquetrees is used as default

It's not.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had copied that from the Sparsepak extension. Gone, now.

@ChrisRackauckas ChrisRackauckas merged commit 736a9cd into SciML:main Aug 12, 2025
113 of 118 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants