diff --git a/Project.toml b/Project.toml index bbf6af778..33cabb12a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,149 +1,149 @@ -name = "LinearSolve" -uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -authors = ["SciML"] -version = "2.28.0" - -[deps] -ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" -DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" -FastLapackInterface = "29a986be-02c6-4525-aec4-84b980013641" -GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -KLU = "ef3ab10e-7fda-4108-b977-705223b18434" -Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" -LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" -Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7" -Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" -PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" -Preferences = "21216c6a-2e73-6563-6e65-726566657250" -RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" -Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" -StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" -UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" - -[weakdeps] -BandedMatrices = "aae01518-5342-5314-be14-df237901396f" -BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e" -Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" -EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" -FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" -HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" -IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" -KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -Metal = "dde4c033-4e86-420c-a63e-0dd931031962" -Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" -RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" - -[extensions] -LinearSolveBandedMatricesExt = "BandedMatrices" -LinearSolveBlockDiagonalsExt = "BlockDiagonals" -LinearSolveCUDAExt = "CUDA" -LinearSolveCUDSSExt = "CUDSS" -LinearSolveEnzymeExt = ["Enzyme", "EnzymeCore"] -LinearSolveFastAlmostBandedMatricesExt = ["FastAlmostBandedMatrices"] -LinearSolveHYPREExt = "HYPRE" -LinearSolveIterativeSolversExt = "IterativeSolvers" -LinearSolveKernelAbstractionsExt = "KernelAbstractions" -LinearSolveKrylovKitExt = "KrylovKit" -LinearSolveMetalExt = "Metal" -LinearSolvePardisoExt = "Pardiso" -LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools" - -[compat] -AllocCheck = "0.1" -Aqua = "0.8" -ArrayInterface = "7.7" -BandedMatrices = "1.5" -BlockDiagonals = "0.1.42" -CUDA = "5" -CUDSS = "0.1" -ChainRulesCore = "1.22" -ConcreteStructs = "0.2.3" -DocStringExtensions = "0.9.3" -EnumX = "1.0.4" -Enzyme = "0.11.15" -EnzymeCore = "0.6.5" -FastAlmostBandedMatrices = "0.1" -FastLapackInterface = "2" -FiniteDiff = "2.22" -ForwardDiff = "0.10.36" -GPUArraysCore = "0.1.6" -HYPRE = "1.4.0" -InteractiveUtils = "1.10" -IterativeSolvers = "0.9.3" -JET = "0.8.28" -KLU = "0.6" -KernelAbstractions = "0.9.16" -Krylov = "0.9" -KrylovKit = "0.6" -LazyArrays = "1.8" -Libdl = "1.10" -LinearAlgebra = "1.10" -MPI = "0.20" -Markdown = "1.10" -Metal = "1" -MultiFloats = "1" -Pardiso = "0.5" -Pkg = "1" -PrecompileTools = "1.2" -Preferences = "1.4" -Random = "1" -RecursiveArrayTools = "3.8" -RecursiveFactorization = "0.2.14" -Reexport = "1" -SafeTestsets = "0.1" -SciMLBase = "2.26.3" -SciMLOperators = "0.3.7" -Setfield = "1" -SparseArrays = "1.10" -Sparspak = "0.3.6" -StableRNGs = "1" -StaticArrays = "1.5" -StaticArraysCore = "1.4.2" -Test = "1" -UnPack = "1" -Zygote = "0.6.69" -julia = "1.10" - -[extras] -AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" -Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" -BandedMatrices = "aae01518-5342-5314-be14-df237901396f" -BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" -Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" -FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" -FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" -JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" -KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" -MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" -Metal = "dde4c033-4e86-420c-a63e-0dd931031962" -MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" -SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" - -[targets] -test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck", "StableRNGs", "Zygote"] +name = "LinearSolve" +uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +authors = ["SciML"] +version = "2.29.0" + +[deps] +ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" +FastLapackInterface = "29a986be-02c6-4525-aec4-84b980013641" +GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +KLU = "ef3ab10e-7fda-4108-b977-705223b18434" +Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" +Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MKL_jll = "856f044c-d86e-5d09-b602-aeab76dc8ba7" +Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" +PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" +RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" +Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" +StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" + +[weakdeps] +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" +BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e" +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" +EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" +FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" +HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" +IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +Metal = "dde4c033-4e86-420c-a63e-0dd931031962" +Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" + +[extensions] +LinearSolveBandedMatricesExt = "BandedMatrices" +LinearSolveBlockDiagonalsExt = "BlockDiagonals" +LinearSolveCUDAExt = "CUDA" +LinearSolveCUDSSExt = "CUDSS" +LinearSolveEnzymeExt = ["Enzyme", "EnzymeCore"] +LinearSolveFastAlmostBandedMatricesExt = ["FastAlmostBandedMatrices"] +LinearSolveHYPREExt = "HYPRE" +LinearSolveIterativeSolversExt = "IterativeSolvers" +LinearSolveKernelAbstractionsExt = "KernelAbstractions" +LinearSolveKrylovKitExt = "KrylovKit" +LinearSolveMetalExt = "Metal" +LinearSolvePardisoExt = "Pardiso" +LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools" + +[compat] +AllocCheck = "0.1" +Aqua = "0.8" +ArrayInterface = "7.7" +BandedMatrices = "1.5" +BlockDiagonals = "0.1.42" +CUDA = "5" +CUDSS = "0.1" +ChainRulesCore = "1.22" +ConcreteStructs = "0.2.3" +DocStringExtensions = "0.9.3" +EnumX = "1.0.4" +Enzyme = "0.11.15" +EnzymeCore = "0.6.5" +FastAlmostBandedMatrices = "0.1" +FastLapackInterface = "2" +FiniteDiff = "2.22" +ForwardDiff = "0.10.36" +GPUArraysCore = "0.1.6" +HYPRE = "1.4.0" +InteractiveUtils = "1.10" +IterativeSolvers = "0.9.3" +JET = "0.8.28" +KLU = "0.6" +KernelAbstractions = "0.9.16" +Krylov = "0.9" +KrylovKit = "0.6" +LazyArrays = "1.8" +Libdl = "1.10" +LinearAlgebra = "1.10" +MPI = "0.20" +Markdown = "1.10" +Metal = "1" +MultiFloats = "1" +Pardiso = "0.5" +Pkg = "1" +PrecompileTools = "1.2" +Preferences = "1.4" +Random = "1" +RecursiveArrayTools = "3.8" +RecursiveFactorization = "0.2.14" +Reexport = "1" +SafeTestsets = "0.1" +SciMLBase = "2.26.3" +SciMLOperators = "0.3.7" +Setfield = "1" +SparseArrays = "1.10" +Sparspak = "0.3.6" +StableRNGs = "1" +StaticArrays = "1.5" +StaticArraysCore = "1.4.2" +Test = "1" +UnPack = "1" +Zygote = "0.6.69" +julia = "1.10" + +[extras] +AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" +BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0" +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" +FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" +JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" +Metal = "dde4c033-4e86-420c-a63e-0dd931031962" +MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[targets] +test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck", "StableRNGs", "Zygote"] diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 57e6d5bd0..5df86c925 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -1,258 +1,265 @@ -module LinearSolve -if isdefined(Base, :Experimental) && - isdefined(Base.Experimental, Symbol("@max_methods")) - @eval Base.Experimental.@max_methods 1 -end - -import PrecompileTools - -PrecompileTools.@recompile_invalidations begin - using ArrayInterface - using RecursiveFactorization - using Base: cache_dependencies, Bool - using LinearAlgebra - using SparseArrays - using SparseArrays: AbstractSparseMatrixCSC, nonzeros, rowvals, getcolptr - using LazyArrays: @~, BroadcastArray - using SciMLBase: AbstractLinearAlgorithm - using SciMLOperators - using SciMLOperators: AbstractSciMLOperator, IdentityOperator - using Setfield - using UnPack - using KLU - using Sparspak - using FastLapackInterface - using DocStringExtensions - using EnumX - using Markdown - using ChainRulesCore - import InteractiveUtils - - import StaticArraysCore: StaticArray, SVector, MVector, SMatrix, MMatrix - - using LinearAlgebra: BlasInt, LU - using LinearAlgebra.LAPACK: require_one_based_indexing, - chkfinite, chkstride1, - @blasfunc, chkargsok - - import GPUArraysCore - import Preferences - import ConcreteStructs: @concrete - - # wrap - import Krylov - using SciMLBase - import Preferences -end - -const CRC = ChainRulesCore - -if Preferences.@load_preference("LoadMKL_JLL", true) - using MKL_jll - const usemkl = MKL_jll.is_available() -else - const usemkl = false -end - -using Reexport -@reexport using SciMLBase -using SciMLBase: _unwrap_val - -abstract type SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm end -abstract type AbstractFactorization <: SciMLLinearSolveAlgorithm end -abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end -abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end - -# Traits - -needs_concrete_A(alg::AbstractFactorization) = true -needs_concrete_A(alg::AbstractKrylovSubspaceMethod) = false -needs_concrete_A(alg::AbstractSolveFunction) = false - -# Util -is_underdetermined(x) = false -is_underdetermined(A::AbstractMatrix) = size(A, 1) < size(A, 2) -is_underdetermined(A::AbstractSciMLOperator) = size(A, 1) < size(A, 2) - -_isidentity_struct(A) = false -_isidentity_struct(λ::Number) = isone(λ) -_isidentity_struct(A::UniformScaling) = isone(A.λ) -_isidentity_struct(::SciMLOperators.IdentityOperator) = true - -# Dispatch Friendly way to check if an extension is loaded -__is_extension_loaded(::Val) = false - -# Check if a sparsity pattern has changed -pattern_changed(fact, A) = false - -function _fast_sym_givens! end - -# Code - -const INCLUDE_SPARSE = Preferences.@load_preference("include_sparse", Base.USE_GPL_LIBS) - -EnumX.@enumx DefaultAlgorithmChoice begin - LUFactorization - QRFactorization - DiagonalFactorization - DirectLdiv! - SparspakFactorization - KLUFactorization - UMFPACKFactorization - KrylovJL_GMRES - GenericLUFactorization - RFLUFactorization - LDLtFactorization - BunchKaufmanFactorization - CHOLMODFactorization - SVDFactorization - CholeskyFactorization - NormalCholeskyFactorization - AppleAccelerateLUFactorization - MKLLUFactorization - QRFactorizationPivoted - KrylovJL_CRAIGMR - KrylovJL_LSMR -end - -struct DefaultLinearSolver <: SciMLLinearSolveAlgorithm - alg::DefaultAlgorithmChoice.T -end - -const BLASELTYPES = Union{Float32, Float64, ComplexF32, ComplexF64} - -include("common.jl") -include("factorization.jl") -include("appleaccelerate.jl") -include("mkl.jl") -include("simplelu.jl") -include("simplegmres.jl") -include("iterative_wrappers.jl") -include("preconditioners.jl") -include("solve_function.jl") -include("default.jl") -include("init.jl") -include("extension_algs.jl") -include("adjoint.jl") -include("deprecated.jl") - -@generated function SciMLBase.solve!(cache::LinearCache, alg::AbstractFactorization; - kwargs...) - quote - if cache.isfresh - fact = do_factorization(alg, cache.A, cache.b, cache.u) - cache.cacheval = fact - cache.isfresh = false - end - y = _ldiv!(cache.u, @get_cacheval(cache, $(Meta.quot(defaultalg_symbol(alg)))), - cache.b) - - #= - retcode = if LinearAlgebra.issuccess(fact) - SciMLBase.ReturnCode.Success - else - SciMLBase.ReturnCode.Failure - end - SciMLBase.build_linear_solution(alg, y, nothing, cache; retcode = retcode) - =# - SciMLBase.build_linear_solution(alg, y, nothing, cache) - end -end - -@static if INCLUDE_SPARSE - include("factorization_sparse.jl") -end - -# Solver Specific Traits -## Needs Square Matrix -""" - needs_square_A(alg) - -Returns `true` if the algorithm requires a square matrix. -""" -needs_square_A(::Nothing) = false # Linear Solve automatically will use a correct alg! -needs_square_A(alg::SciMLLinearSolveAlgorithm) = true -for alg in (:QRFactorization, :FastQRFactorization, :NormalCholeskyFactorization, - :NormalBunchKaufmanFactorization) - @eval needs_square_A(::$(alg)) = false -end -for kralg in (Krylov.lsmr!, Krylov.craigmr!) - @eval needs_square_A(::KrylovJL{$(typeof(kralg))}) = false -end -for alg in (:LUFactorization, :FastLUFactorization, :SVDFactorization, - :GenericFactorization, :GenericLUFactorization, :SimpleLUFactorization, - :RFLUFactorization, :UMFPACKFactorization, :KLUFactorization, :SparspakFactorization, - :DiagonalFactorization, :CholeskyFactorization, :BunchKaufmanFactorization, - :CHOLMODFactorization, :LDLtFactorization, :AppleAccelerateLUFactorization, - :MKLLUFactorization, :MetalLUFactorization) - @eval needs_square_A(::$(alg)) = true -end - -const IS_OPENBLAS = Ref(true) -isopenblas() = IS_OPENBLAS[] - -const HAS_APPLE_ACCELERATE = Ref(false) -appleaccelerate_isavailable() = HAS_APPLE_ACCELERATE[] - -PrecompileTools.@compile_workload begin - A = rand(4, 4) - b = rand(4) - prob = LinearProblem(A, b) - sol = solve(prob) - sol = solve(prob, LUFactorization()) - sol = solve(prob, RFLUFactorization()) - sol = solve(prob, KrylovJL_GMRES()) -end - -@static if INCLUDE_SPARSE - PrecompileTools.@compile_workload begin - A = sprand(4, 4, 0.3) + I - b = rand(4) - prob = LinearProblem(A, b) - sol = solve(prob, KLUFactorization()) - sol = solve(prob, UMFPACKFactorization()) - end -end - -PrecompileTools.@compile_workload begin - A = sprand(4, 4, 0.3) + I - b = rand(4) - prob = LinearProblem(A * A', b) - sol = solve(prob) # in case sparspak is used as default - sol = solve(prob, SparspakFactorization()) -end - -ALREADY_WARNED_CUDSS = Ref{Bool}(false) -error_no_cudss_lu(A) = nothing -cudss_loaded(A) = false - -export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization, - 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, IterativeSolversJL_IDRS, - KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES - -export SimpleGMRES - -export HYPREAlgorithm -export CudaOffloadFactorization -export MKLPardisoFactorize, MKLPardisoIterate -export PardisoJL -export MKLLUFactorization -export AppleAccelerateLUFactorization -export MetalLUFactorization - -export OperatorAssumptions, OperatorCondition - -export LinearSolveAdjoint - -end +module LinearSolve +if isdefined(Base, :Experimental) && + isdefined(Base.Experimental, Symbol("@max_methods")) + @eval Base.Experimental.@max_methods 1 +end + +import PrecompileTools + +PrecompileTools.@recompile_invalidations begin + using ArrayInterface + using RecursiveFactorization + using Base: cache_dependencies, Bool + using LinearAlgebra + using SparseArrays + using SparseArrays: AbstractSparseMatrixCSC, nonzeros, rowvals, getcolptr + using LazyArrays: @~, BroadcastArray + using SciMLBase: AbstractLinearAlgorithm + using SciMLOperators + using SciMLOperators: AbstractSciMLOperator, IdentityOperator + using Setfield + using UnPack + using KLU + using Sparspak + using FastLapackInterface + using DocStringExtensions + using EnumX + using Markdown + using ChainRulesCore + import InteractiveUtils + + import StaticArraysCore: StaticArray, SVector, MVector, SMatrix, MMatrix + + using LinearAlgebra: BlasInt, LU + using LinearAlgebra.LAPACK: require_one_based_indexing, + chkfinite, chkstride1, + @blasfunc, chkargsok + + import GPUArraysCore + import Preferences + import ConcreteStructs: @concrete + + # wrap + import Krylov + using SciMLBase + import Preferences +end + +const CRC = ChainRulesCore + +if Preferences.@load_preference("LoadMKL_JLL", true) + using MKL_jll + const usemkl = MKL_jll.is_available() +else + const usemkl = false +end + +using Reexport +@reexport using SciMLBase +using SciMLBase: _unwrap_val + +abstract type SciMLLinearSolveAlgorithm <: SciMLBase.AbstractLinearAlgorithm end +abstract type AbstractFactorization <: SciMLLinearSolveAlgorithm end +abstract type AbstractKrylovSubspaceMethod <: SciMLLinearSolveAlgorithm end +abstract type AbstractSolveFunction <: SciMLLinearSolveAlgorithm end + +# Traits + +needs_concrete_A(alg::AbstractFactorization) = true +needs_concrete_A(alg::AbstractKrylovSubspaceMethod) = false +needs_concrete_A(alg::AbstractSolveFunction) = false + +# Util +is_underdetermined(x) = false +is_underdetermined(A::AbstractMatrix) = size(A, 1) < size(A, 2) +is_underdetermined(A::AbstractSciMLOperator) = size(A, 1) < size(A, 2) + +_isidentity_struct(A) = false +_isidentity_struct(λ::Number) = isone(λ) +_isidentity_struct(A::UniformScaling) = isone(A.λ) +_isidentity_struct(::SciMLOperators.IdentityOperator) = true + +# Dispatch Friendly way to check if an extension is loaded +__is_extension_loaded(::Val) = false + +# Check if a sparsity pattern has changed +pattern_changed(fact, A) = false + +function _fast_sym_givens! end + +# Code + +const INCLUDE_SPARSE = Preferences.@load_preference("include_sparse", Base.USE_GPL_LIBS) + +EnumX.@enumx DefaultAlgorithmChoice begin + LUFactorization + QRFactorization + DiagonalFactorization + DirectLdiv! + SparspakFactorization + KLUFactorization + UMFPACKFactorization + KrylovJL_GMRES + GenericLUFactorization + RFLUFactorization + LDLtFactorization + BunchKaufmanFactorization + CHOLMODFactorization + SVDFactorization + CholeskyFactorization + NormalCholeskyFactorization + AppleAccelerateLUFactorization + MKLLUFactorization + QRFactorizationPivoted + KrylovJL_CRAIGMR + KrylovJL_LSMR +end + +struct DefaultLinearSolver <: SciMLLinearSolveAlgorithm + alg::DefaultAlgorithmChoice.T +end + +const BLASELTYPES = Union{Float32, Float64, ComplexF32, ComplexF64} + +include("common.jl") +include("factorization.jl") +include("appleaccelerate.jl") +include("mkl.jl") +include("simplelu.jl") +include("simplegmres.jl") +include("iterative_wrappers.jl") +include("preconditioners.jl") +include("solve_function.jl") +include("default.jl") +include("init.jl") +include("extension_algs.jl") +include("adjoint.jl") +include("deprecated.jl") + +@inline function _notsuccessful(F::LinearAlgebra.QRCompactWY) + (m, n) = size(F) + U = view(F.factors, 1:min(m, n), 1:n) + return any(iszero, Iterators.reverse(@view U[diagind(U)])) +end +@inline _notsuccessful(F) = hasmethod(LinearAlgebra.issuccess, (typeof(F),)) ? + !LinearAlgebra.issuccess(F) : false + +@generated function SciMLBase.solve!(cache::LinearCache, alg::AbstractFactorization; + kwargs...) + quote + if cache.isfresh + fact = do_factorization(alg, cache.A, cache.b, cache.u) + cache.cacheval = fact + + # If factorization was not successful, return failure. Don't reset `isfresh` + if _notsuccessful(fact) + return SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Failure) + end + + cache.isfresh = false + end + + y = _ldiv!(cache.u, @get_cacheval(cache, $(Meta.quot(defaultalg_symbol(alg)))), + cache.b) + return SciMLBase.build_linear_solution(alg, y, nothing, cache) + end +end + +@static if INCLUDE_SPARSE + include("factorization_sparse.jl") +end + +# Solver Specific Traits +## Needs Square Matrix +""" + needs_square_A(alg) + +Returns `true` if the algorithm requires a square matrix. +""" +needs_square_A(::Nothing) = false # Linear Solve automatically will use a correct alg! +needs_square_A(alg::SciMLLinearSolveAlgorithm) = true +for alg in (:QRFactorization, :FastQRFactorization, :NormalCholeskyFactorization, + :NormalBunchKaufmanFactorization) + @eval needs_square_A(::$(alg)) = false +end +for kralg in (Krylov.lsmr!, Krylov.craigmr!) + @eval needs_square_A(::KrylovJL{$(typeof(kralg))}) = false +end +for alg in (:LUFactorization, :FastLUFactorization, :SVDFactorization, + :GenericFactorization, :GenericLUFactorization, :SimpleLUFactorization, + :RFLUFactorization, :UMFPACKFactorization, :KLUFactorization, :SparspakFactorization, + :DiagonalFactorization, :CholeskyFactorization, :BunchKaufmanFactorization, + :CHOLMODFactorization, :LDLtFactorization, :AppleAccelerateLUFactorization, + :MKLLUFactorization, :MetalLUFactorization) + @eval needs_square_A(::$(alg)) = true +end + +const IS_OPENBLAS = Ref(true) +isopenblas() = IS_OPENBLAS[] + +const HAS_APPLE_ACCELERATE = Ref(false) +appleaccelerate_isavailable() = HAS_APPLE_ACCELERATE[] + +PrecompileTools.@compile_workload begin + A = rand(4, 4) + b = rand(4) + prob = LinearProblem(A, b) + sol = solve(prob) + sol = solve(prob, LUFactorization()) + sol = solve(prob, RFLUFactorization()) + sol = solve(prob, KrylovJL_GMRES()) +end + +@static if INCLUDE_SPARSE + PrecompileTools.@compile_workload begin + A = sprand(4, 4, 0.3) + I + b = rand(4) + prob = LinearProblem(A, b) + sol = solve(prob, KLUFactorization()) + sol = solve(prob, UMFPACKFactorization()) + end +end + +PrecompileTools.@compile_workload begin + A = sprand(4, 4, 0.3) + I + b = rand(4) + prob = LinearProblem(A * A', b) + sol = solve(prob) # in case sparspak is used as default + sol = solve(prob, SparspakFactorization()) +end + +ALREADY_WARNED_CUDSS = Ref{Bool}(false) +error_no_cudss_lu(A) = nothing +cudss_loaded(A) = false + +export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization, + 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, IterativeSolversJL_IDRS, + KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES + +export SimpleGMRES + +export HYPREAlgorithm +export CudaOffloadFactorization +export MKLPardisoFactorize, MKLPardisoIterate +export PardisoJL +export MKLLUFactorization +export AppleAccelerateLUFactorization +export MetalLUFactorization + +export OperatorAssumptions, OperatorCondition + +export LinearSolveAdjoint + +end diff --git a/test/default_algs.jl b/test/default_algs.jl index 8d3ac261e..e6d69fe52 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -112,3 +112,35 @@ prob = LinearProblem(funcop, b) sol1 = solve(prob) sol2 = solve(prob, LinearSolve.KrylovJL_CRAIGMR()) @test sol1.u == sol2.u + +# Default for Underdetermined problem but the size is a long rectangle +A = [2.0 1.0 + 0.0 0.0 + 0.0 0.0] +b = [1.0, 0.0, 0.0] +prob = LinearProblem(A, b) +sol = solve(prob) + +@test !SciMLBase.successful_retcode(sol.retcode) + +## Show that we cannot select a default alg once by checking the rank, since it might change +## later in the cache +## Common occurrence for iterative nonlinear solvers using linear solve +A = [2.0 1.0 + 1.0 1.0 + 0.0 0.0] +b = [1.0, 1.0, 0.0] +prob = LinearProblem(A, b) + +cache = init(prob) +sol = solve!(cache) + +@test sol.u ≈ [0.0, 1.0] + +cache.A = [2.0 1.0 + 0.0 0.0 + 0.0 0.0] + +sol = solve!(cache) + +@test !SciMLBase.successful_retcode(sol.retcode)