From 16557b9f8842b6a1fa56396bb4481205816f26af Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 7 Jun 2023 14:17:23 +0200 Subject: [PATCH 1/3] workaround Cholesky factorization limitations on v1.6 It won't handle adjoint but that makes it error even when its not used from the default choice, and so this makes it so the defaults don't error. Past v1.7 it's not needed. --- src/factorization.jl | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 11e69b24c..823a5ee50 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -260,13 +260,7 @@ function init_cacheval(alg::CholeskyFactorization, A::Matrix{Float64}, b, u, Pl, PREALLOCATED_CHOLESKY end -function init_cacheval(alg::CholeskyFactorization, A::Diagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) - nothing -end - -function init_cacheval(alg::CholeskyFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, +function init_cacheval(alg::CholeskyFactorization, A::Union{Diagonal,AbstractSciMLOperator}, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) nothing @@ -281,7 +275,7 @@ end end function init_cacheval(alg::CholeskyFactorization, - A::Adjoint{<:Array}, b, u, Pl, Pr, + A::Adjoint{<:Number, <:Array}, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) nothing From c35900cb25b90ac740d23665ce259244c9262b9e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 8 Jun 2023 01:03:50 +0200 Subject: [PATCH 2/3] fix cholmod limitations --- src/factorization.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 823a5ee50..a5c41026c 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -838,10 +838,10 @@ function init_cacheval(alg::CHOLMODFactorization, nothing end -function init_cacheval(alg::CHOLMODFactorization, A::SparseMatrixCSC{Float64, Int}, b, u, +function init_cacheval(alg::CHOLMODFactorization, A::Union{SparseMatrixCSC{T, Int}, Symmetric{T, SparseMatrixCSC{T, Int}}}, b, u, Pl, Pr, maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + verbose::Bool, assumptions::OperatorAssumptions) where {T <: Union{Float32,Float64}} PREALLOCATED_CHOLMOD end From a1fedb81f36f3771038681e9736bf8b348dd8a46 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 8 Jun 2023 01:19:35 +0200 Subject: [PATCH 3/3] format --- docs/make.jl | 32 +- docs/pages.jl | 12 +- docs/src/advanced/developing.md | 2 +- docs/src/basics/FAQ.md | 2 +- ext/LinearSolveCUDAExt.jl | 2 +- ext/LinearSolveHYPREExt.jl | 70 ++-- ext/LinearSolveIterativeSolversExt.jl | 54 +-- ext/LinearSolveKrylovKitExt.jl | 8 +- ext/LinearSolvePardisoExt.jl | 26 +- src/LinearSolve.jl | 22 +- src/common.jl | 78 ++--- src/default.jl | 56 +-- src/deprecated.jl | 4 +- src/extension_algs.jl | 8 +- src/factorization.jl | 483 +++++++++++++------------- src/factorization_sparse.jl | 12 +- src/iterative_wrappers.jl | 24 +- src/simplelu.jl | 2 +- src/solve_function.jl | 2 +- test/basictests.jl | 52 +-- test/default_algs.jl | 2 +- test/hypretests.jl | 5 +- test/nonsquare.jl | 8 +- test/pardiso/pardiso.jl | 6 +- test/resolve.jl | 8 +- 25 files changed, 493 insertions(+), 487 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index dae4c0b62..677921189 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -9,21 +9,21 @@ DocMeta.setdocmeta!(LinearSolve, :DocTestSetup, :(using LinearSolve); recursive include("pages.jl") makedocs(sitename = "LinearSolve.jl", - authors = "Chris Rackauckas", - modules = [LinearSolve, LinearSolve.SciMLBase], - clean = true, doctest = false, linkcheck = true, - strict = [ - :doctest, - :linkcheck, - :parse_error, - :example_block, - # Other available options are - # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block - ], - format = Documenter.HTML(assets = ["assets/favicon.ico"], - canonical = "https://docs.sciml.ai/LinearSolve/stable/"), - pages = pages) + authors = "Chris Rackauckas", + modules = [LinearSolve, LinearSolve.SciMLBase], + clean = true, doctest = false, linkcheck = true, + strict = [ + :doctest, + :linkcheck, + :parse_error, + :example_block, + # Other available options are + # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block + ], + format = Documenter.HTML(assets = ["assets/favicon.ico"], + canonical = "https://docs.sciml.ai/LinearSolve/stable/"), + pages = pages) deploydocs(; - repo = "github.com/SciML/LinearSolve.jl", - push_preview = true) + repo = "github.com/SciML/LinearSolve.jl", + push_preview = true) diff --git a/docs/pages.jl b/docs/pages.jl index 762377418..d8f2e3cda 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,14 +2,14 @@ pages = ["index.md", "Tutorials" => Any["tutorials/linear.md" - "tutorials/caching_interface.md"], + "tutorials/caching_interface.md"], "Basics" => Any["basics/LinearProblem.md", - "basics/common_solver_opts.md", - "basics/OperatorAssumptions.md", - "basics/Preconditioners.md", - "basics/FAQ.md"], + "basics/common_solver_opts.md", + "basics/OperatorAssumptions.md", + "basics/Preconditioners.md", + "basics/FAQ.md"], "Solvers" => Any["solvers/solvers.md"], "Advanced" => Any["advanced/developing.md" - "advanced/custom.md"], + "advanced/custom.md"], "Release Notes" => "release_notes.md", ] diff --git a/docs/src/advanced/developing.md b/docs/src/advanced/developing.md index 62638bc48..f3861173b 100644 --- a/docs/src/advanced/developing.md +++ b/docs/src/advanced/developing.md @@ -18,7 +18,7 @@ basic machinery. A simplified version is: struct MyLUFactorization{P} <: SciMLBase.AbstractLinearAlgorithm end function init_cacheval(alg::MyLUFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, - verbose) + verbose) lu!(convert(AbstractMatrix, A)) end diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index 42bbf36c8..108f748b6 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -80,7 +80,7 @@ b = rand(n) weights = rand(n) realprec = lu(rand(n, n)) # some random preconditioner Pl = LinearSolve.ComposePreconditioner(LinearSolve.InvPreconditioner(Diagonal(weights)), - realprec) + realprec) Pr = Diagonal(weights) prob = LinearProblem(A, b) diff --git a/ext/LinearSolveCUDAExt.jl b/ext/LinearSolveCUDAExt.jl index 7dcd5173b..08acdaa28 100644 --- a/ext/LinearSolveCUDAExt.jl +++ b/ext/LinearSolveCUDAExt.jl @@ -4,7 +4,7 @@ using CUDA, LinearAlgebra, LinearSolve, SciMLBase using SciMLBase: AbstractSciMLOperator function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::CudaOffloadFactorization; - kwargs...) + kwargs...) if cache.isfresh fact = LinearSolve.do_factorization(alg, CUDA.CuArray(cache.A), cache.b, cache.u) cache = LinearSolve.set_cacheval(cache, fact) diff --git a/ext/LinearSolveHYPREExt.jl b/ext/LinearSolveHYPREExt.jl index 5da6b4083..00f969fe3 100644 --- a/ext/LinearSolveHYPREExt.jl +++ b/ext/LinearSolveHYPREExt.jl @@ -4,8 +4,8 @@ using LinearAlgebra using HYPRE.LibHYPRE: HYPRE_Complex using HYPRE: HYPRE, HYPREMatrix, HYPRESolver, HYPREVector using LinearSolve: HYPREAlgorithm, LinearCache, LinearProblem, LinearSolve, - OperatorAssumptions, default_tol, init_cacheval, __issquare, - __conditioning + OperatorAssumptions, default_tol, init_cacheval, __issquare, + __conditioning using SciMLBase: LinearProblem, SciMLBase using UnPack: @unpack using Setfield: @set! @@ -21,8 +21,8 @@ mutable struct HYPRECache end function LinearSolve.init_cacheval(alg::HYPREAlgorithm, A, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) return HYPRECache(nothing, nothing, nothing, nothing, true, true, true) end @@ -54,21 +54,21 @@ end # fill!(similar(b, size(A, 2)), false) since HYPREArrays are not AbstractArrays. function SciMLBase.init(prob::LinearProblem, alg::HYPREAlgorithm, - args...; - alias_A = false, alias_b = false, - # TODO: Implement eltype for HYPREMatrix in HYPRE.jl? Looks useful - # even if it is not AbstractArray. - abstol = default_tol(prob.A isa HYPREMatrix ? HYPRE_Complex : - eltype(prob.A)), - reltol = default_tol(prob.A isa HYPREMatrix ? HYPRE_Complex : - eltype(prob.A)), - # TODO: Implement length() for HYPREVector in HYPRE.jl? - maxiters::Int = prob.b isa HYPREVector ? 1000 : length(prob.b), - verbose::Bool = false, - Pl = LinearAlgebra.I, - Pr = LinearAlgebra.I, - assumptions = OperatorAssumptions(), - kwargs...) + args...; + alias_A = false, alias_b = false, + # TODO: Implement eltype for HYPREMatrix in HYPRE.jl? Looks useful + # even if it is not AbstractArray. + abstol = default_tol(prob.A isa HYPREMatrix ? HYPRE_Complex : + eltype(prob.A)), + reltol = default_tol(prob.A isa HYPREMatrix ? HYPRE_Complex : + eltype(prob.A)), + # TODO: Implement length() for HYPREVector in HYPRE.jl? + maxiters::Int = prob.b isa HYPREVector ? 1000 : length(prob.b), + verbose::Bool = false, + Pl = LinearAlgebra.I, + Pr = LinearAlgebra.I, + assumptions = OperatorAssumptions(), + kwargs...) @unpack A, b, u0, p = prob A = A isa HYPREMatrix ? A : HYPREMatrix(A) @@ -82,23 +82,23 @@ function SciMLBase.init(prob::LinearProblem, alg::HYPREAlgorithm, # Initialize internal alg cache cacheval = init_cacheval(alg, A, b, u0, Pl, Pr, maxiters, abstol, reltol, verbose, - assumptions) + assumptions) Tc = typeof(cacheval) isfresh = true cache = LinearCache{ - typeof(A), typeof(b), typeof(u0), typeof(p), typeof(alg), Tc, - typeof(Pl), typeof(Pr), typeof(reltol), - typeof(__issquare(assumptions)) - }(A, b, u0, p, alg, cacheval, isfresh, Pl, Pr, abstol, reltol, - maxiters, - verbose, assumptions) + typeof(A), typeof(b), typeof(u0), typeof(p), typeof(alg), Tc, + typeof(Pl), typeof(Pr), typeof(reltol), + typeof(__issquare(assumptions)), + }(A, b, u0, p, alg, cacheval, isfresh, Pl, Pr, abstol, reltol, + maxiters, + verbose, assumptions) return cache end # Solvers whose constructor requires passing the MPI communicator const COMM_SOLVERS = Union{HYPRE.BiCGSTAB, HYPRE.FlexGMRES, HYPRE.GMRES, HYPRE.ParaSails, - HYPRE.PCG} + HYPRE.PCG} create_solver(::Type{S}, comm) where {S <: COMM_SOLVERS} = S(comm) # Solvers whose constructor should not be passed the MPI communicator @@ -120,10 +120,10 @@ function create_solver(alg::HYPREAlgorithm, cache::LinearCache) # Construct solver options solver_options = (; - AbsoluteTol = cache.abstol, - MaxIter = cache.maxiters, - PrintLevel = Int(cache.verbose), - Tol = cache.reltol) + AbsoluteTol = cache.abstol, + MaxIter = cache.maxiters, + PrintLevel = Int(cache.verbose), + Tol = cache.reltol) # Preconditioner (uses Pl even though it might not be a *left* preconditioner just *a* # preconditioner) @@ -211,16 +211,16 @@ function SciMLBase.solve!(cache::LinearCache, alg::HYPREAlgorithm, args...; kwar stats = nothing ret = SciMLBase.LinearSolution{T, N, typeof(cache.u), typeof(resid), typeof(alg), - typeof(cache), typeof(stats)}(cache.u, resid, alg, retc, - iters, cache, stats) + typeof(cache), typeof(stats)}(cache.u, resid, alg, retc, + iters, cache, stats) return ret end # HYPREArrays are not AbstractArrays so perform some type-piracy function SciMLBase.LinearProblem(A::HYPREMatrix, b::HYPREVector, - p = SciMLBase.NullParameters(); - u0::Union{HYPREVector, Nothing} = nothing, kwargs...) + p = SciMLBase.NullParameters(); + u0::Union{HYPREVector, Nothing} = nothing, kwargs...) return LinearProblem{true}(A, b, p; u0 = u0, kwargs) end diff --git a/ext/LinearSolveIterativeSolversExt.jl b/ext/LinearSolveIterativeSolversExt.jl index 1cff07958..b111c8f77 100644 --- a/ext/LinearSolveIterativeSolversExt.jl +++ b/ext/LinearSolveIterativeSolversExt.jl @@ -11,31 +11,31 @@ else end function LinearSolve.IterativeSolversJL(args...; - generate_iterator = IterativeSolvers.gmres_iterable!, - gmres_restart = 0, kwargs...) + generate_iterator = IterativeSolvers.gmres_iterable!, + gmres_restart = 0, kwargs...) return IterativeSolversJL(generate_iterator, gmres_restart, - args, kwargs) + args, kwargs) end function LinearSolve.IterativeSolversJL_CG(args...; kwargs...) IterativeSolversJL(args...; - generate_iterator = IterativeSolvers.cg_iterator!, - kwargs...) + generate_iterator = IterativeSolvers.cg_iterator!, + kwargs...) end function LinearSolve.IterativeSolversJL_GMRES(args...; kwargs...) IterativeSolversJL(args...; - generate_iterator = IterativeSolvers.gmres_iterable!, - kwargs...) + generate_iterator = IterativeSolvers.gmres_iterable!, + kwargs...) end function LinearSolve.IterativeSolversJL_BICGSTAB(args...; kwargs...) IterativeSolversJL(args...; - generate_iterator = IterativeSolvers.bicgstabl_iterator!, - kwargs...) + generate_iterator = IterativeSolvers.bicgstabl_iterator!, + kwargs...) end function LinearSolve.IterativeSolversJL_MINRES(args...; kwargs...) IterativeSolversJL(args...; - generate_iterator = IterativeSolvers.minres_iterable!, - kwargs...) + generate_iterator = IterativeSolvers.minres_iterable!, + kwargs...) end LinearSolve._isidentity_struct(::IterativeSolvers.Identity) = true @@ -43,33 +43,33 @@ LinearSolve.default_alias_A(::IterativeSolversJL, ::Any, ::Any) = true LinearSolve.default_alias_b(::IterativeSolversJL, ::Any, ::Any) = true function LinearSolve.init_cacheval(alg::IterativeSolversJL, A, b, u, Pl, Pr, maxiters::Int, - abstol, - reltol, - verbose::Bool, assumptions::OperatorAssumptions) + abstol, + reltol, + verbose::Bool, assumptions::OperatorAssumptions) restart = (alg.gmres_restart == 0) ? min(20, size(A, 1)) : alg.gmres_restart kwargs = (abstol = abstol, reltol = reltol, maxiter = maxiters, - alg.kwargs...) + alg.kwargs...) iterable = if alg.generate_iterator === IterativeSolvers.cg_iterator! !LinearSolve._isidentity_struct(Pr) && @warn "$(alg.generate_iterator) doesn't support right preconditioning" alg.generate_iterator(u, A, b, Pl; - kwargs...) + kwargs...) elseif alg.generate_iterator === IterativeSolvers.gmres_iterable! alg.generate_iterator(u, A, b; Pl = Pl, Pr = Pr, restart = restart, - kwargs...) + kwargs...) elseif alg.generate_iterator === IterativeSolvers.bicgstabl_iterator! !!LinearSolve._isidentity_struct(Pr) && @warn "$(alg.generate_iterator) doesn't support right preconditioning" alg.generate_iterator(u, A, b, alg.args...; Pl = Pl, - abstol = abstol, reltol = reltol, - max_mv_products = maxiters * 2, - alg.kwargs...) + abstol = abstol, reltol = reltol, + max_mv_products = maxiters * 2, + alg.kwargs...) else # minres, qmr alg.generate_iterator(u, A, b, alg.args...; - abstol = abstol, reltol = reltol, maxiter = maxiters, - alg.kwargs...) + abstol = abstol, reltol = reltol, maxiter = maxiters, + alg.kwargs...) end return iterable end @@ -77,10 +77,10 @@ end function SciMLBase.solve!(cache::LinearCache, alg::IterativeSolversJL; kwargs...) if cache.isfresh || !(typeof(alg) <: IterativeSolvers.GMRESIterable) solver = LinearSolve.init_cacheval(alg, cache.A, cache.b, cache.u, cache.Pl, - cache.Pr, - cache.maxiters, cache.abstol, cache.reltol, - cache.verbose, - cache.assumptions) + cache.Pr, + cache.maxiters, cache.abstol, cache.reltol, + cache.verbose, + cache.assumptions) cache.cacheval = solver cache.isfresh = false end @@ -111,7 +111,7 @@ function purge_history!(iter::IterativeSolvers.GMRESIterable, x, b) iter.b = b iter.residual.current = IterativeSolvers.init!(iter.arnoldi, iter.x, iter.b, iter.Pl, - iter.Ax, initially_zero = true) + iter.Ax, initially_zero = true) IterativeSolvers.init_residual!(iter.residual, iter.residual.current) iter.β = iter.residual.current nothing diff --git a/ext/LinearSolveKrylovKitExt.jl b/ext/LinearSolveKrylovKitExt.jl index 7603ba0e1..71aa5dddf 100644 --- a/ext/LinearSolveKrylovKitExt.jl +++ b/ext/LinearSolveKrylovKitExt.jl @@ -4,8 +4,8 @@ using LinearSolve, KrylovKit, LinearAlgebra using LinearSolve: LinearCache function LinearSolve.KrylovKitJL(args...; - KrylovAlg = KrylovKit.GMRES, gmres_restart = 0, - kwargs...) + KrylovAlg = KrylovKit.GMRES, gmres_restart = 0, + kwargs...) return KrylovKitJL(KrylovAlg, gmres_restart, args, kwargs) end @@ -28,7 +28,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovKitJL; kwargs...) krylovdim = (alg.gmres_restart == 0) ? min(20, size(cache.A, 1)) : alg.gmres_restart kwargs = (atol = atol, rtol = rtol, maxiter = maxiter, verbosity = verbosity, - krylovdim = krylovdim, alg.kwargs...) + krylovdim = krylovdim, alg.kwargs...) x, info = KrylovKit.linsolve(cache.A, cache.b, cache.u; kwargs...) @@ -37,7 +37,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovKitJL; kwargs...) retcode = info.converged == 1 ? ReturnCode.Default : ReturnCode.ConvergenceFailure iters = info.numiter return SciMLBase.build_linear_solution(alg, cache.u, resid, cache; retcode = retcode, - iters = iters) + iters = iters) end end diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index 8f7e38a4c..afc486a8f 100644 --- a/ext/LinearSolvePardisoExt.jl +++ b/ext/LinearSolvePardisoExt.jl @@ -12,16 +12,16 @@ LinearSolve.needs_concrete_A(alg::PardisoJL) = true # TODO schur complement functionality function LinearSolve.init_cacheval(alg::PardisoJL, - A, - b, - u, - Pl, - Pr, - maxiters::Int, - abstol, - reltol, - verbose::Bool, - assumptions::LinearSolve.OperatorAssumptions) + A, + b, + u, + Pl, + Pr, + maxiters::Int, + abstol, + reltol, + verbose::Bool, + assumptions::LinearSolve.OperatorAssumptions) @unpack nprocs, solver_type, matrix_type, iparm, dparm = alg A = convert(AbstractMatrix, A) @@ -93,9 +93,9 @@ function LinearSolve.init_cacheval(alg::PardisoJL, end Pardiso.pardiso(solver, - u, - SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), - b) + u, + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), + b) return solver end diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 2bc9a28ba..ee4a10acb 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -91,7 +91,7 @@ include("extension_algs.jl") include("deprecated.jl") @generated function SciMLBase.solve!(cache::LinearCache, alg::AbstractFactorization; - kwargs...) + kwargs...) quote if cache.isfresh fact = do_factorization(alg, cache.A, cache.b, cache.u) @@ -99,7 +99,7 @@ include("deprecated.jl") cache.isfresh = false end y = _ldiv!(cache.u, @get_cacheval(cache, $(Meta.quot(defaultalg_symbol(alg)))), - cache.b) + cache.b) #= retcode = if LinearAlgebra.issuccess(fact) @@ -162,19 +162,19 @@ PrecompileTools.@compile_workload begin end export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization, - GenericLUFactorization, SimpleLUFactorization, RFLUFactorization, - NormalCholeskyFactorization, NormalBunchKaufmanFactorization, - UMFPACKFactorization, KLUFactorization, FastLUFactorization, FastQRFactorization, - SparspakFactorization, DiagonalFactorization, CholeskyFactorization, - BunchKaufmanFactorization, CHOLMODFactorization, LDLtFactorization + GenericLUFactorization, SimpleLUFactorization, RFLUFactorization, + NormalCholeskyFactorization, NormalBunchKaufmanFactorization, + UMFPACKFactorization, KLUFactorization, FastLUFactorization, FastQRFactorization, + SparspakFactorization, DiagonalFactorization, CholeskyFactorization, + BunchKaufmanFactorization, CHOLMODFactorization, LDLtFactorization export LinearSolveFunction, DirectLdiv! export KrylovJL, KrylovJL_CG, KrylovJL_MINRES, KrylovJL_GMRES, - KrylovJL_BICGSTAB, KrylovJL_LSMR, KrylovJL_CRAIGMR, - IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES, - IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES, - KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES + KrylovJL_BICGSTAB, KrylovJL_LSMR, KrylovJL_CRAIGMR, + IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES, + IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES, + KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES export HYPREAlgorithm export CudaOffloadFactorization diff --git a/src/common.jl b/src/common.jl index 262fa2011..953a82d3c 100644 --- a/src/common.jl +++ b/src/common.jl @@ -59,7 +59,7 @@ struct OperatorAssumptions{T} end function OperatorAssumptions(issquare = nothing; - condition::OperatorCondition.T = OperatorCondition.IllConditioned) + condition::OperatorCondition.T = OperatorCondition.IllConditioned) OperatorAssumptions{typeof(issquare)}(issquare, condition) end __issquare(assump::OperatorAssumptions) = assump.issq @@ -112,17 +112,17 @@ default_alias_A(::AbstractKrylovSubspaceMethod, ::Any, ::Any) = true default_alias_b(::AbstractKrylovSubspaceMethod, ::Any, ::Any) = true function SciMLBase.init(prob::LinearProblem, alg::SciMLLinearSolveAlgorithm, - args...; - alias_A = default_alias_A(alg, prob.A, prob.b), - alias_b = default_alias_b(alg, prob.A, prob.b), - abstol = default_tol(eltype(prob.A)), - reltol = default_tol(eltype(prob.A)), - maxiters::Int = length(prob.b), - verbose::Bool = false, - Pl = IdentityOperator(size(prob.A)[1]), - Pr = IdentityOperator(size(prob.A)[2]), - assumptions = OperatorAssumptions(issquare(prob.A)), - kwargs...) + args...; + alias_A = default_alias_A(alg, prob.A, prob.b), + alias_b = default_alias_b(alg, prob.A, prob.b), + abstol = default_tol(eltype(prob.A)), + reltol = default_tol(eltype(prob.A)), + maxiters::Int = length(prob.b), + verbose::Bool = false, + Pl = IdentityOperator(size(prob.A)[1]), + Pr = IdentityOperator(size(prob.A)[2]), + assumptions = OperatorAssumptions(issquare(prob.A)), + kwargs...) @unpack A, b, u0, p = prob A = if alias_A @@ -151,35 +151,35 @@ function SciMLBase.init(prob::LinearProblem, alg::SciMLLinearSolveAlgorithm, end cacheval = init_cacheval(alg, A, b, u0, Pl, Pr, maxiters, abstol, reltol, verbose, - assumptions) + assumptions) isfresh = true Tc = typeof(cacheval) cache = LinearCache{ - typeof(A), - typeof(b), - typeof(u0), - typeof(p), - typeof(alg), - Tc, - typeof(Pl), - typeof(Pr), - typeof(reltol), - typeof(assumptions.issq) - }(A, - b, - u0, - p, - alg, - cacheval, - isfresh, - Pl, - Pr, - abstol, - reltol, - maxiters, - verbose, - assumptions) + typeof(A), + typeof(b), + typeof(u0), + typeof(p), + typeof(alg), + Tc, + typeof(Pl), + typeof(Pr), + typeof(reltol), + typeof(assumptions.issq), + }(A, + b, + u0, + p, + alg, + cacheval, + isfresh, + Pl, + Pr, + abstol, + reltol, + maxiters, + verbose, + assumptions) return cache end @@ -188,8 +188,8 @@ function SciMLBase.solve(prob::LinearProblem, args...; kwargs...) end function SciMLBase.solve(prob::LinearProblem, - alg::Union{SciMLLinearSolveAlgorithm, Nothing}, - args...; kwargs...) + alg::Union{SciMLLinearSolveAlgorithm, Nothing}, + args...; kwargs...) solve!(init(prob, alg, args...; kwargs...)) end diff --git a/src/default.jl b/src/default.jl index 257545766..3eb67627a 100644 --- a/src/default.jl +++ b/src/default.jl @@ -1,6 +1,6 @@ needs_concrete_A(alg::DefaultLinearSolver) = true mutable struct DefaultLinearSolverInit{T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, - T13, T14, T15, T16} + T13, T14, T15, T16} LUFactorization::T1 QRFactorization::T2 DiagonalFactorization::T3 @@ -24,7 +24,7 @@ end defaultalg(A, b) = defaultalg(A, b, OperatorAssumptions(true)) function defaultalg(A::Union{DiffEqArrayOperator, MatrixOperator}, b, - assump::OperatorAssumptions) + assump::OperatorAssumptions) defaultalg(A.A, b, assump) end @@ -67,7 +67,7 @@ function defaultalg(A::Symmetric{<:Number, <:SparseMatrixCSC}, b, ::OperatorAssu end function defaultalg(A::AbstractSparseMatrixCSC{Tv, Ti}, b, - assump::OperatorAssumptions) where {Tv, Ti} + assump::OperatorAssumptions) where {Tv, Ti} if assump.issq DefaultLinearSolver(DefaultAlgorithmChoice.SparspakFactorization) else @@ -77,7 +77,7 @@ end @static if INCLUDE_SPARSE function defaultalg(A::AbstractSparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b, - assump::OperatorAssumptions) where {Ti} + assump::OperatorAssumptions) where {Ti} if assump.issq if length(b) <= 10_000 DefaultLinearSolver(DefaultAlgorithmChoice.KLUFactorization) @@ -117,7 +117,7 @@ end # Ambiguity handling function defaultalg(A::GPUArraysCore.AbstractGPUArray, b::GPUArraysCore.AbstractGPUArray, - assump::OperatorAssumptions) + assump::OperatorAssumptions) if assump.condition === OperatorCondition.IllConditioned || !assump.issq DefaultLinearSolver(DefaultAlgorithmChoice.QRFactorization) else @@ -130,7 +130,7 @@ function defaultalg(A::GPUArraysCore.AbstractGPUArray, b::GPUArraysCore.Abstract end function defaultalg(A::SciMLBase.AbstractSciMLOperator, b, - assump::OperatorAssumptions) + assump::OperatorAssumptions) if has_ldiv!(A) return DefaultLinearSolver(DefaultAlgorithmChoice.DirectLdiv!) elseif !assump.issq @@ -240,33 +240,33 @@ end ## Catch high level interface function SciMLBase.init(prob::LinearProblem, alg::Nothing, - args...; - assumptions = OperatorAssumptions(issquare(prob.A)), - kwargs...) + args...; + assumptions = OperatorAssumptions(issquare(prob.A)), + kwargs...) alg = defaultalg(prob.A, prob.b, assumptions) SciMLBase.init(prob, alg, args...; assumptions, kwargs...) end function SciMLBase.solve!(cache::LinearCache, alg::Nothing, - args...; assump::OperatorAssumptions = OperatorAssumptions(), - kwargs...) + args...; assump::OperatorAssumptions = OperatorAssumptions(), + kwargs...) @unpack A, b = cache SciMLBase.solve!(cache, defaultalg(A, b, assump), args...; kwargs...) end function init_cacheval(alg::Nothing, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, - verbose::Bool, assump::OperatorAssumptions) + verbose::Bool, assump::OperatorAssumptions) init_cacheval(defaultalg(A, b, assump), A, b, u, Pl, Pr, maxiters, abstol, reltol, - verbose, - assump) + verbose, + assump) end """ cache.cacheval = NamedTuple(LUFactorization = cache of LUFactorization, ...) """ @generated function init_cacheval(alg::DefaultLinearSolver, A, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, - verbose::Bool, assump::OperatorAssumptions) + abstol, reltol, + verbose::Bool, assump::OperatorAssumptions) caches = map(first.(EnumX.symbol_map(DefaultAlgorithmChoice.T))) do alg if alg === :KrylovJL_GMRES quote @@ -274,17 +274,17 @@ cache.cacheval = NamedTuple(LUFactorization = cache of LUFactorization, ...) nothing else init_cacheval($(algchoice_to_alg(alg)), A, b, u, Pl, Pr, maxiters, - abstol, reltol, - verbose, - assump) + abstol, reltol, + verbose, + assump) end end else quote init_cacheval($(algchoice_to_alg(alg)), A, b, u, Pl, Pr, maxiters, abstol, - reltol, - verbose, - assump) + reltol, + verbose, + assump) end end end @@ -304,20 +304,20 @@ else end """ @generated function SciMLBase.solve!(cache::LinearCache, alg::DefaultLinearSolver, - args...; - assump::OperatorAssumptions = OperatorAssumptions(), - kwargs...) + args...; + assump::OperatorAssumptions = OperatorAssumptions(), + kwargs...) ex = :() for alg in first.(EnumX.symbol_map(DefaultAlgorithmChoice.T)) newex = quote sol = SciMLBase.solve!(cache, $(algchoice_to_alg(alg)), args...; kwargs...) SciMLBase.build_linear_solution(alg, sol.u, sol.resid, sol.cache; - retcode = sol.retcode, - iters = sol.iters, stats = sol.stats) + retcode = sol.retcode, + iters = sol.iters, stats = sol.stats) end ex = if ex == :() Expr(:elseif, :(Symbol(alg.alg) === $(Meta.quot(alg))), newex, - :(error("Algorithm Choice not Allowed"))) + :(error("Algorithm Choice not Allowed"))) else Expr(:elseif, :(Symbol(alg.alg) === $(Meta.quot(alg))), newex, ex) end diff --git a/src/deprecated.jl b/src/deprecated.jl index a1ac1b92b..d11d9a29e 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -81,8 +81,8 @@ end @static if VERSION >= v"1.7" @deprecate SciMLBase.solve(cache::LinearCache, args...; kwargs...) SciMLBase.solve!(cache::LinearCache, - args...; - kwargs...) false + args...; + kwargs...) false else function SciMLBase.solve(cache::LinearCache, args...; kwargs...) @warn "SciMLBase.solve(cache::LinearCache, args...; kwargs...) is deprecated for SciMLBase.solve!(cache::LinearCache, args...; kwargs...) to follow the CommonSolve.jl interface." diff --git a/src/extension_algs.jl b/src/extension_algs.jl index a8ef47b2d..43d5a25b8 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -163,10 +163,10 @@ MKLPardisoIterate(; kwargs...) = PardisoJL(; solver_type = 1, kwargs...) dparm::Union{Vector{Tuple{Int, Int}}, Nothing} function PardisoJL(; nprocs::Union{Int, Nothing} = nothing, - solver_type = nothing, - matrix_type = nothing, - iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, - dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) + solver_type = nothing, + matrix_type = nothing, + iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, + dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) ext = Base.get_extension(@__MODULE__, :LinearSolvePardisoExt) if ext === nothing error("PardisoJL requires that Pardiso is loaded, i.e. `using Pardiso`") diff --git a/src/factorization.jl b/src/factorization.jl index a5c41026c..ffbc9dcca 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -19,7 +19,7 @@ end #RF Bad fallback: will fail if `A` is just a stand-in # This should instead just create the factorization type. function init_cacheval(alg::AbstractFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, - reltol, verbose::Bool, assumptions::OperatorAssumptions) + reltol, verbose::Bool, assumptions::OperatorAssumptions) do_factorization(alg, convert(AbstractMatrix, A), b, u) end @@ -98,32 +98,32 @@ function do_factorization(alg::GenericLUFactorization, A, b, u) end function init_cacheval(alg::Union{LUFactorization, GenericLUFactorization}, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.lu_instance(convert(AbstractMatrix, A)) end const PREALLOCATED_LU = ArrayInterface.lu_instance(rand(1, 1)) function init_cacheval(alg::Union{LUFactorization, GenericLUFactorization}, - A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Matrix{Float64}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) PREALLOCATED_LU end function init_cacheval(alg::Union{LUFactorization, GenericLUFactorization}, - A::AbstractSciMLOperator, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end @static if VERSION < v"1.7-" function init_cacheval(alg::Union{LUFactorization, GenericLUFactorization}, - A::Union{Diagonal, SymTridiagonal}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{Diagonal, SymTridiagonal}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end end @@ -168,30 +168,30 @@ function do_factorization(alg::QRFactorization, A, b, u) end function init_cacheval(alg::QRFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.qr_instance(convert(AbstractMatrix, A)) end const PREALLOCATED_QR = ArrayInterface.qr_instance(rand(1, 1)) function init_cacheval(alg::QRFactorization, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) PREALLOCATED_QR end function init_cacheval(alg::QRFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end @static if VERSION < v"1.7-" function init_cacheval(alg::QRFactorization, - A::Union{Diagonal, SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{Diagonal, SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end end @@ -241,8 +241,8 @@ function do_factorization(alg::CholeskyFactorization, A, b, u) end function init_cacheval(alg::CholeskyFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.cholesky_instance(convert(AbstractMatrix, A), alg.pivot) end @@ -255,29 +255,30 @@ end const PREALLOCATED_CHOLESKY = ArrayInterface.cholesky_instance(rand(1, 1), cholpivot) function init_cacheval(alg::CholeskyFactorization, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) PREALLOCATED_CHOLESKY end -function init_cacheval(alg::CholeskyFactorization, A::Union{Diagonal,AbstractSciMLOperator}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) +function init_cacheval(alg::CholeskyFactorization, + A::Union{Diagonal, AbstractSciMLOperator}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end @static if VERSION < v"1.7beta" function init_cacheval(alg::CholeskyFactorization, - A::Union{SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end function init_cacheval(alg::CholeskyFactorization, - A::Adjoint{<:Number, <:Array}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Adjoint{<:Number, <:Array}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end end @@ -304,14 +305,14 @@ function do_factorization(alg::LDLtFactorization, A, b, u) end function init_cacheval(alg::LDLtFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) nothing end function init_cacheval(alg::LDLtFactorization, A::SymTridiagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.ldlt_instance(convert(AbstractMatrix, A)) end @@ -340,30 +341,30 @@ function do_factorization(alg::SVDFactorization, A, b, u) end function init_cacheval(alg::SVDFactorization, A::Matrix, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.svd_instance(convert(AbstractMatrix, A)) end const PREALLOCATED_SVD = ArrayInterface.svd_instance(rand(1, 1)) function init_cacheval(alg::SVDFactorization, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) PREALLOCATED_SVD end function init_cacheval(alg::SVDFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end @static if VERSION < v"1.7-" function init_cacheval(alg::SVDFactorization, - A::Union{Diagonal, SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{Diagonal, SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end end @@ -391,24 +392,25 @@ function do_factorization(alg::BunchKaufmanFactorization, A, b, u) end function init_cacheval(alg::BunchKaufmanFactorization, A::Symmetric{<:Number, <:Matrix}, b, - u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.bunchkaufman_instance(convert(AbstractMatrix, A)) end -const PREALLOCATED_BUNCHKAUFMAN = ArrayInterface.bunchkaufman_instance(Symmetric(rand(1, 1))) +const PREALLOCATED_BUNCHKAUFMAN = ArrayInterface.bunchkaufman_instance(Symmetric(rand(1, + 1))) function init_cacheval(alg::BunchKaufmanFactorization, - A::Symmetric{Float64, Matrix{Float64}}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Symmetric{Float64, Matrix{Float64}}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) PREALLOCATED_BUNCHKAUFMAN end function init_cacheval(alg::BunchKaufmanFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end @@ -446,165 +448,165 @@ function do_factorization(alg::GenericFactorization, A, b, u) end function init_cacheval(alg::GenericFactorization{typeof(lu)}, A, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.lu_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::GenericFactorization{typeof(lu!)}, A, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.lu_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::GenericFactorization{typeof(lu)}, - A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.lu_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(lu!)}, - A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.lu_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(lu)}, A::Diagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end function init_cacheval(alg::GenericFactorization{typeof(lu)}, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.lu_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(lu!)}, A::Diagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end function init_cacheval(alg::GenericFactorization{typeof(lu!)}, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.lu_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(qr)}, A, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.qr_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::GenericFactorization{typeof(qr!)}, A, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.qr_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::GenericFactorization{typeof(qr)}, - A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.qr_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(qr!)}, - A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.qr_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(qr)}, A::Diagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end function init_cacheval(alg::GenericFactorization{typeof(qr)}, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.qr_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(qr!)}, A::Diagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end function init_cacheval(alg::GenericFactorization{typeof(qr!)}, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.qr_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(svd)}, A, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.svd_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::GenericFactorization{typeof(svd!)}, A, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.svd_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::GenericFactorization{typeof(svd)}, - A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.svd_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(svd!)}, - A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.svd_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(svd)}, A::Diagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end function init_cacheval(alg::GenericFactorization{typeof(svd)}, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.svd_instance(A) end function init_cacheval(alg::GenericFactorization{typeof(svd!)}, A::Diagonal, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end function init_cacheval(alg::GenericFactorization{typeof(svd!)}, A::Tridiagonal, b, u, Pl, - Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.svd_instance(A) end function init_cacheval(alg::GenericFactorization, A::Diagonal, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end function init_cacheval(alg::GenericFactorization, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.lu_instance(A) end function init_cacheval(alg::GenericFactorization, A::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(A) end function init_cacheval(alg::Union{GenericFactorization{typeof(bunchkaufman!)}, - GenericFactorization{typeof(bunchkaufman)}}, - A::Union{Hermitian, Symmetric}, b, u, Pl, Pr, maxiters::Int, abstol, - reltol, verbose::Bool, assumptions::OperatorAssumptions) + GenericFactorization{typeof(bunchkaufman)}}, + A::Union{Hermitian, Symmetric}, b, u, Pl, Pr, maxiters::Int, abstol, + reltol, verbose::Bool, assumptions::OperatorAssumptions) BunchKaufman(A.data, Array(1:size(A, 1)), A.uplo, true, false, 0) end function init_cacheval(alg::Union{GenericFactorization{typeof(bunchkaufman!)}, - GenericFactorization{typeof(bunchkaufman)}}, - A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + GenericFactorization{typeof(bunchkaufman)}}, + A::StridedMatrix{<:LinearAlgebra.BlasFloat}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) if eltype(A) <: Complex return bunchkaufman!(Hermitian(A)) else @@ -617,18 +619,18 @@ end # Cholesky needs the posdef matrix, for GenericFactorization assume structure is needed function init_cacheval(alg::Union{GenericFactorization{typeof(cholesky)}, - GenericFactorization{typeof(cholesky!)}}, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + GenericFactorization{typeof(cholesky!)}}, A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) newA = copy(convert(AbstractMatrix, A)) do_factorization(alg, newA, b, u) end function init_cacheval(alg::Union{GenericFactorization}, - A::Union{Hermitian{T, <:SparseMatrixCSC}, - Symmetric{T, <:SparseMatrixCSC}}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T} + A::Union{Hermitian{T, <:SparseMatrixCSC}, + Symmetric{T, <:SparseMatrixCSC}}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T} newA = copy(convert(AbstractMatrix, A)) do_factorization(alg, newA, b, u) end @@ -657,46 +659,47 @@ end @static if VERSION < v"1.9.0-DEV.1622" const PREALLOCATED_UMFPACK = SuiteSparse.UMFPACK.UmfpackLU(C_NULL, C_NULL, 0, 0, - [0], Int64[], Float64[], 0) + [0], Int64[], Float64[], 0) finalizer(SuiteSparse.UMFPACK.umfpack_free_symbolic, PREALLOCATED_UMFPACK) else const PREALLOCATED_UMFPACK = SuiteSparse.UMFPACK.UmfpackLU(SparseMatrixCSC(0, 0, [1], - Int64[], - Float64[])) + Int64[], + Float64[])) end function init_cacheval(alg::UMFPACKFactorization, - A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) nothing end function init_cacheval(alg::UMFPACKFactorization, A::SparseMatrixCSC{Float64, Int}, b, u, - Pl, Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) PREALLOCATED_UMFPACK end -function init_cacheval(alg::UMFPACKFactorization, A::AbstractSparseArray, b, u, Pl, Pr, maxiters::Int, abstol, - reltol, - verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval(alg::UMFPACKFactorization, A::AbstractSparseArray, b, u, Pl, Pr, + maxiters::Int, abstol, + reltol, + verbose::Bool, assumptions::OperatorAssumptions) A = convert(AbstractMatrix, A) zerobased = SparseArrays.getcolptr(A)[1] == 0 @static if VERSION < v"1.9.0-DEV.1622" res = SuiteSparse.UMFPACK.UmfpackLU(C_NULL, C_NULL, size(A, 1), size(A, 2), - zerobased ? - copy(SparseArrays.getcolptr(A)) : - SuiteSparse.decrement(SparseArrays.getcolptr(A)), - zerobased ? copy(rowvals(A)) : - SuiteSparse.decrement(rowvals(A)), - copy(nonzeros(A)), 0) + zerobased ? + copy(SparseArrays.getcolptr(A)) : + SuiteSparse.decrement(SparseArrays.getcolptr(A)), + zerobased ? copy(rowvals(A)) : + SuiteSparse.decrement(rowvals(A)), + copy(nonzeros(A)), 0) finalizer(SuiteSparse.UMFPACK.umfpack_free_symbolic, res) return res else return SuiteSparse.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), - rowvals(A), nonzeros(A))) + rowvals(A), nonzeros(A))) end end @@ -711,11 +714,11 @@ function SciMLBase.solve!(cache::LinearCache, alg::UMFPACKFactorization; kwargs. SuiteSparse.decrement(SparseArrays.getrowval(A)) == @get_cacheval(cache, :UMFPACKFactorization).rowval) fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + nonzeros(A))) else fact = lu!(@get_cacheval(cache, :UMFPACKFactorization), - SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) end else fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))) @@ -746,28 +749,29 @@ Base.@kwdef struct KLUFactorization <: AbstractFactorization end const PREALLOCATED_KLU = KLU.KLUFactorization(SparseMatrixCSC(0, 0, [1], Int64[], - Float64[])) + Float64[])) function init_cacheval(alg::KLUFactorization, - A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) nothing end function init_cacheval(alg::KLUFactorization, A::SparseMatrixCSC{Float64, Int}, b, u, Pl, - Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + Pr, + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) PREALLOCATED_KLU end -function init_cacheval(alg::KLUFactorization, A::AbstractSparseArray, b, u, Pl, Pr, maxiters::Int, abstol, - reltol, - verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval(alg::KLUFactorization, A::AbstractSparseArray, b, u, Pl, Pr, + maxiters::Int, abstol, + reltol, + verbose::Bool, assumptions::OperatorAssumptions) A = convert(AbstractMatrix, A) return KLU.KLUFactorization(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + nonzeros(A))) end function SciMLBase.solve!(cache::LinearCache, alg::KLUFactorization; kwargs...) @@ -781,7 +785,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::KLUFactorization; kwargs...) cacheval.colptr && SuiteSparse.decrement(SparseArrays.getrowval(A)) == cacheval.rowval) fact = KLU.klu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + nonzeros(A))) else # If we have a cacheval already, run umfpack_symbolic to ensure the symbolic factorization exists # This won't recompute if it does. @@ -791,14 +795,14 @@ function SciMLBase.solve!(cache::LinearCache, alg::KLUFactorization; kwargs...) KLU.klu_factor!(cacheval) end fact = KLU.klu!(cacheval, - SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) end else # New fact each time since the sparsity pattern can change # and thus it needs to reallocate fact = KLU.klu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + nonzeros(A))) end cache.cacheval = fact cache.isfresh = false @@ -832,16 +836,17 @@ end const PREALLOCATED_CHOLMOD = cholesky(SparseMatrixCSC(0, 0, [1], Int64[], Float64[])) function init_cacheval(alg::CHOLMODFactorization, - A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) nothing end -function init_cacheval(alg::CHOLMODFactorization, A::Union{SparseMatrixCSC{T, Int}, Symmetric{T, SparseMatrixCSC{T, Int}}}, b, u, - Pl, Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) where {T <: Union{Float32,Float64}} +function init_cacheval(alg::CHOLMODFactorization, + A::Union{SparseMatrixCSC{T, Int}, Symmetric{T, SparseMatrixCSC{T, Int}}}, b, u, + Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) where {T <: Union{Float32, Float64}} PREALLOCATED_CHOLMOD end @@ -883,36 +888,36 @@ function RFLUFactorization(; pivot = Val(true), thread = Val(true)) end function init_cacheval(alg::RFLUFactorization, A, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ipiv = Vector{LinearAlgebra.BlasInt}(undef, min(size(A)...)) ArrayInterface.lu_instance(convert(AbstractMatrix, A)), ipiv end function init_cacheval(alg::RFLUFactorization, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ipiv = Vector{LinearAlgebra.BlasInt}(undef, 0) PREALLOCATED_LU, ipiv end function init_cacheval(alg::RFLUFactorization, - A::Union{AbstractSparseArray, AbstractSciMLOperator}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::Union{AbstractSparseArray, AbstractSciMLOperator}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) nothing, nothing end @static if VERSION < v"1.7-" function init_cacheval(alg::RFLUFactorization, - A::Union{Diagonal, SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + A::Union{Diagonal, SymTridiagonal, Tridiagonal}, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) nothing, nothing end end function SciMLBase.solve!(cache::LinearCache, alg::RFLUFactorization{P, T}; - kwargs...) where {P, T} + kwargs...) where {P, T} A = cache.A A = convert(AbstractMatrix, A) fact, ipiv = @get_cacheval(cache, :RFLUFactorization) @@ -966,40 +971,40 @@ else end const PREALLOCATED_NORMALCHOLESKY = ArrayInterface.cholesky_instance(rand(1, 1), - normcholpivot) + normcholpivot) function init_cacheval(alg::NormalCholeskyFactorization, - A::Union{AbstractSparseArray, - Symmetric{<:Number, <:AbstractSparseArray}}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{AbstractSparseArray, + Symmetric{<:Number, <:AbstractSparseArray}}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.cholesky_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::NormalCholeskyFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.cholesky_instance(convert(AbstractMatrix, A), alg.pivot) end function init_cacheval(alg::NormalCholeskyFactorization, - A::Union{Diagonal, AbstractSciMLOperator}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{Diagonal, AbstractSciMLOperator}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end @static if VERSION < v"1.7-" function init_cacheval(alg::NormalCholeskyFactorization, - A::Union{Tridiagonal, SymTridiagonal}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{Tridiagonal, SymTridiagonal}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end end function SciMLBase.solve!(cache::LinearCache, alg::NormalCholeskyFactorization; - kwargs...) + kwargs...) A = cache.A A = convert(AbstractMatrix, A) if cache.isfresh @@ -1045,13 +1050,13 @@ default_alias_A(::NormalBunchKaufmanFactorization, ::Any, ::Any) = true default_alias_b(::NormalBunchKaufmanFactorization, ::Any, ::Any) = true function init_cacheval(alg::NormalBunchKaufmanFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.bunchkaufman_instance(convert(AbstractMatrix, A)) end function SciMLBase.solve!(cache::LinearCache, alg::NormalBunchKaufmanFactorization; - kwargs...) + kwargs...) A = cache.A A = convert(AbstractMatrix, A) if cache.isfresh @@ -1073,12 +1078,12 @@ A special implementation only for solving `Diagonal` matrices fast. struct DiagonalFactorization <: AbstractFactorization end function init_cacheval(alg::DiagonalFactorization, A, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) nothing end function SciMLBase.solve!(cache::LinearCache, alg::DiagonalFactorization; - kwargs...) + kwargs...) A = convert(AbstractMatrix, cache.A) if cache.u isa Vector && cache.b isa Vector @simd ivdep for i in eachindex(cache.u) @@ -1109,8 +1114,8 @@ this version does not allow for choice of pivoting method. struct FastLUFactorization <: AbstractFactorization end function init_cacheval(::FastLUFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ws = LUWs(A) return WorkspaceAndFactors(ws, ArrayInterface.lu_instance(convert(AbstractMatrix, A))) end @@ -1123,7 +1128,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::FastLUFactorization; kwargs.. # we will fail here if A is a different *size* than in a previous version of the same cache. # it may instead be desirable to resize the workspace. @set! ws_and_fact.factors = LinearAlgebra.LU(LAPACK.getrf!(ws_and_fact.workspace, - A)...) + A)...) cache.cacheval = ws_and_fact cache.isfresh = false end @@ -1153,39 +1158,39 @@ end @static if VERSION < v"1.7beta" function init_cacheval(alg::FastQRFactorization{Val{false}}, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ws = QRWYWs(A; blocksize = alg.blocksize) return WorkspaceAndFactors(ws, - ArrayInterface.qr_instance(convert(AbstractMatrix, A))) + ArrayInterface.qr_instance(convert(AbstractMatrix, A))) end function init_cacheval(::FastQRFactorization{Val{true}}, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ws = QRpWs(A) return WorkspaceAndFactors(ws, - ArrayInterface.qr_instance(convert(AbstractMatrix, A))) + ArrayInterface.qr_instance(convert(AbstractMatrix, A))) end else function init_cacheval(alg::FastQRFactorization{NoPivot}, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ws = QRWYWs(A; blocksize = alg.blocksize) return WorkspaceAndFactors(ws, - ArrayInterface.qr_instance(convert(AbstractMatrix, A))) + ArrayInterface.qr_instance(convert(AbstractMatrix, A))) end function init_cacheval(::FastQRFactorization{ColumnNorm}, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ws = QRpWs(A) return WorkspaceAndFactors(ws, - ArrayInterface.qr_instance(convert(AbstractMatrix, A))) + ArrayInterface.qr_instance(convert(AbstractMatrix, A))) end end function SciMLBase.solve!(cache::LinearCache, alg::FastQRFactorization{P}; - kwargs...) where {P} + kwargs...) where {P} A = cache.A A = convert(AbstractMatrix, A) ws_and_fact = @get_cacheval(cache, :FastQRFactorization) @@ -1199,10 +1204,10 @@ function SciMLBase.solve!(cache::LinearCache, alg::FastQRFactorization{P}; end if P === nopivot @set! ws_and_fact.factors = LinearAlgebra.QRCompactWY(LAPACK.geqrt!(ws_and_fact.workspace, - A)...) + A)...) else @set! ws_and_fact.factors = LinearAlgebra.QRPivoted(LAPACK.geqp3!(ws_and_fact.workspace, - A)...) + A)...) end cache.cacheval = ws_and_fact cache.isfresh = false @@ -1234,33 +1239,33 @@ Base.@kwdef struct SparspakFactorization <: AbstractFactorization end const PREALLOCATED_SPARSEPAK = sparspaklu(SparseMatrixCSC(0, 0, [1], Int64[], Float64[]), - factorize = false) + factorize = false) function init_cacheval(alg::SparspakFactorization, - A::Union{Matrix, Nothing, AbstractSciMLOperator}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + A::Union{Matrix, Nothing, AbstractSciMLOperator}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) nothing end function init_cacheval(::SparspakFactorization, A::SparseMatrixCSC{Float64, Int}, b, u, Pl, - Pr, maxiters::Int, abstol, - reltol, - verbose::Bool, assumptions::OperatorAssumptions) + Pr, maxiters::Int, abstol, + reltol, + verbose::Bool, assumptions::OperatorAssumptions) PREALLOCATED_SPARSEPAK end function init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, - reltol, - verbose::Bool, assumptions::OperatorAssumptions) + reltol, + verbose::Bool, assumptions::OperatorAssumptions) A = convert(AbstractMatrix, A) if typeof(A) <: SparseArrays.AbstractSparseArray return sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A)), - factorize = false) + nonzeros(A)), + factorize = false) else return sparspaklu(SparseMatrixCSC(0, 0, [1], Int64[], eltype(A)[]), - factorize = false) + factorize = false) end end @@ -1269,11 +1274,11 @@ function SciMLBase.solve!(cache::LinearCache, alg::SparspakFactorization; kwargs if cache.isfresh if cache.cacheval !== nothing && alg.reuse_symbolic fact = sparspaklu!(@get_cacheval(cache, :SparspakFactorization), - SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + nonzeros(A))) else fact = sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + nonzeros(A))) end cache.cacheval = fact cache.isfresh = false @@ -1287,7 +1292,7 @@ for alg in InteractiveUtils.subtypes(AbstractFactorization) maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) init_cacheval(alg, A.A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) end -end \ No newline at end of file +end diff --git a/src/factorization_sparse.jl b/src/factorization_sparse.jl index 57c5a4849..b633796dd 100644 --- a/src/factorization_sparse.jl +++ b/src/factorization_sparse.jl @@ -1,15 +1,15 @@ # Specialize QR for the non-square case # Missing ldiv! definitions: https://github.com/JuliaSparse/SparseArrays.jl/issues/242 function _ldiv!(x::Vector, - A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY, - SuiteSparse.SPQR.QRSparse, - SuiteSparse.CHOLMOD.Factor}, b::Vector) + A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY, + SuiteSparse.SPQR.QRSparse, + SuiteSparse.CHOLMOD.Factor}, b::Vector) x .= A \ b end function _ldiv!(x::AbstractVector, - A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY, - SuiteSparse.SPQR.QRSparse, - SuiteSparse.CHOLMOD.Factor}, b::AbstractVector) + A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY, + SuiteSparse.SPQR.QRSparse, + SuiteSparse.CHOLMOD.Factor}, b::AbstractVector) x .= A \ b end diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl index 211fc1ab4..a0e6018a6 100644 --- a/src/iterative_wrappers.jl +++ b/src/iterative_wrappers.jl @@ -19,10 +19,10 @@ struct KrylovJL{F, I, A, K} <: AbstractKrylovSubspaceMethod end function KrylovJL(args...; KrylovAlg = Krylov.gmres!, - gmres_restart = 0, window = 0, - kwargs...) + gmres_restart = 0, window = 0, + kwargs...) return KrylovJL(KrylovAlg, gmres_restart, window, - args, kwargs) + args, kwargs) end default_alias_A(::KrylovJL, ::Any, ::Any) = true @@ -170,7 +170,7 @@ end # zeroinit allows for init_cacheval to start by initing with A (0,0) function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions; zeroinit = true) + verbose::Bool, assumptions::OperatorAssumptions; zeroinit = true) KS = get_KrylovJL_solver(alg.KrylovAlg) if zeroinit @@ -225,8 +225,8 @@ end function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...) if cache.isfresh solver = init_cacheval(alg, cache.A, cache.b, cache.u, cache.Pl, cache.Pr, - cache.maxiters, cache.abstol, cache.reltol, cache.verbose, - cache.assumptions, zeroinit = false) + cache.maxiters, cache.abstol, cache.reltol, cache.verbose, + cache.assumptions, zeroinit = false) cache.cacheval = solver cache.isfresh = false end @@ -245,24 +245,24 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...) args = (@get_cacheval(cache, :KrylovJL_GMRES), cache.A, cache.b) kwargs = (atol = atol, rtol = rtol, itmax = itmax, verbose = verbose, - ldiv = true, history = true, alg.kwargs...) + ldiv = true, history = true, alg.kwargs...) if cache.cacheval isa Krylov.CgSolver N !== I && @warn "$(alg.KrylovAlg) doesn't support right preconditioning." Krylov.solve!(args...; M = M, - kwargs...) + kwargs...) elseif cache.cacheval isa Krylov.GmresSolver Krylov.solve!(args...; M = M, N = N, - kwargs...) + kwargs...) elseif cache.cacheval isa Krylov.BicgstabSolver Krylov.solve!(args...; M = M, N = N, - kwargs...) + kwargs...) elseif cache.cacheval isa Krylov.MinresSolver N !== I && @warn "$(alg.KrylovAlg) doesn't support right preconditioning." Krylov.solve!(args...; M = M, - kwargs...) + kwargs...) else Krylov.solve!(args...; kwargs...) end @@ -271,5 +271,5 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...) resid = stats.residuals |> last return SciMLBase.build_linear_solution(alg, cache.u, resid, cache; - iters = stats.niter) + iters = stats.niter) end diff --git a/src/simplelu.jl b/src/simplelu.jl index 66fc38aa9..8049b243d 100644 --- a/src/simplelu.jl +++ b/src/simplelu.jl @@ -143,6 +143,6 @@ function SciMLBase.solve!(cache::LinearCache, alg::SimpleLUFactorization; kwargs end function init_cacheval(alg::SimpleLUFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, - reltol, verbose::Bool, assumptions::OperatorAssumptions) + reltol, verbose::Bool, assumptions::OperatorAssumptions) LUSolver(convert(AbstractMatrix, A)) end diff --git a/src/solve_function.jl b/src/solve_function.jl index 567d4dfe5..e26c7c621 100644 --- a/src/solve_function.jl +++ b/src/solve_function.jl @@ -4,7 +4,7 @@ struct LinearSolveFunction{F} <: AbstractSolveFunction end function SciMLBase.solve!(cache::LinearCache, alg::LinearSolveFunction, - args...; kwargs...) + args...; kwargs...) @unpack A, b, u, p, isfresh, Pl, Pr, cacheval = cache @unpack solve_func = alg diff --git a/test/basictests.jl b/test/basictests.jl index 33b919083..a6efeac85 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -202,10 +202,10 @@ end @testset "Concrete Factorizations" begin for alg in (LUFactorization(), - QRFactorization(), - SVDFactorization(), - RFLUFactorization(), - LinearSolve.defaultalg(prob1.A, prob1.b)) + QRFactorization(), + SVDFactorization(), + RFLUFactorization(), + LinearSolve.defaultalg(prob1.A, prob1.b)) @testset "$alg" begin test_interface(alg, prob1, prob2) end @@ -214,14 +214,14 @@ end @testset "Generic Factorizations" begin for fact_alg in (lu, lu!, - qr, qr!, - cholesky, - #cholesky!, - #ldlt, ldlt!, - bunchkaufman, bunchkaufman!, - lq, lq!, - svd, svd!, - LinearAlgebra.factorize) + qr, qr!, + cholesky, + #cholesky!, + #ldlt, ldlt!, + bunchkaufman, bunchkaufman!, + lq, lq!, + svd, svd!, + LinearAlgebra.factorize) @testset "fact_alg = $fact_alg" begin alg = GenericFactorization(fact_alg = fact_alg) test_interface(alg, prob1, prob2) @@ -232,10 +232,10 @@ end @testset "KrylovJL" begin kwargs = (; gmres_restart = 5) for alg in (("Default", KrylovJL(kwargs...)), - ("CG", KrylovJL_CG(kwargs...)), - ("GMRES", KrylovJL_GMRES(kwargs...)), - # ("BICGSTAB",KrylovJL_BICGSTAB(kwargs...)), - ("MINRES", KrylovJL_MINRES(kwargs...))) + ("CG", KrylovJL_CG(kwargs...)), + ("GMRES", KrylovJL_GMRES(kwargs...)), + # ("BICGSTAB",KrylovJL_BICGSTAB(kwargs...)), + ("MINRES", KrylovJL_MINRES(kwargs...))) @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) end @@ -246,11 +246,11 @@ end @testset "IterativeSolversJL" begin kwargs = (; gmres_restart = 5) for alg in (("Default", IterativeSolversJL(kwargs...)), - ("CG", IterativeSolversJL_CG(kwargs...)), - ("GMRES", IterativeSolversJL_GMRES(kwargs...)) - # ("BICGSTAB",IterativeSolversJL_BICGSTAB(kwargs...)), - # ("MINRES",IterativeSolversJL_MINRES(kwargs...)), - ) + ("CG", IterativeSolversJL_CG(kwargs...)), + ("GMRES", IterativeSolversJL_GMRES(kwargs...)) + # ("BICGSTAB",IterativeSolversJL_BICGSTAB(kwargs...)), + # ("MINRES",IterativeSolversJL_MINRES(kwargs...)), + ) @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) end @@ -261,8 +261,8 @@ end @testset "KrylovKit" begin kwargs = (; gmres_restart = 5) for alg in (("Default", KrylovKitJL(kwargs...)), - ("CG", KrylovKitJL_CG(kwargs...)), - ("GMRES", KrylovKitJL_GMRES(kwargs...))) + ("CG", KrylovKitJL_CG(kwargs...)), + ("GMRES", KrylovKitJL_GMRES(kwargs...))) @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) end @@ -378,7 +378,7 @@ end @testset "LinearSolveFunction" begin function sol_func(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, - kwargs...) + kwargs...) if verbose == true println("out-of-place solve") end @@ -386,7 +386,7 @@ end end function sol_func!(A, b, u, p, newA, Pl, Pr, solverdata; verbose = true, - kwargs...) + kwargs...) if verbose == true println("in-place solve") end @@ -397,7 +397,7 @@ end prob2 = LinearProblem(A1, b1; u0 = x1) for alg in (LinearSolveFunction(sol_func), - LinearSolveFunction(sol_func!)) + LinearSolveFunction(sol_func!)) test_interface(alg, prob1, prob2) end end diff --git a/test/default_algs.jl b/test/default_algs.jl index 0735976e1..671dbe1f4 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -9,7 +9,7 @@ using LinearSolve, LinearAlgebra, SparseArrays, Test, JET LinearSolve.DefaultAlgorithmChoice.DiagonalFactorization @test LinearSolve.defaultalg(nothing, zeros(5), - LinearSolve.OperatorAssumptions(false)).alg === + LinearSolve.OperatorAssumptions(false)).alg === LinearSolve.DefaultAlgorithmChoice.QRFactorization @test LinearSolve.defaultalg(sprand(1000, 1000, 0.01), zeros(1000)).alg === diff --git a/test/hypretests.jl b/test/hypretests.jl index 512dda6e0..3cceb6765 100644 --- a/test/hypretests.jl +++ b/test/hypretests.jl @@ -1,6 +1,7 @@ using HYPRE -using HYPRE.LibHYPRE: HYPRE_BigInt, HYPRE_Complex, HYPRE_IJMatrixGetValues, - HYPRE_IJVectorGetValues, HYPRE_Int +using HYPRE.LibHYPRE: HYPRE_BigInt, + HYPRE_Complex, HYPRE_IJMatrixGetValues, + HYPRE_IJVectorGetValues, HYPRE_Int using LinearAlgebra using LinearSolve using MPI diff --git a/test/nonsquare.jl b/test/nonsquare.jl index 763004310..c678c17b3 100644 --- a/test/nonsquare.jl +++ b/test/nonsquare.jl @@ -38,8 +38,8 @@ solve(LinearProblem(A, b)).u; solve(LinearProblem(A, b), (LinearSolve.NormalCholeskyFactorization())).u; solve(LinearProblem(A, b), (LinearSolve.NormalBunchKaufmanFactorization())).u; solve(LinearProblem(A, b), - assumptions = (OperatorAssumptions(false; - condition = OperatorCondition.WellConditioned))).u; + assumptions = (OperatorAssumptions(false; + condition = OperatorCondition.WellConditioned))).u; A = sprandn(5000, 100, 0.1) b = randn(5000) @@ -47,5 +47,5 @@ b = randn(5000) solve(LinearProblem(A, b)).u; solve(LinearProblem(A, b), (LinearSolve.NormalCholeskyFactorization())).u; solve(LinearProblem(A, b), - assumptions = (OperatorAssumptions(false; - condition = OperatorCondition.WellConditioned))).u; + assumptions = (OperatorAssumptions(false; + condition = OperatorCondition.WellConditioned))).u; diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index 512eed18e..efc1f4c23 100644 --- a/test/pardiso/pardiso.jl +++ b/test/pardiso/pardiso.jl @@ -2,9 +2,9 @@ using LinearSolve, SparseArrays, Random import Pardiso A1 = sparse([1.0 0 -2 3 - 0 5 1 2 - -2 1 4 -7 - 3 2 -7 5]) + 0 5 1 2 + -2 1 4 -7 + 3 2 -7 5]) b1 = rand(4) prob1 = LinearProblem(A1, b1) diff --git a/test/resolve.jl b/test/resolve.jl index 8648aa44d..936014459 100644 --- a/test/resolve.jl +++ b/test/resolve.jl @@ -42,25 +42,25 @@ linsolve.A = A @test solve!(linsolve).u ≈ [1.0, 0.5] A = Symmetric([1.0 2.0 - 2.0 1.0]) + 2.0 1.0]) b = [1.0, 2.0] prob = LinearProblem(A, b) linsolve = init(prob, BunchKaufmanFactorization(), alias_A = false, alias_b = false) @test solve!(linsolve).u ≈ [1.0, 0.0] @test solve!(linsolve).u ≈ [1.0, 0.0] A = Symmetric([1.0 2.0 - 2.0 1.0]) + 2.0 1.0]) linsolve.A = A @test solve!(linsolve).u ≈ [1.0, 0.0] A = Symmetric([1.0 2.0 - 2.0 1.0]) + 2.0 1.0]) b = [1.0, 2.0] prob = LinearProblem(A, b) linsolve = init(prob, CholeskyFactorization(), alias_A = false, alias_b = false) @test solve!(linsolve).u ≈ [1.0, 0.0] @test solve!(linsolve).u ≈ [1.0, 0.0] A = Symmetric([1.0 2.0 - 2.0 1.0]) + 2.0 1.0]) b = [1.0, 2.0] @test solve!(linsolve).u ≈ [1.0, 0.0]