From 47cb4028d61c11c8bb6eb4d649111f293c5e9ba4 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Mon, 18 Dec 2023 15:24:22 -0500 Subject: [PATCH 01/55] Patch the Static Arrays tests --- Project.toml | 3 ++- test/static_arrays.jl | 16 +++++++++------- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index ba9907272..77c065549 100644 --- a/Project.toml +++ b/Project.toml @@ -131,8 +131,9 @@ 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" [targets] -test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck"] +test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck", "StableRNGs"] diff --git a/test/static_arrays.jl b/test/static_arrays.jl index 72676d482..300c86c25 100644 --- a/test/static_arrays.jl +++ b/test/static_arrays.jl @@ -1,8 +1,10 @@ -using LinearSolve, StaticArrays, LinearAlgebra, Test +using LinearSolve, StaticArrays, LinearAlgebra, Test, StableRNGs using AllocCheck -A = SMatrix{5, 5}(Hermitian(rand(5, 5) + I)) -b = SVector{5}(rand(5)) +rng = StableRNG(0) + +A = SMatrix{5, 5}(Hermitian(rand(rng, 5, 5) + I)) +b = SVector{5}(rand(rng, 5)) @check_allocs __solve_no_alloc(A, b, alg) = solve(LinearProblem(A, b), alg) @@ -27,16 +29,16 @@ for alg in (nothing, LUFactorization(), SVDFactorization(), CholeskyFactorizatio @test norm(A * sol .- b) < 1e-10 end -A = SMatrix{7, 5}(rand(7, 5)) -b = SVector{7}(rand(7)) +A = SMatrix{7, 5}(rand(rng, 7, 5)) +b = SVector{7}(rand(rng, 7)) for alg in (nothing, SVDFactorization(), KrylovJL_LSMR()) @inferred solve(LinearProblem(A, b), alg) @test_nowarn solve(LinearProblem(A, b), alg) end -A = SMatrix{5, 7}(rand(5, 7)) -b = SVector{5}(rand(5)) +A = SMatrix{5, 7}(rand(rng, 5, 7)) +b = SVector{5}(rand(rng, 5)) for alg in (nothing, SVDFactorization(), KrylovJL_LSMR()) @inferred solve(LinearProblem(A, b), alg) From ff28502a08123cfcb50d27b709c8920a65ab9f10 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Dec 2023 02:05:38 -0500 Subject: [PATCH 02/55] Update factorization.jl --- src/factorization.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/factorization.jl b/src/factorization.jl index 695909ae1..6a6e91889 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -13,6 +13,11 @@ _ldiv!(x, A, b) = ldiv!(x, A, b) _ldiv!(x, A, b::SVector) = (x .= A \ b) _ldiv!(::SVector, A, b::SVector) = (A \ b) _ldiv!(::SVector, A, b) = (A \ b) +_ldiv!(::StaticArraysCore.SVector, + A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, + b::AbstractVector) = (A \ b) +_ldiv!(::SVector, A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, + b::SVector) = (A \ b) function _ldiv!(x::Vector, A::Factorization, b::Vector) # workaround https://github.com/JuliaLang/julia/issues/43507 From 8337af2618086cdb1d35b5829dc482ddc23e6aa0 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Dec 2023 02:13:28 -0500 Subject: [PATCH 03/55] Update factorization.jl --- src/factorization.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 6a6e91889..695909ae1 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -13,11 +13,6 @@ _ldiv!(x, A, b) = ldiv!(x, A, b) _ldiv!(x, A, b::SVector) = (x .= A \ b) _ldiv!(::SVector, A, b::SVector) = (A \ b) _ldiv!(::SVector, A, b) = (A \ b) -_ldiv!(::StaticArraysCore.SVector, - A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, - b::AbstractVector) = (A \ b) -_ldiv!(::SVector, A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, - b::SVector) = (A \ b) function _ldiv!(x::Vector, A::Factorization, b::Vector) # workaround https://github.com/JuliaLang/julia/issues/43507 From c2374ddaf9e69c77ce95501e51c278a10cdfdc40 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 20 Dec 2023 02:13:50 -0500 Subject: [PATCH 04/55] Update factorization_sparse.jl --- src/factorization_sparse.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/factorization_sparse.jl b/src/factorization_sparse.jl index 357e46446..1bc0ac261 100644 --- a/src/factorization_sparse.jl +++ b/src/factorization_sparse.jl @@ -13,3 +13,10 @@ function _ldiv!(x::AbstractVector, SparseArrays.CHOLMOD.Factor}, b::AbstractVector) x .= A \ b end + +# Ambiguity removal +_ldiv!(::SVector, + A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, + b::AbstractVector) = (A \ b) +_ldiv!(::SVector, A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, + b::SVector) = (A \ b) From 593cd43e264ff940a00b612a5cc109256ba63898 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 21 Dec 2023 14:54:49 -0500 Subject: [PATCH 05/55] Update simplegmres.jl --- src/simplegmres.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/simplegmres.jl b/src/simplegmres.jl index 03b8014f2..76682e547 100644 --- a/src/simplegmres.jl +++ b/src/simplegmres.jl @@ -37,9 +37,16 @@ struct SimpleGMRES{UBD} <: AbstractKrylovSubspaceMethod blocksize::Int warm_start::Bool + function SimpleGMRES{UBD}(; restart::Bool = true, blocksize::Int = 0, + warm_start::Bool = false, memory::Int = 20) where {UBD} + UBD && @assert blocksize > 0 + return new{UBD}(restart, memory, blocksize, warm_start) + end + function SimpleGMRES(; restart::Bool = true, blocksize::Int = 0, warm_start::Bool = false, memory::Int = 20) - return new{blocksize > 0}(restart, memory, blocksize, warm_start) + return SimpleGMRES{blocksize > 0}(; restart, memory, blocksize, + warm_start) end end From 58df085ad33e2ab8e38661bd2eb3d5bfab8881df Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Fri, 22 Dec 2023 10:29:57 -0500 Subject: [PATCH 06/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ba9907272..9cdfba710 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.21.1" +version = "2.21.2" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From 28e8a720ae9bb36b1b7489df09d74b4643a199ad Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Fri, 22 Dec 2023 16:39:51 +0100 Subject: [PATCH 07/55] downgrade CI --- .github/workflows/Downgrade.yml | 31 +++++++++++++++++ Project.toml | 60 ++++++++++++++++----------------- src/LinearSolve.jl | 1 - src/init.jl | 18 +--------- 4 files changed, 61 insertions(+), 49 deletions(-) create mode 100644 .github/workflows/Downgrade.yml diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml new file mode 100644 index 000000000..d365424e3 --- /dev/null +++ b/.github/workflows/Downgrade.yml @@ -0,0 +1,31 @@ +name: Downgrade +on: + pull_request: + branches: + - main + paths-ignore: + - 'docs/**' + push: + branches: + - main + paths-ignore: + - 'docs/**' +jobs: + test: + runs-on: ubuntu-latest + strategy: + matrix: + version: ['1'] + group: + - Core + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + - uses: cjdoris/julia-downgrade-compat-action@v1 +# if: ${{ matrix.version == '1.6' }} + with: + skip: Pkg,TOML + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 \ No newline at end of file diff --git a/Project.toml b/Project.toml index ba9907272..d92da9de9 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,6 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" Preferences = "21216c6a-2e73-6563-6e65-726566657250" RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -Requires = "ae029012-a4dd-5104-9daa-d747884805df" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" @@ -61,51 +60,50 @@ LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools" [compat] AllocCheck = "0.1" Aqua = "0.8" -ArrayInterface = "7.4.11" -BandedMatrices = "1" -BlockDiagonals = "0.1" +ArrayInterface = "7.7" +BandedMatrices = "1.4" +BlockDiagonals = "0.1.42" CUDA = "5" -ConcreteStructs = "0.2" -DocStringExtensions = "0.9" +ConcreteStructs = "0.2.3" +DocStringExtensions = "0.9.3" EnumX = "1" -Enzyme = "0.11" -EnzymeCore = "0.6" +Enzyme = "0.11.11" +EnzymeCore = "0.6.4" FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" -FiniteDiff = "2" -ForwardDiff = "0.10" -GPUArraysCore = "0.1" -HYPRE = "1.4.0" -InteractiveUtils = "1.6" +FiniteDiff = "2.21" +ForwardDiff = "0.10.36" +GPUArraysCore = "0.1.5" +HYPRE = "1.5.0" +InteractiveUtils = "1.9" IterativeSolvers = "0.9.3" -JET = "0.8" -KLU = "0.3.0, 0.4" +JET = "0.8.21" +KLU = "0.4.1" KernelAbstractions = "0.9" -Krylov = "0.9" +Krylov = "0.9.5" KrylovKit = "0.6" -Libdl = "1.6" +Libdl = "1.9" LinearAlgebra = "1.9" -MPI = "0.20" +MPI = "0.20.19" Metal = "0.5" MultiFloats = "1" Pardiso = "0.5" Pkg = "1" -PrecompileTools = "1" -Preferences = "1" +PrecompileTools = "1.2" +Preferences = "1.4" Random = "1" -RecursiveArrayTools = "2, 3" -RecursiveFactorization = "0.2.8" -Reexport = "1" -Requires = "1" +RecursiveArrayTools = "3.2" +RecursiveFactorization = "0.2.21" +Reexport = "1.2" SafeTestsets = "0.1" -SciMLBase = "2" -SciMLOperators = "0.3" -Setfield = "1" +SciMLBase = "2.11" +SciMLOperators = "0.3.7" +Setfield = "1.1" SparseArrays = "1.9" -Sparspak = "0.3.6" -StaticArrays = "1" -StaticArraysCore = "1" -Test = "1" +Sparspak = "0.3.9" +StaticArrays = "1.8" +StaticArraysCore = "1.4" +Test = "1.9" UnPack = "1" julia = "1.9" diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 572e310e9..fc039d13f 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -23,7 +23,6 @@ PrecompileTools.@recompile_invalidations begin using FastLapackInterface using DocStringExtensions using EnumX - using Requires import InteractiveUtils import StaticArraysCore: StaticArray, SVector, MVector, SMatrix, MMatrix diff --git a/src/init.jl b/src/init.jl index 1291b4556..e10960f1d 100644 --- a/src/init.jl +++ b/src/init.jl @@ -1,19 +1,3 @@ function __init__() - @static if VERSION < v"1.7beta" - blas = BLAS.vendor() - IS_OPENBLAS[] = blas == :openblas64 || blas == :openblas - else - IS_OPENBLAS[] = occursin("openblas", BLAS.get_config().loaded_libs[1].libname) - end - @static if !isdefined(Base, :get_extension) - @require IterativeSolvers="b77e0a4c-d291-57a0-90e8-8db25a27a240" begin - include("../ext/LinearSolveIterativeSolversExt.jl") - end - @require KrylovKit="0b1a1467-8014-51b9-945f-bf0ae24f4b77" begin - include("../ext/LinearSolveKrylovKitExt.jl") - end - @require Enzyme="7da242da-08ed-463a-9acd-ee780be4f1d9" begin - include("../ext/LinearSolveEnzymeExt.jl") - end - end + IS_OPENBLAS[] = occursin("openblas", BLAS.get_config().loaded_libs[1].libname) end From 4cc2075395f5ac27fd2f4bef86aa442f65cac564 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Fri, 22 Dec 2023 18:44:03 +0100 Subject: [PATCH 08/55] loosen bounds --- Project.toml | 52 ++++++++++++++++++++++++++-------------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index d92da9de9..1ba11e509 100644 --- a/Project.toml +++ b/Project.toml @@ -60,50 +60,50 @@ LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools" [compat] AllocCheck = "0.1" Aqua = "0.8" -ArrayInterface = "7.7" -BandedMatrices = "1.4" +ArrayInterface = "7.5" +BandedMatrices = "1" BlockDiagonals = "0.1.42" CUDA = "5" -ConcreteStructs = "0.2.3" -DocStringExtensions = "0.9.3" +ConcreteStructs = "0.2" +DocStringExtensions = "0.9" EnumX = "1" -Enzyme = "0.11.11" -EnzymeCore = "0.6.4" +Enzyme = "0.11.10" +EnzymeCore = "0.6.2" FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" -FiniteDiff = "2.21" -ForwardDiff = "0.10.36" -GPUArraysCore = "0.1.5" -HYPRE = "1.5.0" +FiniteDiff = "2.18" +ForwardDiff = "0.10.13" +GPUArraysCore = "0.1" +HYPRE = "1.4.0" InteractiveUtils = "1.9" IterativeSolvers = "0.9.3" -JET = "0.8.21" -KLU = "0.4.1" +JET = "0.8" +KLU = "0.3.0, 0.4" KernelAbstractions = "0.9" -Krylov = "0.9.5" +Krylov = "0.9" KrylovKit = "0.6" Libdl = "1.9" LinearAlgebra = "1.9" -MPI = "0.20.19" +MPI = "0.20" Metal = "0.5" MultiFloats = "1" Pardiso = "0.5" Pkg = "1" -PrecompileTools = "1.2" +PrecompileTools = "1.1" Preferences = "1.4" Random = "1" -RecursiveArrayTools = "3.2" -RecursiveFactorization = "0.2.21" -Reexport = "1.2" +RecursiveArrayTools = "2.38, 3" +RecursiveFactorization = "0.2.14" +Reexport = "1" SafeTestsets = "0.1" -SciMLBase = "2.11" -SciMLOperators = "0.3.7" -Setfield = "1.1" +SciMLBase = "2" +SciMLOperators = "0.3" +Setfield = "1" SparseArrays = "1.9" -Sparspak = "0.3.9" -StaticArrays = "1.8" -StaticArraysCore = "1.4" -Test = "1.9" +Sparspak = "0.3.6" +StaticArrays = "1.5" +StaticArraysCore = "1.2" +Test = "1" UnPack = "1" julia = "1.9" @@ -133,4 +133,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck"] +test = ["Aqua", "Test", "IterativeSolvers", "InteractiveUtils", "JET", "KrylovKit", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff", "HYPRE", "MPI", "BlockDiagonals", "Enzyme", "FiniteDiff", "BandedMatrices", "FastAlmostBandedMatrices", "StaticArrays", "AllocCheck"] \ No newline at end of file From 80af62653a43d321167af35741d6d05518ceac9f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 10:25:52 -0500 Subject: [PATCH 09/55] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 6cffe438c..547f611bb 100644 --- a/Project.toml +++ b/Project.toml @@ -101,6 +101,7 @@ SciMLOperators = "0.3" Setfield = "1" SparseArrays = "1.9" Sparspak = "0.3.6" +StableRNGs = "1" StaticArrays = "1.5" StaticArraysCore = "1.2" Test = "1" From f10bc7088409a22f3742a55653fc08886386b32e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 09:52:05 -0500 Subject: [PATCH 10/55] Fix potential uninitialized in Symmetric building of NormalChol cache --- src/factorization.jl | 2 +- test/default_algs.jl | 41 +++++++++++++++++++++-------------------- 2 files changed, 22 insertions(+), 21 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 695909ae1..7c261ec6b 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -1010,7 +1010,7 @@ function init_cacheval(alg::NormalCholeskyFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) A_ = convert(AbstractMatrix, A) - return ArrayInterface.cholesky_instance(Symmetric((A)' * A), alg.pivot) + return ArrayInterface.cholesky_instance(Symmetric(Matrix{eltype(A)}(undef,0,0)), alg.pivot) end function init_cacheval(alg::NormalCholeskyFactorization, diff --git a/test/default_algs.jl b/test/default_algs.jl index 431e2d155..11d536136 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -45,24 +45,25 @@ solve(prob) prob = LinearProblem(sprand(11000, 11000, 0.5), zeros(11000)) solve(prob) -@static if VERSION >= v"v1.9-" - # Test inference - A = rand(4, 4) - b = rand(4) - prob = LinearProblem(A, b) - VERSION ≥ v"1.10-" && JET.@test_opt init(prob, nothing) - JET.@test_opt solve(prob, LUFactorization()) - JET.@test_opt solve(prob, GenericLUFactorization()) - @test_skip JET.@test_opt solve(prob, QRFactorization()) - JET.@test_opt solve(prob, DiagonalFactorization()) - #JET.@test_opt solve(prob, SVDFactorization()) - #JET.@test_opt solve(prob, KrylovJL_GMRES()) +# Test inference +A = rand(4, 4) +b = rand(4) +prob = LinearProblem(A, b) +VERSION ≥ v"1.10-" && JET.@test_opt init(prob, nothing) +JET.@test_opt solve(prob, LUFactorization()) +JET.@test_opt solve(prob, GenericLUFactorization()) +@test_skip JET.@test_opt solve(prob, QRFactorization()) +JET.@test_opt solve(prob, DiagonalFactorization()) +#JET.@test_opt solve(prob, SVDFactorization()) +#JET.@test_opt solve(prob, KrylovJL_GMRES()) - prob = LinearProblem(sparse(A), b) - #JET.@test_opt solve(prob, UMFPACKFactorization()) - #JET.@test_opt solve(prob, KLUFactorization()) - #JET.@test_opt solve(prob, SparspakFactorization()) - #JET.@test_opt solve(prob) - @inferred solve(prob) - @inferred init(prob, nothing) -end +prob = LinearProblem(sparse(A), b) +#JET.@test_opt solve(prob, UMFPACKFactorization()) +#JET.@test_opt solve(prob, KLUFactorization()) +#JET.@test_opt solve(prob, SparspakFactorization()) +#JET.@test_opt solve(prob) +@inferred solve(prob) +@inferred init(prob, nothing) + +prob = LinearProblem(big.(rand(10, 10)), big.(zeros(10))) +solve(prob) From 30ff6eb5e21e14f71830abdc123abbfc2098f39c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 26 Dec 2023 12:49:56 -0500 Subject: [PATCH 11/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 547f611bb..4abc081ea 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.21.2" +version = "2.22.0" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From 40743e101edc7abbaa90d066e5cb82ae7409c969 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 29 Dec 2023 00:40:17 -0500 Subject: [PATCH 12/55] Fix and test solvers for non-square operators Fixes https://github.com/SciML/LinearSolve.jl/issues/414 --- src/LinearSolve.jl | 2 ++ src/default.jl | 16 +++++++++------- src/iterative_wrappers.jl | 18 ++++++++++++++++-- test/default_algs.jl | 39 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 66 insertions(+), 9 deletions(-) diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index fc039d13f..74f0efa7d 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -103,6 +103,8 @@ EnumX.@enumx DefaultAlgorithmChoice begin AppleAccelerateLUFactorization MKLLUFactorization QRFactorizationPivoted + KrylovJL_CRAIGMR + KrylovJL_LSMR end struct DefaultLinearSolver <: SciMLLinearSolveAlgorithm diff --git a/src/default.jl b/src/default.jl index 0d2daf19c..8972b2a38 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, T17, T18, T19} + T13, T14, T15, T16, T17, T18, T19, T20, T21} LUFactorization::T1 QRFactorization::T2 DiagonalFactorization::T3 @@ -20,6 +20,8 @@ mutable struct DefaultLinearSolverInit{T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, AppleAccelerateLUFactorization::T17 MKLLUFactorization::T18 QRFactorizationPivoted::T19 + KrylovJL_CRAIGMR::T20 + KrylovJL_LSMR::T21 end # Legacy fallback @@ -254,11 +256,11 @@ function algchoice_to_alg(alg::Symbol) elseif alg === :AppleAccelerateLUFactorization AppleAccelerateLUFactorization() elseif alg === :QRFactorizationPivoted - @static if VERSION ≥ v"1.7beta" - QRFactorization(ColumnNorm()) - else - QRFactorization(Val(true)) - end + QRFactorization(ColumnNorm()) + elseif alg === :KrylovJL_CRAIGMR + KrylovJL_CRAIGMR() + elseif alg === :KrylovJL_LSMR + KrylovJL_LSMR() else error("Algorithm choice symbol $alg not allowed in the default") end @@ -387,7 +389,7 @@ end quote getproperty(cache.cacheval,$(Meta.quot(alg)))' \ dy end - elseif alg in Symbol.((DefaultAlgorithmChoice.KrylovJL_GMRES,)) + elseif alg in Symbol.((DefaultAlgorithmChoice.KrylovJL_GMRES,DefaultAlgorithmChoice.KrylovJL_LSMR, DefaultAlgorithmChoice.KrylovJL_CRAIGMR)) quote invprob = LinearSolve.LinearProblem(transpose(cache.A), dy) solve(invprob, cache.alg; diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl index 294cfe7f1..319d373e2 100644 --- a/src/iterative_wrappers.jl +++ b/src/iterative_wrappers.jl @@ -243,7 +243,21 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...) itmax = cache.maxiters verbose = cache.verbose ? 1 : 0 - args = (@get_cacheval(cache, :KrylovJL_GMRES), cache.A, cache.b) + cacheval = if cache.alg isa DefaultLinearSolver + if alg.KrylovAlg === Krylov.gmres! + @get_cacheval(cache, :KrylovJL_GMRES) + elseif alg.KrylovAlg === Krylov.craigmr! + @get_cacheval(cache, :KrylovJL_CRAIGMR) + elseif alg.KrylovAlg === Krylov.lsmr! + @get_cacheval(cache, :KrylovJL_LSMR) + else + error("Default linear solver can only be these three choices! Report this bug!") + end + else + cache.cacheval + end + + args = (cacheval, cache.A, cache.b) kwargs = (atol = atol, rtol = rtol, itmax = itmax, verbose = verbose, ldiv = true, history = true, alg.kwargs...) @@ -268,7 +282,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...) end stats = @get_cacheval(cache, :KrylovJL_GMRES).stats - resid = stats.residuals |> last + resid = !isempty(stats.residuals) ? last(stats.residuals) : zero(eltype(stats.residuals)) retcode = if !stats.solved if stats.status == "maximum number of iterations exceeded" diff --git a/test/default_algs.jl b/test/default_algs.jl index 11d536136..5aad1a30c 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -67,3 +67,42 @@ prob = LinearProblem(sparse(A), b) prob = LinearProblem(big.(rand(10, 10)), big.(zeros(10))) solve(prob) + +## Operator defaults +## https://github.com/SciML/LinearSolve.jl/issues/414 + +m, n = 2, 2 +A = rand(m, n) +b = rand(m) +x = rand(n) +f = (du, u, p, t) -> mul!(du, A, u) +fadj = (du, u, p, t) -> mul!(du, A', u) +fo = FunctionOperator(f, x, b; op_adjoint = fadj) +prob = LinearProblem(fo, b) +sol1 = solve(prob) +sol2 = solve(prob, LinearSolve.KrylovJL_GMRES()) +@test sol1.u == sol2.u + +m, n = 3, 2 +A = rand(m, n) +b = rand(m) +x = rand(n) +f = (du, u, p, t) -> mul!(du, A, u) +fadj = (du, u, p, t) -> mul!(du, A', u) +fo = FunctionOperator(f, x, b; op_adjoint = fadj) +prob = LinearProblem(fo, b) +sol1 = solve(prob) +sol2 = solve(prob, LinearSolve.KrylovJL_LSMR()) +@test sol1.u == sol2.u + +m, n = 2, 3 +A = rand(m, n) +b = rand(m) +x = rand(n) +f = (du, u, p, t) -> mul!(du, A, u) +fadj = (du, u, p, t) -> mul!(du, A', u) +fo = FunctionOperator(f, x, b; op_adjoint = fadj) +prob = LinearProblem(fo, b) +sol1 = solve(prob) +sol2 = solve(prob, LinearSolve.KrylovJL_CRAIGMR()) +@test sol1.u == sol2.u From dfa74a22289a944e8e3ec30b680d5529e8286915 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 7 Jan 2024 12:27:38 -0500 Subject: [PATCH 13/55] Disable error throwing in RFLU All are supposed to be run with check=false --- src/factorization.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/factorization.jl b/src/factorization.jl index 7c261ec6b..9d9b64804 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -957,7 +957,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::RFLUFactorization{P, T}; if length(ipiv) != min(size(A)...) ipiv = Vector{LinearAlgebra.BlasInt}(undef, min(size(A)...)) end - fact = RecursiveFactorization.lu!(A, ipiv, Val(P), Val(T)) + fact = RecursiveFactorization.lu!(A, ipiv, Val(P), Val(T), check=false) cache.cacheval = (fact, ipiv) cache.isfresh = false end From 49183f7b68833686eaa7636b13a1bd49e381446b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 7 Jan 2024 12:42:40 -0500 Subject: [PATCH 14/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4abc081ea..cf0dfe020 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.22.0" +version = "2.22.1" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From 25b92bf1124b5718236bec2f02fbde3fff97aa7c Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 8 Jan 2024 05:20:50 +0000 Subject: [PATCH 15/55] Bump crate-ci/typos from 1.16.23 to 1.17.0 Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.16.23 to 1.17.0. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.16.23...v1.17.0) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 74af4eff7..0decc88b4 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.16.23 \ No newline at end of file + uses: crate-ci/typos@v1.17.0 \ No newline at end of file From 999904484df26d92281b23d2faae99744f4fede5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 8 Jan 2024 21:29:21 -0500 Subject: [PATCH 16/55] Fix 32-bit OS compatibility with Krylov wrappers Fixes https://github.com/SciML/DelayDiffEq.jl/actions/runs/7455784891/job/20285397505?pr=279 downstream --- src/iterative_wrappers.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl index 319d373e2..dbdef5b9b 100644 --- a/src/iterative_wrappers.jl +++ b/src/iterative_wrappers.jl @@ -181,7 +181,7 @@ function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters::Int, abstol, re alg.KrylovAlg === Krylov.gpmr! || alg.KrylovAlg === Krylov.fom!) if A isa SparseMatrixCSC - KS(SparseMatrixCSC(0, 0, [1], Int64[], eltype(A)[]), eltype(b)[], 1) + KS(SparseMatrixCSC(0, 0, [1], Int[], eltype(A)[]), eltype(b)[], 1) elseif A isa Matrix KS(Matrix{eltype(A)}(undef, 0, 0), eltype(b)[], 1) else @@ -189,7 +189,7 @@ function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters::Int, abstol, re end else if A isa SparseMatrixCSC - KS(SparseMatrixCSC(0, 0, [1], Int64[], eltype(A)[]), eltype(b)[]) + KS(SparseMatrixCSC(0, 0, [1], Int[], eltype(A)[]), eltype(b)[]) elseif A isa Matrix KS(Matrix{eltype(A)}(undef, 0, 0), eltype(b)[]) else From a78fb8b5401ba612b724cedbc395e768a7e0b171 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 22 Jan 2024 05:35:59 +0000 Subject: [PATCH 17/55] Bump actions/cache from 3 to 4 Bumps [actions/cache](https://github.com/actions/cache) from 3 to 4. - [Release notes](https://github.com/actions/cache/releases) - [Changelog](https://github.com/actions/cache/blob/main/RELEASES.md) - [Commits](https://github.com/actions/cache/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/cache dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 794935d72..82df8524f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -32,7 +32,7 @@ jobs: - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} - - uses: actions/cache@v3 + - uses: actions/cache@v4 env: cache-name: cache-artifacts with: From 5f5abb5bdccb91c028a79ad3bec20e42cdb0d6c4 Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 22 Jan 2024 22:01:16 +0100 Subject: [PATCH 18/55] Fix typo. --- docs/src/basics/Preconditioners.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/basics/Preconditioners.md b/docs/src/basics/Preconditioners.md index 926b5dba7..1325694ad 100644 --- a/docs/src/basics/Preconditioners.md +++ b/docs/src/basics/Preconditioners.md @@ -83,13 +83,13 @@ The following preconditioners match the interface of LinearSolve.jl. - [Preconditioners.CholeskyPreconditioner(A, i)](https://github.com/JuliaLinearAlgebra/Preconditioners.jl): An incomplete Cholesky preconditioner with cut-off level `i`. Requires `A` as a `AbstractMatrix` and positive semi-definite. - - [AlgebraicMultiGrid](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl): + - [AlgebraicMultigrid](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl): Implementations of the algebraic multigrid method. Must be converted to a - preconditioner via `AlgebraicMultiGrid.aspreconditioner(AlgebraicMultiGrid.precmethod(A))`. + preconditioner via `AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.precmethod(A))`. Requires `A` as a `AbstractMatrix`. Provides the following methods: - + `AlgebraicMultiGrid.ruge_stuben(A)` - + `AlgebraicMultiGrid.smoothed_aggregation(A)` + + `AlgebraicMultigrid.ruge_stuben(A)` + + `AlgebraicMultigrid.smoothed_aggregation(A)` - [PyAMG](https://github.com/cortner/PyAMG.jl): Implementations of the algebraic multigrid method. Must be converted to a preconditioner via `PyAMG.aspreconditioner(PyAMG.precmethod(A))`. From eb659438187a23a8b3f5bded25fd683230d0552d Mon Sep 17 00:00:00 2001 From: termi-official Date: Mon, 22 Jan 2024 22:14:08 +0100 Subject: [PATCH 19/55] add KrylovPreconditioners.jl --- docs/src/basics/Preconditioners.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/src/basics/Preconditioners.md b/docs/src/basics/Preconditioners.md index 1325694ad..05faedb06 100644 --- a/docs/src/basics/Preconditioners.md +++ b/docs/src/basics/Preconditioners.md @@ -111,3 +111,9 @@ The following preconditioners match the interface of LinearSolve.jl. preconditioners which supports distributed computing via MPI. These can be written using the LinearSolve.jl interface choosing algorithms like `HYPRE.ILU` and `HYPRE.BoomerAMG`. + - [KrylovPreconditioners.jl](https://github.com/JuliaSmoothOptimizers/KrylovPreconditioners.jl/): Provides GPU-ready + preconditioners via KernelAbstractions.jl. At the time of writing the package provides the following methods: + + + Incomplete Cholesky decomposition `KrylovPreconditioners.kp_ic0(A)` + + Incomplete LU decomposition `KrylovPreconditioners.kp_ilu0(A)` + + Block Jacobi `KrylovPreconditioners.BlockJacobiPreconditioner(A, nblocks, device)` From b98cf5290447c89aa60deb1f1de5b96f33a84268 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 5 Feb 2024 05:49:07 +0000 Subject: [PATCH 20/55] Bump crate-ci/typos from 1.17.0 to 1.18.0 Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.17.0 to 1.18.0. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.17.0...v1.18.0) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 0decc88b4..9246edd2a 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.17.0 \ No newline at end of file + uses: crate-ci/typos@v1.18.0 \ No newline at end of file From 9328fa4dfe730e5621466807db64a0c473b889a7 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Tue, 6 Feb 2024 04:39:11 +0100 Subject: [PATCH 21/55] [skip ci] Update dependabot.yml --- .github/dependabot.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/dependabot.yml b/.github/dependabot.yml index 1e8a051e2..ec3b005a0 100644 --- a/.github/dependabot.yml +++ b/.github/dependabot.yml @@ -7,4 +7,4 @@ updates: interval: "weekly" ignore: - dependency-name: "crate-ci/typos" - update-types: ["version-update:semver-patch"] + update-types: ["version-update:semver-patch", "version-update:semver-minor"] From 09d7f844f3401c62c4e9b34572ebd8a4f1bf5487 Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Sat, 10 Feb 2024 00:10:49 +0000 Subject: [PATCH 22/55] CompatHelper: bump compat for KLU to 0.5, (keep existing compat) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index cf0dfe020..3bd31cc8a 100644 --- a/Project.toml +++ b/Project.toml @@ -78,7 +78,7 @@ HYPRE = "1.4.0" InteractiveUtils = "1.9" IterativeSolvers = "0.9.3" JET = "0.8" -KLU = "0.3.0, 0.4" +KLU = "0.3.0, 0.4, 0.5" KernelAbstractions = "0.9" Krylov = "0.9" KrylovKit = "0.6" From 5525694f9a98feb781929f1ba9390813f25736b8 Mon Sep 17 00:00:00 2001 From: Keno Fischer Date: Sat, 10 Feb 2024 21:07:37 +0000 Subject: [PATCH 23/55] Don't use type-printing from inside a generated function Base's type printing code uses the REPL's currently active module to determine how to print types. Unfortunately, this means it is illegal to use this in generated functions. Switch to grabbing the symbol directly from the typename, since this now errors on julia master: ``` Got exception outside of a @test LoadError: MethodError: no method matching active_module() You may have intended to import Base.active_module The applicable method may be too new: running in world age 26446, while current world is 26576. Closest candidates are: active_module(::REPL.LineEditREPL) (method too new to be called from this world context.) @ REPL ~/julia/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:590 active_module() (method too new to be called from this world context.) @ REPL ~/julia/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:586 active_module(::REPL.REPLDisplay) (method too new to be called from this world context.) @ REPL ~/julia/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:592 ... Stacktrace: [1] #invokelatest#2 @ ./essentials.jl:1025 [inlined] [2] invokelatest @ ./essentials.jl:1022 [inlined] [3] active_module @ ./show.jl:512 [inlined] [4] show_type_name(io::IOBuffer, tn::Core.TypeName) @ Base ./show.jl:1051 [5] _show_type(io::IOBuffer, x::Type) @ Base ./show.jl:966 [6] show(io::IOBuffer, x::Type) @ Base ./show.jl:963 [7] print(io::IOBuffer, x::Type) @ Base ./strings/io.jl:35 [8] print_to_string(xs::Type) @ Base ./strings/io.jl:148 [9] string @ ./strings/io.jl:189 [inlined] [10] defaultalg_symbol(::Type{LinearSolve.LUFactorization{LinearAlgebra.RowMaximum}}) @ LinearSolve ~/.julia/packages/LinearSolve/4aydZ/src/default.jl:323 ``` --- src/default.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/default.jl b/src/default.jl index 8972b2a38..31c24a020 100644 --- a/src/default.jl +++ b/src/default.jl @@ -320,7 +320,7 @@ cache.cacheval = NamedTuple(LUFactorization = cache of LUFactorization, ...) end function defaultalg_symbol(::Type{T}) where {T} - Symbol(split(string(SciMLBase.parameterless_type(T)), ".")[end]) + Base.typename(SciMLBase.parameterless_type(T)).name end defaultalg_symbol(::Type{<:GenericFactorization{typeof(ldlt!)}}) = :LDLtFactorization From d20c688547131342bcbd35665e6766dc6b64855c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:01:20 -0500 Subject: [PATCH 24/55] Fix downgrade CI --- Project.toml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 3bd31cc8a..104758750 100644 --- a/Project.toml +++ b/Project.toml @@ -75,15 +75,15 @@ FiniteDiff = "2.18" ForwardDiff = "0.10.13" GPUArraysCore = "0.1" HYPRE = "1.4.0" -InteractiveUtils = "1.9" +InteractiveUtils = "1.10" IterativeSolvers = "0.9.3" JET = "0.8" KLU = "0.3.0, 0.4, 0.5" KernelAbstractions = "0.9" Krylov = "0.9" KrylovKit = "0.6" -Libdl = "1.9" -LinearAlgebra = "1.9" +Libdl = "1.10" +LinearAlgebra = "1.10" MPI = "0.20" Metal = "0.5" MultiFloats = "1" @@ -99,14 +99,14 @@ SafeTestsets = "0.1" SciMLBase = "2" SciMLOperators = "0.3" Setfield = "1" -SparseArrays = "1.9" +SparseArrays = "1.10" Sparspak = "0.3.6" StableRNGs = "1" StaticArrays = "1.5" StaticArraysCore = "1.2" Test = "1" UnPack = "1" -julia = "1.9" +julia = "1.10" [extras] AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a" From 5f91cd4144be59de67150ff0087c05a21931c21c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:12:01 -0500 Subject: [PATCH 25/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 104758750..793b5b8cd 100644 --- a/Project.toml +++ b/Project.toml @@ -96,7 +96,7 @@ RecursiveArrayTools = "2.38, 3" RecursiveFactorization = "0.2.14" Reexport = "1" SafeTestsets = "0.1" -SciMLBase = "2" +SciMLBase = "2.24.0" SciMLOperators = "0.3" Setfield = "1" SparseArrays = "1.10" From 7388015609bf3c36fc1131d98a99493fd7156bc5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:14:48 -0500 Subject: [PATCH 26/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 793b5b8cd..230d648e9 100644 --- a/Project.toml +++ b/Project.toml @@ -97,7 +97,7 @@ RecursiveFactorization = "0.2.14" Reexport = "1" SafeTestsets = "0.1" SciMLBase = "2.24.0" -SciMLOperators = "0.3" +SciMLOperators = "0.3.7" Setfield = "1" SparseArrays = "1.10" Sparspak = "0.3.6" From 510737bda9695d93fd060938a712d51408828404 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:14:59 -0500 Subject: [PATCH 27/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 230d648e9..c9ee0d455 100644 --- a/Project.toml +++ b/Project.toml @@ -96,7 +96,7 @@ RecursiveArrayTools = "2.38, 3" RecursiveFactorization = "0.2.14" Reexport = "1" SafeTestsets = "0.1" -SciMLBase = "2.24.0" +SciMLBase = "2.23.0" SciMLOperators = "0.3.7" Setfield = "1" SparseArrays = "1.10" From e6bb2df6774d1d46d08951b673acf222c7553df5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:20:11 -0500 Subject: [PATCH 28/55] Update Project.toml --- Project.toml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index c9ee0d455..959594eb1 100644 --- a/Project.toml +++ b/Project.toml @@ -60,11 +60,11 @@ LinearSolveRecursiveArrayToolsExt = "RecursiveArrayTools" [compat] AllocCheck = "0.1" Aqua = "0.8" -ArrayInterface = "7.5" -BandedMatrices = "1" +ArrayInterface = "7.7" +BandedMatrices = "1.5" BlockDiagonals = "0.1.42" CUDA = "5" -ConcreteStructs = "0.2" +ConcreteStructs = "0.2.3" DocStringExtensions = "0.9" EnumX = "1" Enzyme = "0.11.10" From 7842bb2fdef7cc7dbd31078df71b690b12fe661f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:21:52 -0500 Subject: [PATCH 29/55] Update Project.toml --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 959594eb1..8c012da93 100644 --- a/Project.toml +++ b/Project.toml @@ -78,7 +78,7 @@ HYPRE = "1.4.0" InteractiveUtils = "1.10" IterativeSolvers = "0.9.3" JET = "0.8" -KLU = "0.3.0, 0.4, 0.5" +KLU = "0.5" KernelAbstractions = "0.9" Krylov = "0.9" KrylovKit = "0.6" @@ -92,7 +92,7 @@ Pkg = "1" PrecompileTools = "1.1" Preferences = "1.4" Random = "1" -RecursiveArrayTools = "2.38, 3" +RecursiveArrayTools = "3.8" RecursiveFactorization = "0.2.14" Reexport = "1" SafeTestsets = "0.1" From b1235ac9a49a10eefc395da7b04cc9986b073411 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:23:41 -0500 Subject: [PATCH 30/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8c012da93..e28c229e1 100644 --- a/Project.toml +++ b/Project.toml @@ -103,7 +103,7 @@ SparseArrays = "1.10" Sparspak = "0.3.6" StableRNGs = "1" StaticArrays = "1.5" -StaticArraysCore = "1.2" +StaticArraysCore = "1.4.2" Test = "1" UnPack = "1" julia = "1.10" From 98ed289fb76cf8f8b6c548baf59002cefe448734 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:40:18 -0500 Subject: [PATCH 31/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e28c229e1..8e3f7c41b 100644 --- a/Project.toml +++ b/Project.toml @@ -89,7 +89,7 @@ Metal = "0.5" MultiFloats = "1" Pardiso = "0.5" Pkg = "1" -PrecompileTools = "1.1" +PrecompileTools = "1.2" Preferences = "1.4" Random = "1" RecursiveArrayTools = "3.8" From 4a2ebb9142a84d9dcf0688958e4c3501fd38e77a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:47:24 -0500 Subject: [PATCH 32/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8e3f7c41b..ce82b7fc8 100644 --- a/Project.toml +++ b/Project.toml @@ -73,7 +73,7 @@ FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" FiniteDiff = "2.18" ForwardDiff = "0.10.13" -GPUArraysCore = "0.1" +GPUArraysCore = "0.1.6" HYPRE = "1.4.0" InteractiveUtils = "1.10" IterativeSolvers = "0.9.3" From 06bdd371db93f96b742b860c19de721bf00724ca Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:50:16 -0500 Subject: [PATCH 33/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ce82b7fc8..394fe324c 100644 --- a/Project.toml +++ b/Project.toml @@ -72,7 +72,7 @@ EnzymeCore = "0.6.2" FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" FiniteDiff = "2.18" -ForwardDiff = "0.10.13" +ForwardDiff = "0.10.36" GPUArraysCore = "0.1.6" HYPRE = "1.4.0" InteractiveUtils = "1.10" From 399705be52a1aec96ae729173dcccc268ed438e0 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:54:15 -0500 Subject: [PATCH 34/55] Update Project.toml --- Project.toml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 394fe324c..3f3b27db7 100644 --- a/Project.toml +++ b/Project.toml @@ -65,21 +65,21 @@ BandedMatrices = "1.5" BlockDiagonals = "0.1.42" CUDA = "5" ConcreteStructs = "0.2.3" -DocStringExtensions = "0.9" -EnumX = "1" +DocStringExtensions = "0.9.3" +EnumX = "1.0.4" Enzyme = "0.11.10" EnzymeCore = "0.6.2" FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" -FiniteDiff = "2.18" +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" +JET = "0.8.28" KLU = "0.5" -KernelAbstractions = "0.9" +KernelAbstractions = "0.9.16" Krylov = "0.9" KrylovKit = "0.6" Libdl = "1.10" From f679966f5b1ab086f1e4cd0eade6fb07f1fccfe6 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:56:05 -0500 Subject: [PATCH 35/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3f3b27db7..6648fb41c 100644 --- a/Project.toml +++ b/Project.toml @@ -67,7 +67,7 @@ CUDA = "5" ConcreteStructs = "0.2.3" DocStringExtensions = "0.9.3" EnumX = "1.0.4" -Enzyme = "0.11.10" +Enzyme = "0.11.14" EnzymeCore = "0.6.2" FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" From ac69ea31ca925a453a3417a2a05f4bad1923890e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 21:58:12 -0500 Subject: [PATCH 36/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6648fb41c..915f339b3 100644 --- a/Project.toml +++ b/Project.toml @@ -68,7 +68,7 @@ ConcreteStructs = "0.2.3" DocStringExtensions = "0.9.3" EnumX = "1.0.4" Enzyme = "0.11.14" -EnzymeCore = "0.6.2" +EnzymeCore = "0.6.5" FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" FiniteDiff = "2.22" From fa5278e922bd80ed94c54a7d7b332e1629c07bcd Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 10 Feb 2024 22:37:40 -0500 Subject: [PATCH 37/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 915f339b3..9b79023e2 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.22.1" +version = "2.23.0" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From 908605158423b1b4939b32dee351ef6e0a68c35f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 11 Feb 2024 15:26:11 -0500 Subject: [PATCH 38/55] Avoid constructing any Krylov cache with default dense --- Project.toml | 2 +- src/default.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 9b79023e2..8efa1c37d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.23.0" +version = "2.23.1" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" diff --git a/src/default.jl b/src/default.jl index 31c24a020..48e8a9b1f 100644 --- a/src/default.jl +++ b/src/default.jl @@ -121,7 +121,7 @@ function defaultalg(A::Nothing, b::GPUArraysCore.AnyGPUArray, assump::OperatorAs end # Ambiguity handling -function defaultalg(A::GPUArraysCore.AnyGPUArray, b::GPUArraysCore.AbstractGPUArray, +function defaultalg(A::GPUArraysCore.AnyGPUArray, b::GPUArraysCore.AnyGPUArray, assump::OperatorAssumptions{Bool}) if assump.condition === OperatorCondition.IllConditioned || !assump.issq DefaultLinearSolver(DefaultAlgorithmChoice.QRFactorization) @@ -152,7 +152,7 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) # Special case on Arrays: avoid BLAS for RecursiveFactorization.jl when # it makes sense according to the benchmarks, which is dependent on # whether MKL or OpenBLAS is being used - if (A === nothing && !(b isa GPUArraysCore.AbstractGPUArray)) || A isa Matrix + if (A === nothing && !(b isa GPUArraysCore.AnyGPUArray)) || A isa Matrix if (A === nothing || eltype(A) <: BLASELTYPES) && ArrayInterface.can_setindex(b) && @@ -296,7 +296,7 @@ cache.cacheval = NamedTuple(LUFactorization = cache of LUFactorization, ...) abstol, reltol, verbose::Bool, assump::OperatorAssumptions) caches = map(first.(EnumX.symbol_map(DefaultAlgorithmChoice.T))) do alg - if alg === :KrylovJL_GMRES + if alg === :KrylovJL_GMRES || alg === :KrylovJL_CRAIGMR || alg === :KrylovJL_LSMR quote if A isa Matrix || A isa SparseMatrixCSC nothing From 91fe4d42e807aecd62c149f4f87cfd96bb215a93 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 11 Feb 2024 15:30:14 -0500 Subject: [PATCH 39/55] fix spell check --- test/default_algs.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/default_algs.jl b/test/default_algs.jl index 5aad1a30c..2dedc13a5 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -77,8 +77,8 @@ b = rand(m) x = rand(n) f = (du, u, p, t) -> mul!(du, A, u) fadj = (du, u, p, t) -> mul!(du, A', u) -fo = FunctionOperator(f, x, b; op_adjoint = fadj) -prob = LinearProblem(fo, b) +funcop = FunctionOperator(f, x, b; op_adjoint = fadj) +prob = LinearProblem(funcop, b) sol1 = solve(prob) sol2 = solve(prob, LinearSolve.KrylovJL_GMRES()) @test sol1.u == sol2.u @@ -89,8 +89,8 @@ b = rand(m) x = rand(n) f = (du, u, p, t) -> mul!(du, A, u) fadj = (du, u, p, t) -> mul!(du, A', u) -fo = FunctionOperator(f, x, b; op_adjoint = fadj) -prob = LinearProblem(fo, b) +funcop = FunctionOperator(f, x, b; op_adjoint = fadj) +prob = LinearProblem(funcop, b) sol1 = solve(prob) sol2 = solve(prob, LinearSolve.KrylovJL_LSMR()) @test sol1.u == sol2.u @@ -101,8 +101,8 @@ b = rand(m) x = rand(n) f = (du, u, p, t) -> mul!(du, A, u) fadj = (du, u, p, t) -> mul!(du, A', u) -fo = FunctionOperator(f, x, b; op_adjoint = fadj) -prob = LinearProblem(fo, b) +funcop = FunctionOperator(f, x, b; op_adjoint = fadj) +prob = LinearProblem(funcop, b) sol1 = solve(prob) sol2 = solve(prob, LinearSolve.KrylovJL_CRAIGMR()) @test sol1.u == sol2.u From 8becdf48e44f43aa373bf065aacc9ea67c4978a7 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Sun, 11 Feb 2024 20:36:28 -0500 Subject: [PATCH 40/55] Check appleaccelerate is available in __init__ --- Project.toml | 2 +- src/LinearSolve.jl | 3 +++ src/appleaccelerate.jl | 20 ++++++++++++-------- src/init.jl | 2 ++ 4 files changed, 18 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 8efa1c37d..4613921b9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.23.1" +version = "2.23.2" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 74f0efa7d..7ffaadf9c 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -182,6 +182,9 @@ 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) diff --git a/src/appleaccelerate.jl b/src/appleaccelerate.jl index b6b089f0f..44bb07dcb 100644 --- a/src/appleaccelerate.jl +++ b/src/appleaccelerate.jl @@ -14,16 +14,20 @@ to avoid allocations and does not require libblastrampoline. """ struct AppleAccelerateLUFactorization <: AbstractFactorization end -function appleaccelerate_isavailable() - libacc_hdl = Libdl.dlopen_e(libacc) - if libacc_hdl == C_NULL - return false - end +@static if !Sys.isapple() + __appleaccelerate_isavailable() = false +else + function __appleaccelerate_isavailable() + libacc_hdl = Libdl.dlopen_e(libacc) + if libacc_hdl == C_NULL + return false + end - if dlsym_e(libacc_hdl, "dgetrf_") == C_NULL - return false + if dlsym_e(libacc_hdl, "dgetrf_") == C_NULL + return false + end + return true end - return true end function aa_getrf!(A::AbstractMatrix{<:ComplexF64}; diff --git a/src/init.jl b/src/init.jl index e10960f1d..12df4865d 100644 --- a/src/init.jl +++ b/src/init.jl @@ -1,3 +1,5 @@ function __init__() IS_OPENBLAS[] = occursin("openblas", BLAS.get_config().loaded_libs[1].libname) + + HAS_APPLE_ACCELERATE[] = __appleaccelerate_isavailable() end From 8e5e41e2440d6902ad399579e0d642788a27e053 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 13 Feb 2024 12:11:46 -0500 Subject: [PATCH 41/55] Preallocate more caches --- Project.toml | 2 +- src/default.jl | 3 ++- src/factorization.jl | 25 ++++++++++++++++++++----- 3 files changed, 23 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 4613921b9..f5831acef 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.23.2" +version = "2.23.3" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" diff --git a/src/default.jl b/src/default.jl index 48e8a9b1f..13697bd14 100644 --- a/src/default.jl +++ b/src/default.jl @@ -159,7 +159,8 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) (__conditioning(assump) === OperatorCondition.IllConditioned || __conditioning(assump) === OperatorCondition.WellConditioned) if length(b) <= 10 - DefaultAlgorithmChoice.GenericLUFactorization + # DefaultAlgorithmChoice.GenericLUFactorization + DefaultAlgorithmChoice.RFLUFactorization elseif appleaccelerate_isavailable() DefaultAlgorithmChoice.AppleAccelerateLUFactorization elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500) || diff --git a/src/factorization.jl b/src/factorization.jl index 9d9b64804..60fe39c66 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -173,18 +173,24 @@ function init_cacheval(alg::QRFactorization, A, b, u, Pl, Pr, ArrayInterface.qr_instance(convert(AbstractMatrix, A), alg.pivot) end +const PREALLOCATED_QR_ColumnNorm = ArrayInterface.qr_instance(rand(1, 1), ColumnNorm()) + +function init_cacheval(alg::QRFactorization{ColumnNorm}, A::Matrix{Float64}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + return PREALLOCATED_QR_ColumnNorm +end + function init_cacheval(alg::QRFactorization, A::Union{<:Adjoint, <:Transpose}, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) A isa GPUArraysCore.AnyGPUArray && return qr(A) return qr(A, alg.pivot) end -const PREALLOCATED_QR = ArrayInterface.qr_instance(rand(1, 1)) +const PREALLOCATED_QR_NoPivot = ArrayInterface.qr_instance(rand(1, 1)) function init_cacheval(alg::QRFactorization{NoPivot}, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) - PREALLOCATED_QR + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + return PREALLOCATED_QR_NoPivot end function init_cacheval(alg::QRFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, @@ -1010,7 +1016,16 @@ function init_cacheval(alg::NormalCholeskyFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) A_ = convert(AbstractMatrix, A) - return ArrayInterface.cholesky_instance(Symmetric(Matrix{eltype(A)}(undef,0,0)), alg.pivot) + MType = ArrayInterface.parameterless_type(A_) + return ArrayInterface.cholesky_instance(Symmetric(MType{eltype(A), 2}(undef,0,0)), alg.pivot) +end + +const PREALLOCATED_NORMALCHOLESKY = ArrayInterface.cholesky_instance( + Symmetric(rand(1, 1)), NoPivot()) + +function init_cacheval(alg::NormalCholeskyFactorization, A::Matrix, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + return PREALLOCATED_NORMALCHOLESKY end function init_cacheval(alg::NormalCholeskyFactorization, From 6fd247c4bd32cd33446adb149b0d93f200391412 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 13 Feb 2024 12:16:44 -0500 Subject: [PATCH 42/55] Restrict to Float64 --- src/default.jl | 3 +-- src/factorization.jl | 6 +++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/default.jl b/src/default.jl index 13697bd14..48e8a9b1f 100644 --- a/src/default.jl +++ b/src/default.jl @@ -159,8 +159,7 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) (__conditioning(assump) === OperatorCondition.IllConditioned || __conditioning(assump) === OperatorCondition.WellConditioned) if length(b) <= 10 - # DefaultAlgorithmChoice.GenericLUFactorization - DefaultAlgorithmChoice.RFLUFactorization + DefaultAlgorithmChoice.GenericLUFactorization elseif appleaccelerate_isavailable() DefaultAlgorithmChoice.AppleAccelerateLUFactorization elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500) || diff --git a/src/factorization.jl b/src/factorization.jl index 60fe39c66..51431647b 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -1020,12 +1020,12 @@ function init_cacheval(alg::NormalCholeskyFactorization, A, b, u, Pl, Pr, return ArrayInterface.cholesky_instance(Symmetric(MType{eltype(A), 2}(undef,0,0)), alg.pivot) end -const PREALLOCATED_NORMALCHOLESKY = ArrayInterface.cholesky_instance( +const PREALLOCATED_NORMALCHOLESKY_SYMMETRIC = ArrayInterface.cholesky_instance( Symmetric(rand(1, 1)), NoPivot()) -function init_cacheval(alg::NormalCholeskyFactorization, A::Matrix, b, u, Pl, Pr, +function init_cacheval(alg::NormalCholeskyFactorization, A::Matrix{Float64}, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) - return PREALLOCATED_NORMALCHOLESKY + return PREALLOCATED_NORMALCHOLESKY_SYMMETRIC end function init_cacheval(alg::NormalCholeskyFactorization, From b6bd5f581593785b2eae203e3784e047d95f8c07 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 13 Feb 2024 15:14:56 -0500 Subject: [PATCH 43/55] EnumX --> Symbol is slow --- src/default.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/default.jl b/src/default.jl index 48e8a9b1f..d3d06f1a0 100644 --- a/src/default.jl +++ b/src/default.jl @@ -218,6 +218,7 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) DefaultLinearSolver(alg) end +@inline algchoice_to_alg(::Val{alg}) where {alg} = algchoice_to_alg(alg) function algchoice_to_alg(alg::Symbol) if alg === :SVDFactorization SVDFactorization(false, LinearAlgebra.QRIteration()) @@ -345,11 +346,12 @@ end retcode = sol.retcode, iters = sol.iters, stats = sol.stats) end + alg_enum = getproperty(LinearSolve.DefaultAlgorithmChoice, alg) ex = if ex == :() - Expr(:elseif, :(Symbol(alg.alg) === $(Meta.quot(alg))), newex, + Expr(:elseif, :(alg.alg == $(alg_enum)), newex, :(error("Algorithm Choice not Allowed"))) else - Expr(:elseif, :(Symbol(alg.alg) === $(Meta.quot(alg))), newex, ex) + Expr(:elseif, :(alg.alg == $(alg_enum)), newex, ex) end end ex = Expr(:if, ex.args...) From f49b8d865ac883ec03963eb4cab701e38ff50978 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 13 Feb 2024 15:26:56 -0500 Subject: [PATCH 44/55] Remove more uses of Symbol --- src/common.jl | 3 ++- src/default.jl | 18 +++++++++++++++++- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/common.jl b/src/common.jl index b206598d5..1f4d50e2e 100644 --- a/src/common.jl +++ b/src/common.jl @@ -90,7 +90,8 @@ function Base.setproperty!(cache::LinearCache, name::Symbol, x) update_cacheval!(cache, :b, x) elseif name === :cacheval && cache.alg isa DefaultLinearSolver @assert cache.cacheval isa DefaultLinearSolverInit - return setfield!(cache.cacheval, Symbol(cache.alg.alg), x) + return __setfield!(cache.cacheval, cache.alg, x) + # return setfield!(cache.cacheval, Symbol(cache.alg.alg), x) end setfield!(cache, name, x) end diff --git a/src/default.jl b/src/default.jl index d3d06f1a0..f3e522c00 100644 --- a/src/default.jl +++ b/src/default.jl @@ -24,6 +24,23 @@ mutable struct DefaultLinearSolverInit{T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, KrylovJL_LSMR::T21 end +@generated function __setfield!(cache::DefaultLinearSolverInit, alg::DefaultLinearSolver, v) + ex = :() + for alg in first.(EnumX.symbol_map(DefaultAlgorithmChoice.T)) + newex = quote + setfield!(cache, $(Meta.quot(alg)), v) + end + alg_enum = getproperty(LinearSolve.DefaultAlgorithmChoice, alg) + ex = if ex == :() + Expr(:elseif, :(alg.alg == $(alg_enum)), newex, + :(error("Algorithm Choice not Allowed"))) + else + Expr(:elseif, :(alg.alg == $(alg_enum)), newex, ex) + end + end + ex = Expr(:if, ex.args...) +end + # Legacy fallback # For SciML algorithms already using `defaultalg`, all assume square matrix. defaultalg(A, b) = defaultalg(A, b, OperatorAssumptions(true)) @@ -218,7 +235,6 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) DefaultLinearSolver(alg) end -@inline algchoice_to_alg(::Val{alg}) where {alg} = algchoice_to_alg(alg) function algchoice_to_alg(alg::Symbol) if alg === :SVDFactorization SVDFactorization(false, LinearAlgebra.QRIteration()) From 085ba9de147a451ca0942cddf99f659937fa1369 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Tue, 13 Feb 2024 16:09:41 -0500 Subject: [PATCH 45/55] Fix tests --- src/default.jl | 2 +- src/factorization.jl | 3 +-- test/default_algs.jl | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/default.jl b/src/default.jl index f3e522c00..9a76403bb 100644 --- a/src/default.jl +++ b/src/default.jl @@ -176,7 +176,7 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) (__conditioning(assump) === OperatorCondition.IllConditioned || __conditioning(assump) === OperatorCondition.WellConditioned) if length(b) <= 10 - DefaultAlgorithmChoice.GenericLUFactorization + DefaultAlgorithmChoice.RFLUFactorization elseif appleaccelerate_isavailable() DefaultAlgorithmChoice.AppleAccelerateLUFactorization elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500) || diff --git a/src/factorization.jl b/src/factorization.jl index 51431647b..1f244b3c4 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -1016,8 +1016,7 @@ function init_cacheval(alg::NormalCholeskyFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) A_ = convert(AbstractMatrix, A) - MType = ArrayInterface.parameterless_type(A_) - return ArrayInterface.cholesky_instance(Symmetric(MType{eltype(A), 2}(undef,0,0)), alg.pivot) + return ArrayInterface.cholesky_instance(Symmetric(Matrix{eltype(A)}(undef,0,0)), alg.pivot) end const PREALLOCATED_NORMALCHOLESKY_SYMMETRIC = ArrayInterface.cholesky_instance( diff --git a/test/default_algs.jl b/test/default_algs.jl index 2dedc13a5..2ffc03b65 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -1,6 +1,6 @@ using LinearSolve, LinearAlgebra, SparseArrays, Test, JET @test LinearSolve.defaultalg(nothing, zeros(3)).alg === - LinearSolve.DefaultAlgorithmChoice.GenericLUFactorization + LinearSolve.DefaultAlgorithmChoice.RFLUFactorization prob = LinearProblem(rand(3, 3), rand(3)) solve(prob) From 9591654598f1bf9d5ff13f0925b78ab7d4007d4b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 13 Feb 2024 17:12:18 -0500 Subject: [PATCH 46/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f5831acef..b2a5b60f7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.23.3" +version = "2.23.4" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From be1a3211462c85ca695042b749d62d8b32d4c0e7 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 13 Feb 2024 17:19:18 -0500 Subject: [PATCH 47/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index b2a5b60f7..f5831acef 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.23.4" +version = "2.23.3" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From a842602c4a138ee33fe460014f5cce30a884f367 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 15 Feb 2024 11:12:30 -0500 Subject: [PATCH 48/55] Add downstream testing for NonlinearSolve --- .github/workflows/Downstream.yml | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index c72a26c11..620c447cd 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -20,7 +20,7 @@ jobs: - {user: SciML, repo: ModelingToolkit.jl, group: All} - {user: SciML, repo: SciMLSensitivity.jl, group: Core1} - {user: SciML, repo: BoundaryValueDiffEq.jl, group: All} - + - {user: SciML, repo: NonlinearSolve.jl, group: All} steps: - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 @@ -50,6 +50,9 @@ jobs: @info "Not compatible with this release. No problem." exception=err exit(0) # Exit immediately, as a success end + env: + RETESTITEMS_NWORKERS: 4 + RETESTITEMS_NWORKER_THREADS: 2 - uses: julia-actions/julia-processcoverage@v1 - uses: codecov/codecov-action@v3 with: From b9f4c73ff89b8e086b90f1200ab16b1c112daca8 Mon Sep 17 00:00:00 2001 From: Avik Pal Date: Thu, 15 Feb 2024 11:39:10 -0500 Subject: [PATCH 49/55] Bump enzyme compat --- .github/workflows/Downstream.yml | 1 + Project.toml | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 620c447cd..1820e609a 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -12,6 +12,7 @@ jobs: env: GROUP: ${{ matrix.package.group }} strategy: + fail-fast: false matrix: julia-version: [1] os: [ubuntu-latest] diff --git a/Project.toml b/Project.toml index f5831acef..2b3422754 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.23.3" +version = "2.23.4" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" @@ -67,7 +67,7 @@ CUDA = "5" ConcreteStructs = "0.2.3" DocStringExtensions = "0.9.3" EnumX = "1.0.4" -Enzyme = "0.11.14" +Enzyme = "0.11.15" EnzymeCore = "0.6.5" FastAlmostBandedMatrices = "0.1" FastLapackInterface = "2" From a50ee388c2aff8fe0d6e29e5131ad7e13f7b9df3 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Sun, 18 Feb 2024 13:25:42 +0100 Subject: [PATCH 50/55] [skip ci] update downgrade CI repo --- .github/workflows/Downgrade.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index d365424e3..2ca56bc51 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -23,9 +23,9 @@ jobs: - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} - - uses: cjdoris/julia-downgrade-compat-action@v1 + - uses: julia-actions/julia-downgrade-compat@v1 # if: ${{ matrix.version == '1.6' }} with: skip: Pkg,TOML - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-runtest@v1 \ No newline at end of file + - uses: julia-actions/julia-runtest@v1 From c8104cf3020f58fe5a1cfb6c7a23e740c2edcac4 Mon Sep 17 00:00:00 2001 From: Arno Strouwen Date: Thu, 22 Feb 2024 04:43:51 +0100 Subject: [PATCH 51/55] reapplying formatting --- .JuliaFormatter.toml | 1 + benchmarks/applelu.jl | 2 +- benchmarks/lu.jl | 2 +- benchmarks/sparselu.jl | 5 +- docs/pages.jl | 6 +- docs/src/advanced/developing.md | 2 +- ext/LinearSolveBandedMatricesExt.jl | 11 +- ext/LinearSolveBlockDiagonalsExt.jl | 2 +- ext/LinearSolveCUDAExt.jl | 6 +- ext/LinearSolveEnzymeExt.jl | 35 +- ext/LinearSolveFastAlmostBandedMatricesExt.jl | 3 +- ext/LinearSolveHYPREExt.jl | 44 +- ext/LinearSolveIterativeSolversExt.jl | 17 +- ext/LinearSolveKrylovKitExt.jl | 4 +- ext/LinearSolveMetalExt.jl | 6 +- ext/LinearSolvePardisoExt.jl | 20 +- ext/LinearSolveRecursiveArrayToolsExt.jl | 2 +- src/LinearSolve.jl | 24 +- src/appleaccelerate.jl | 69 +- src/common.jl | 30 +- src/default.jl | 74 ++- src/deprecated.jl | 3 +- src/extension_algs.jl | 97 ++- src/factorization.jl | 601 +++++++++--------- src/factorization_sparse.jl | 29 +- src/iterative_wrappers.jl | 27 +- src/mkl.jl | 71 ++- src/simplegmres.jl | 34 +- src/simplelu.jl | 4 +- src/solve_function.jl | 2 +- test/basictests.jl | 11 +- test/default_algs.jl | 20 +- test/enzyme.jl | 57 +- test/gpu/cuda.jl | 10 +- test/hypretests.jl | 4 +- test/pardiso/pardiso.jl | 8 +- test/resolve.jl | 20 +- test/sparse_vector.jl | 4 +- 38 files changed, 718 insertions(+), 649 deletions(-) diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 0f93ea574..3b82401e8 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1,3 +1,4 @@ style = "sciml" format_markdown = true +format_docstrings = true annotate_untyped_fields_with_any = false \ No newline at end of file diff --git a/benchmarks/applelu.jl b/benchmarks/applelu.jl index c346ff0c9..c0250a687 100644 --- a/benchmarks/applelu.jl +++ b/benchmarks/applelu.jl @@ -23,7 +23,7 @@ algs = [ GenericLUFactorization(), RFLUFactorization(), AppleAccelerateLUFactorization(), - MetalLUFactorization(), + MetalLUFactorization() ] res = [Float32[] for i in 1:length(algs)] diff --git a/benchmarks/lu.jl b/benchmarks/lu.jl index 7ab624b45..b8d76a517 100644 --- a/benchmarks/lu.jl +++ b/benchmarks/lu.jl @@ -24,7 +24,7 @@ algs = [ RFLUFactorization(), MKLLUFactorization(), FastLUFactorization(), - SimpleLUFactorization(), + SimpleLUFactorization() ] res = [Float64[] for i in 1:length(algs)] diff --git a/benchmarks/sparselu.jl b/benchmarks/sparselu.jl index db917b1cc..aa714f771 100644 --- a/benchmarks/sparselu.jl +++ b/benchmarks/sparselu.jl @@ -36,7 +36,7 @@ algs = [ UMFPACKFactorization(), KLUFactorization(), MKLPardisoFactorize(), - SparspakFactorization(), + SparspakFactorization() ] cols = [:red, :blue, :green, :magenta, :turqoise] # one color per alg lst = [:dash, :solid, :dashdot] # one line style per dim @@ -65,7 +65,8 @@ function run_and_plot(; dims = [1, 2, 3], kmax = 12) u0 = rand(rng, n) for j in 1:length(algs) - bt = @belapsed solve(prob, $(algs[j])).u setup=(prob = LinearProblem(copy($A), + bt = @belapsed solve(prob, $(algs[j])).u setup=(prob = LinearProblem( + copy($A), copy($b); u0 = copy($u0), alias_A = true, diff --git a/docs/pages.jl b/docs/pages.jl index d8f2e3cda..efb625df1 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,7 +2,7 @@ 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", @@ -10,6 +10,6 @@ pages = ["index.md", "basics/FAQ.md"], "Solvers" => Any["solvers/solvers.md"], "Advanced" => Any["advanced/developing.md" - "advanced/custom.md"], - "Release Notes" => "release_notes.md", + "advanced/custom.md"], + "Release Notes" => "release_notes.md" ] diff --git a/docs/src/advanced/developing.md b/docs/src/advanced/developing.md index f3861173b..c15c8b70e 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/ext/LinearSolveBandedMatricesExt.jl b/ext/LinearSolveBandedMatricesExt.jl index c6d851894..3941b0864 100644 --- a/ext/LinearSolveBandedMatricesExt.jl +++ b/ext/LinearSolveBandedMatricesExt.jl @@ -2,7 +2,8 @@ module LinearSolveBandedMatricesExt using BandedMatrices, LinearAlgebra, LinearSolve import LinearSolve: defaultalg, - do_factorization, init_cacheval, DefaultLinearSolver, DefaultAlgorithmChoice + do_factorization, init_cacheval, DefaultLinearSolver, + DefaultAlgorithmChoice # Defaults for BandedMatrices function defaultalg(A::BandedMatrix, b, oa::OperatorAssumptions{Bool}) @@ -35,14 +36,14 @@ for alg in (:SVDFactorization, :MKLLUFactorization, :DiagonalFactorization, :AppleAccelerateLUFactorization, :CholeskyFactorization) @eval begin function init_cacheval(::$(alg), ::BandedMatrix, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) return nothing end end end function init_cacheval(::LUFactorization, A::BandedMatrix, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) return lu(similar(A, 0, 0)) end @@ -54,8 +55,8 @@ for alg in (:SVDFactorization, :MKLLUFactorization, :DiagonalFactorization, :AppleAccelerateLUFactorization, :QRFactorization, :LUFactorization) @eval begin function init_cacheval(::$(alg), ::Symmetric{<:Number, <:BandedMatrix}, b, u, Pl, - Pr, maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + Pr, maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) return nothing end end diff --git a/ext/LinearSolveBlockDiagonalsExt.jl b/ext/LinearSolveBlockDiagonalsExt.jl index 4b1827200..1109bc4db 100644 --- a/ext/LinearSolveBlockDiagonalsExt.jl +++ b/ext/LinearSolveBlockDiagonalsExt.jl @@ -3,7 +3,7 @@ module LinearSolveBlockDiagonalsExt using LinearSolve, BlockDiagonals function LinearSolve.init_cacheval(alg::SimpleGMRES{false}, A::BlockDiagonal, b, args...; - kwargs...) + kwargs...) @assert ndims(A)==2 "ndims(A) == $(ndims(A)). `A` must have ndims == 2." # We need to perform this check even when `zeroinit == true`, since the type of the # cache is dependent on whether we are able to use the specialized dispatch. diff --git a/ext/LinearSolveCUDAExt.jl b/ext/LinearSolveCUDAExt.jl index c036d53e0..1cfb0888e 100644 --- a/ext/LinearSolveCUDAExt.jl +++ b/ext/LinearSolveCUDAExt.jl @@ -6,7 +6,7 @@ using LinearSolve.LinearAlgebra, LinearSolve.SciMLBase, LinearSolve.ArrayInterfa using SciMLBase: AbstractSciMLOperator function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::CudaOffloadFactorization; - kwargs...) + kwargs...) if cache.isfresh fact = qr(CUDA.CuArray(cache.A)) cache.cacheval = fact @@ -18,8 +18,8 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::CudaOffloadFactor end function LinearSolve.init_cacheval(alg::CudaOffloadFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) qr(CUDA.CuArray(A)) end diff --git a/ext/LinearSolveEnzymeExt.jl b/ext/LinearSolveEnzymeExt.jl index 304e302d2..c8d89e874 100644 --- a/ext/LinearSolveEnzymeExt.jl +++ b/ext/LinearSolveEnzymeExt.jl @@ -4,12 +4,13 @@ using LinearSolve using LinearSolve.LinearAlgebra isdefined(Base, :get_extension) ? (import Enzyme) : (import ..Enzyme) - using Enzyme using EnzymeCore -function EnzymeCore.EnzymeRules.forward(func::Const{typeof(LinearSolve.init)}, ::Type{RT}, prob::EnzymeCore.Annotation{LP}, alg::Const; kwargs...) where {RT, LP <: LinearSolve.LinearProblem} +function EnzymeCore.EnzymeRules.forward( + func::Const{typeof(LinearSolve.init)}, ::Type{RT}, prob::EnzymeCore.Annotation{LP}, + alg::Const; kwargs...) where {RT, LP <: LinearSolve.LinearProblem} @assert !(prob isa Const) res = func.val(prob.val, alg.val; kwargs...) if RT <: Const @@ -26,11 +27,13 @@ function EnzymeCore.EnzymeRules.forward(func::Const{typeof(LinearSolve.init)}, : error("Unsupported return type $RT") end -function EnzymeCore.EnzymeRules.forward(func::Const{typeof(LinearSolve.solve!)}, ::Type{RT}, linsolve::EnzymeCore.Annotation{LP}; kwargs...) where {RT, LP <: LinearSolve.LinearCache} +function EnzymeCore.EnzymeRules.forward(func::Const{typeof(LinearSolve.solve!)}, + ::Type{RT}, linsolve::EnzymeCore.Annotation{LP}; + kwargs...) where {RT, LP <: LinearSolve.LinearCache} @assert !(linsolve isa Const) res = func.val(linsolve.val; kwargs...) - + if RT <: Const return res end @@ -56,7 +59,10 @@ function EnzymeCore.EnzymeRules.forward(func::Const{typeof(LinearSolve.solve!)}, return Duplicated(res, dres) end -function EnzymeCore.EnzymeRules.augmented_primal(config, func::Const{typeof(LinearSolve.init)}, ::Type{RT}, prob::EnzymeCore.Annotation{LP}, alg::Const; kwargs...) where {RT, LP <: LinearSolve.LinearProblem} +function EnzymeCore.EnzymeRules.augmented_primal( + config, func::Const{typeof(LinearSolve.init)}, + ::Type{RT}, prob::EnzymeCore.Annotation{LP}, alg::Const; + kwargs...) where {RT, LP <: LinearSolve.LinearProblem} res = func.val(prob.val, alg.val; kwargs...) dres = if EnzymeRules.width(config) == 1 func.val(prob.dval, alg.val; kwargs...) @@ -77,7 +83,6 @@ function EnzymeCore.EnzymeRules.augmented_primal(config, func::Const{typeof(Line (dval.b for dval in dres) end - prob_d_A = if EnzymeRules.width(config) == 1 prob.dval.A else @@ -92,7 +97,10 @@ function EnzymeCore.EnzymeRules.augmented_primal(config, func::Const{typeof(Line return EnzymeCore.EnzymeRules.AugmentedReturn(res, dres, (d_A, d_b, prob_d_A, prob_d_b)) end -function EnzymeCore.EnzymeRules.reverse(config, func::Const{typeof(LinearSolve.init)}, ::Type{RT}, cache, prob::EnzymeCore.Annotation{LP}, alg::Const; kwargs...) where {RT, LP <: LinearSolve.LinearProblem} +function EnzymeCore.EnzymeRules.reverse( + config, func::Const{typeof(LinearSolve.init)}, ::Type{RT}, + cache, prob::EnzymeCore.Annotation{LP}, alg::Const; + kwargs...) where {RT, LP <: LinearSolve.LinearProblem} d_A, d_b, prob_d_A, prob_d_b = cache if EnzymeRules.width(config) == 1 @@ -105,7 +113,7 @@ function EnzymeCore.EnzymeRules.reverse(config, func::Const{typeof(LinearSolve.i d_b .= 0 end else - for (_prob_d_A,_d_A,_prob_d_b, _d_b) in zip(prob_d_A, d_A, prob_d_b, d_b) + for (_prob_d_A, _d_A, _prob_d_b, _d_b) in zip(prob_d_A, d_A, prob_d_b, d_b) if _d_A !== _prob_d_A _prob_d_A .+= _d_A _d_A .= 0 @@ -123,7 +131,10 @@ end # y=inv(A) B # dA −= z y^T # dB += z, where z = inv(A^T) dy -function EnzymeCore.EnzymeRules.augmented_primal(config, func::Const{typeof(LinearSolve.solve!)}, ::Type{RT}, linsolve::EnzymeCore.Annotation{LP}; kwargs...) where {RT, LP <: LinearSolve.LinearCache} +function EnzymeCore.EnzymeRules.augmented_primal( + config, func::Const{typeof(LinearSolve.solve!)}, + ::Type{RT}, linsolve::EnzymeCore.Annotation{LP}; + kwargs...) where {RT, LP <: LinearSolve.LinearCache} res = func.val(linsolve.val; kwargs...) dres = if EnzymeRules.width(config) == 1 @@ -176,7 +187,9 @@ function EnzymeCore.EnzymeRules.augmented_primal(config, func::Const{typeof(Line return EnzymeCore.EnzymeRules.AugmentedReturn(res, dres, cache) end -function EnzymeCore.EnzymeRules.reverse(config, func::Const{typeof(LinearSolve.solve!)}, ::Type{RT}, cache, linsolve::EnzymeCore.Annotation{LP}; kwargs...) where {RT, LP <: LinearSolve.LinearCache} +function EnzymeCore.EnzymeRules.reverse(config, func::Const{typeof(LinearSolve.solve!)}, + ::Type{RT}, cache, linsolve::EnzymeCore.Annotation{LP}; + kwargs...) where {RT, LP <: LinearSolve.LinearCache} y, dys, _linsolve, dAs, dbs = cache @assert !(linsolve isa Const) @@ -202,7 +215,7 @@ function EnzymeCore.EnzymeRules.reverse(config, func::Const{typeof(LinearSolve.s LinearSolve.defaultalg_adjoint_eval(_linsolve, dy) else error("Algorithm $(_linsolve.alg) is currently not supported by Enzyme rules on LinearSolve.jl. Please open an issue on LinearSolve.jl detailing which algorithm is missing the adjoint handling") - end + end dA .-= z * transpose(y) db .+= z diff --git a/ext/LinearSolveFastAlmostBandedMatricesExt.jl b/ext/LinearSolveFastAlmostBandedMatricesExt.jl index 080ec30f5..1ceff10c5 100644 --- a/ext/LinearSolveFastAlmostBandedMatricesExt.jl +++ b/ext/LinearSolveFastAlmostBandedMatricesExt.jl @@ -2,7 +2,8 @@ module LinearSolveFastAlmostBandedMatricesExt using FastAlmostBandedMatrices, LinearAlgebra, LinearSolve import LinearSolve: defaultalg, - do_factorization, init_cacheval, DefaultLinearSolver, DefaultAlgorithmChoice + do_factorization, init_cacheval, DefaultLinearSolver, + DefaultAlgorithmChoice function defaultalg(A::AlmostBandedMatrix, b, oa::OperatorAssumptions{Bool}) if oa.issq diff --git a/ext/LinearSolveHYPREExt.jl b/ext/LinearSolveHYPREExt.jl index 00f969fe3..cef49b7a8 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) @@ -89,7 +89,7 @@ function SciMLBase.init(prob::LinearProblem, alg::HYPREAlgorithm, cache = LinearCache{ typeof(A), typeof(b), typeof(u0), typeof(p), typeof(alg), Tc, typeof(Pl), typeof(Pr), typeof(reltol), - typeof(__issquare(assumptions)), + typeof(__issquare(assumptions)) }(A, b, u0, p, alg, cacheval, isfresh, Pl, Pr, abstol, reltol, maxiters, verbose, assumptions) @@ -219,8 +219,8 @@ 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 3d73197eb..507e75ad2 100644 --- a/ext/LinearSolveIterativeSolversExt.jl +++ b/ext/LinearSolveIterativeSolversExt.jl @@ -11,8 +11,8 @@ 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) end @@ -49,9 +49,9 @@ 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 s = :idrs_s in keys(alg.kwargs) ? alg.kwargs.idrs_s : 4 # shadow space @@ -69,10 +69,10 @@ function LinearSolve.init_cacheval(alg::IterativeSolversJL, A, b, u, Pl, Pr, max elseif alg.generate_iterator === IterativeSolvers.idrs_iterable! !!LinearSolve._isidentity_struct(Pr) && @warn "$(alg.generate_iterator) doesn't support right preconditioning" - history = IterativeSolvers.ConvergenceHistory(partial=true) + history = IterativeSolvers.ConvergenceHistory(partial = true) history[:abstol] = abstol history[:reltol] = reltol - IterativeSolvers.idrs_iterable!(history, u, A, b, s, Pl, abstol, reltol, maxiters; + IterativeSolvers.idrs_iterable!(history, u, A, b, s, Pl, abstol, reltol, maxiters; alg.kwargs...) elseif alg.generate_iterator === IterativeSolvers.bicgstabl_iterator! !!LinearSolve._isidentity_struct(Pr) && @@ -110,7 +110,8 @@ function SciMLBase.solve!(cache::LinearCache, alg::IterativeSolversJL; kwargs... end cache.verbose && println() - resid = cache.cacheval isa IterativeSolvers.IDRSIterable ? cache.cacheval.R : cache.cacheval.residual + resid = cache.cacheval isa IterativeSolvers.IDRSIterable ? cache.cacheval.R : + cache.cacheval.residual if resid isa IterativeSolvers.Residual resid = resid.current end diff --git a/ext/LinearSolveKrylovKitExt.jl b/ext/LinearSolveKrylovKitExt.jl index 71aa5dddf..a26e26688 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 diff --git a/ext/LinearSolveMetalExt.jl b/ext/LinearSolveMetalExt.jl index 2e963f514..036ffa9cd 100644 --- a/ext/LinearSolveMetalExt.jl +++ b/ext/LinearSolveMetalExt.jl @@ -9,13 +9,13 @@ default_alias_A(::MetalLUFactorization, ::Any, ::Any) = false default_alias_b(::MetalLUFactorization, ::Any, ::Any) = false function LinearSolve.init_cacheval(alg::MetalLUFactorization, 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 SciMLBase.solve!(cache::LinearCache, alg::MetalLUFactorization; - kwargs...) + kwargs...) A = cache.A A = convert(AbstractMatrix, A) if cache.isfresh diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl index da7feaf61..60312ce66 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) diff --git a/ext/LinearSolveRecursiveArrayToolsExt.jl b/ext/LinearSolveRecursiveArrayToolsExt.jl index 98da9583c..468f53e84 100644 --- a/ext/LinearSolveRecursiveArrayToolsExt.jl +++ b/ext/LinearSolveRecursiveArrayToolsExt.jl @@ -5,7 +5,7 @@ import LinearSolve: init_cacheval # Krylov.jl tries to init with `ArrayPartition(undef, ...)`. Avoid hitting that! function init_cacheval(alg::LinearSolve.KrylovJL, A, b::ArrayPartition, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, ::LinearSolve.OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, ::LinearSolve.OperatorAssumptions) return nothing end diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index 7ffaadf9c..efec6be98 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -29,8 +29,8 @@ PrecompileTools.@recompile_invalidations begin using LinearAlgebra: BlasInt, LU using LinearAlgebra.LAPACK: require_one_based_indexing, - chkfinite, chkstride1, - @blasfunc, chkargsok + chkfinite, chkstride1, + @blasfunc, chkargsok import GPUArraysCore import Preferences @@ -128,7 +128,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) @@ -214,19 +214,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, IterativeSolversJL_IDRS, - KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_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 diff --git a/src/appleaccelerate.jl b/src/appleaccelerate.jl index 44bb07dcb..917ad9c4a 100644 --- a/src/appleaccelerate.jl +++ b/src/appleaccelerate.jl @@ -31,9 +31,9 @@ else end function aa_getrf!(A::AbstractMatrix{<:ComplexF64}; - ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), - info = Ref{Cint}(), - check = false) + ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), + info = Ref{Cint}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -51,9 +51,9 @@ function aa_getrf!(A::AbstractMatrix{<:ComplexF64}; end function aa_getrf!(A::AbstractMatrix{<:ComplexF32}; - ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), - info = Ref{Cint}(), - check = false) + ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), + info = Ref{Cint}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -71,9 +71,9 @@ function aa_getrf!(A::AbstractMatrix{<:ComplexF32}; end function aa_getrf!(A::AbstractMatrix{<:Float64}; - ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), - info = Ref{Cint}(), - check = false) + ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), + info = Ref{Cint}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -91,9 +91,9 @@ function aa_getrf!(A::AbstractMatrix{<:Float64}; end function aa_getrf!(A::AbstractMatrix{<:Float32}; - ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), - info = Ref{Cint}(), - check = false) + ipiv = similar(A, Cint, min(size(A, 1), size(A, 2))), + info = Ref{Cint}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -112,10 +112,10 @@ function aa_getrf!(A::AbstractMatrix{<:Float32}; end function aa_getrs!(trans::AbstractChar, - A::AbstractMatrix{<:ComplexF64}, - ipiv::AbstractVector{Cint}, - B::AbstractVecOrMat{<:ComplexF64}; - info = Ref{Cint}()) + A::AbstractMatrix{<:ComplexF64}, + ipiv::AbstractVector{Cint}, + B::AbstractVecOrMat{<:ComplexF64}; + info = Ref{Cint}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -136,10 +136,10 @@ function aa_getrs!(trans::AbstractChar, end function aa_getrs!(trans::AbstractChar, - A::AbstractMatrix{<:ComplexF32}, - ipiv::AbstractVector{Cint}, - B::AbstractVecOrMat{<:ComplexF32}; - info = Ref{Cint}()) + A::AbstractMatrix{<:ComplexF32}, + ipiv::AbstractVector{Cint}, + B::AbstractVecOrMat{<:ComplexF32}; + info = Ref{Cint}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -161,10 +161,10 @@ function aa_getrs!(trans::AbstractChar, end function aa_getrs!(trans::AbstractChar, - A::AbstractMatrix{<:Float64}, - ipiv::AbstractVector{Cint}, - B::AbstractVecOrMat{<:Float64}; - info = Ref{Cint}()) + A::AbstractMatrix{<:Float64}, + ipiv::AbstractVector{Cint}, + B::AbstractVecOrMat{<:Float64}; + info = Ref{Cint}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -186,10 +186,10 @@ function aa_getrs!(trans::AbstractChar, end function aa_getrs!(trans::AbstractChar, - A::AbstractMatrix{<:Float32}, - ipiv::AbstractVector{Cint}, - B::AbstractVecOrMat{<:Float32}; - info = Ref{Cint}()) + A::AbstractMatrix{<:Float32}, + ipiv::AbstractVector{Cint}, + B::AbstractVecOrMat{<:Float32}; + info = Ref{Cint}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -220,21 +220,22 @@ const PREALLOCATED_APPLE_LU = begin end function LinearSolve.init_cacheval(alg::AppleAccelerateLUFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) PREALLOCATED_APPLE_LU end -function LinearSolve.init_cacheval(alg::AppleAccelerateLUFactorization, A::AbstractMatrix{<:Union{Float32,ComplexF32,ComplexF64}}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) +function LinearSolve.init_cacheval(alg::AppleAccelerateLUFactorization, + A::AbstractMatrix{<:Union{Float32, ComplexF32, ComplexF64}}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) A = rand(eltype(A), 0, 0) luinst = ArrayInterface.lu_instance(A) LU(luinst.factors, similar(A, Cint, 0), luinst.info), Ref{Cint}() end function SciMLBase.solve!(cache::LinearCache, alg::AppleAccelerateLUFactorization; - kwargs...) + kwargs...) A = cache.A A = convert(AbstractMatrix, A) if cache.isfresh diff --git a/src/common.jl b/src/common.jl index 1f4d50e2e..a14e976a8 100644 --- a/src/common.jl +++ b/src/common.jl @@ -7,7 +7,7 @@ can be dependent on the condition number of the matrix. The condition number can ```julia using LinearAlgebra -cond(rand(100,100)) +cond(rand(100, 100)) ``` However, in practice this computation is very expensive and thus not possible for most practical cases. @@ -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 @@ -128,17 +128,17 @@ end __init_u0_from_Ab(::SMatrix{S1, S2}, b) where {S1, S2} = zeros(SVector{S2, eltype(b)}) 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(real(eltype(prob.b))), - reltol = default_tol(real(eltype(prob.b))), - 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(real(eltype(prob.b))), + reltol = default_tol(real(eltype(prob.b))), + 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 || A isa SMatrix @@ -181,8 +181,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 9a76403bb..4621edcff 100644 --- a/src/default.jl +++ b/src/default.jl @@ -30,7 +30,7 @@ end newex = quote setfield!(cache, $(Meta.quot(alg)), v) end - alg_enum = getproperty(LinearSolve.DefaultAlgorithmChoice, alg) + alg_enum = getproperty(LinearSolve.DefaultAlgorithmChoice, alg) ex = if ex == :() Expr(:elseif, :(alg.alg == $(alg_enum)), newex, :(error("Algorithm Choice not Allowed"))) @@ -46,7 +46,7 @@ end defaultalg(A, b) = defaultalg(A, b, OperatorAssumptions(true)) function defaultalg(A::Union{DiffEqArrayOperator, MatrixOperator}, b, - assump::OperatorAssumptions{Bool}) + assump::OperatorAssumptions{Bool}) defaultalg(A.A, b, assump) end @@ -92,12 +92,13 @@ function defaultalg(A::Symmetric{<:Number, <:Array}, b, ::OperatorAssumptions{Bo DefaultLinearSolver(DefaultAlgorithmChoice.BunchKaufmanFactorization) end -function defaultalg(A::Symmetric{<:Number, <:SparseMatrixCSC}, b, ::OperatorAssumptions{Bool}) +function defaultalg( + A::Symmetric{<:Number, <:SparseMatrixCSC}, b, ::OperatorAssumptions{Bool}) DefaultLinearSolver(DefaultAlgorithmChoice.CHOLMODFactorization) end function defaultalg(A::AbstractSparseMatrixCSC{Tv, Ti}, b, - assump::OperatorAssumptions{Bool}) where {Tv, Ti} + assump::OperatorAssumptions{Bool}) where {Tv, Ti} if assump.issq DefaultLinearSolver(DefaultAlgorithmChoice.SparspakFactorization) else @@ -107,7 +108,7 @@ end @static if INCLUDE_SPARSE function defaultalg(A::AbstractSparseMatrixCSC{<:Union{Float64, ComplexF64}, Ti}, b, - assump::OperatorAssumptions{Bool}) where {Ti} + assump::OperatorAssumptions{Bool}) where {Ti} if assump.issq if length(b) <= 10_000 && length(nonzeros(A)) / length(A) < 2e-4 DefaultLinearSolver(DefaultAlgorithmChoice.KLUFactorization) @@ -129,7 +130,8 @@ function defaultalg(A::GPUArraysCore.AnyGPUArray, b, assump::OperatorAssumptions end # A === nothing case -function defaultalg(A::Nothing, b::GPUArraysCore.AnyGPUArray, assump::OperatorAssumptions{Bool}) +function defaultalg( + A::Nothing, b::GPUArraysCore.AnyGPUArray, assump::OperatorAssumptions{Bool}) if assump.condition === OperatorCondition.IllConditioned || !assump.issq DefaultLinearSolver(DefaultAlgorithmChoice.QRFactorization) else @@ -139,7 +141,7 @@ end # Ambiguity handling function defaultalg(A::GPUArraysCore.AnyGPUArray, b::GPUArraysCore.AnyGPUArray, - assump::OperatorAssumptions{Bool}) + assump::OperatorAssumptions{Bool}) if assump.condition === OperatorCondition.IllConditioned || !assump.issq DefaultLinearSolver(DefaultAlgorithmChoice.QRFactorization) else @@ -148,7 +150,7 @@ function defaultalg(A::GPUArraysCore.AnyGPUArray, b::GPUArraysCore.AnyGPUArray, end function defaultalg(A::SciMLBase.AbstractSciMLOperator, b, - assump::OperatorAssumptions{Bool}) + assump::OperatorAssumptions{Bool}) if has_ldiv!(A) return DefaultLinearSolver(DefaultAlgorithmChoice.DirectLdiv!) elseif !assump.issq @@ -180,12 +182,12 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) elseif appleaccelerate_isavailable() DefaultAlgorithmChoice.AppleAccelerateLUFactorization elseif (length(b) <= 100 || (isopenblas() && length(b) <= 500) || - (usemkl && length(b) <= 200)) && + (usemkl && length(b) <= 200)) && (A === nothing ? eltype(b) <: Union{Float32, Float64} : eltype(A) <: Union{Float32, Float64}) DefaultAlgorithmChoice.RFLUFactorization - #elseif A === nothing || A isa Matrix - # alg = FastLUFactorization() + #elseif A === nothing || A isa Matrix + # alg = FastLUFactorization() elseif usemkl DefaultAlgorithmChoice.MKLLUFactorization else @@ -196,7 +198,7 @@ function defaultalg(A, b, assump::OperatorAssumptions{Bool}) elseif __conditioning(assump) === OperatorCondition.SuperIllConditioned DefaultAlgorithmChoice.SVDFactorization elseif usemkl && (A === nothing ? eltype(b) <: BLASELTYPES : - eltype(A) <: BLASELTYPES) + eltype(A) <: BLASELTYPES) DefaultAlgorithmChoice.MKLLUFactorization else DefaultAlgorithmChoice.LUFactorization @@ -286,21 +288,22 @@ end ## Catch high level interface function SciMLBase.init(prob::LinearProblem, alg::Nothing, - args...; - assumptions = OperatorAssumptions(issquare(prob.A)), - kwargs...) - SciMLBase.init(prob, defaultalg(prob.A, prob.b, assumptions), args...; assumptions, kwargs...) + args...; + assumptions = OperatorAssumptions(issquare(prob.A)), + kwargs...) + SciMLBase.init( + prob, defaultalg(prob.A, prob.b, assumptions), 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) @@ -310,8 +313,8 @@ 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 || alg === :KrylovJL_CRAIGMR || alg === :KrylovJL_LSMR quote @@ -345,15 +348,15 @@ defaultalg_symbol(::Type{<:QRFactorization{ColumnNorm}}) = :QRFactorizationPivot """ if alg.alg === DefaultAlgorithmChoice.LUFactorization - SciMLBase.solve!(cache, LUFactorization(), args...; kwargs...)) +SciMLBase.solve!(cache, LUFactorization(), args...; kwargs...)) 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 @@ -362,7 +365,7 @@ end retcode = sol.retcode, iters = sol.iters, stats = sol.stats) end - alg_enum = getproperty(LinearSolve.DefaultAlgorithmChoice, alg) + alg_enum = getproperty(LinearSolve.DefaultAlgorithmChoice, alg) ex = if ex == :() Expr(:elseif, :(alg.alg == $(alg_enum)), newex, :(error("Algorithm Choice not Allowed"))) @@ -389,7 +392,7 @@ end DefaultAlgorithmChoice.AppleAccelerateLUFactorization, DefaultAlgorithmChoice.RFLUFactorization)) quote - getproperty(cache.cacheval,$(Meta.quot(alg)))[1]' \ dy + getproperty(cache.cacheval, $(Meta.quot(alg)))[1]' \ dy end elseif alg in Symbol.((DefaultAlgorithmChoice.LUFactorization, DefaultAlgorithmChoice.QRFactorization, @@ -405,9 +408,11 @@ end DefaultAlgorithmChoice.QRFactorizationPivoted, DefaultAlgorithmChoice.GenericLUFactorization)) quote - getproperty(cache.cacheval,$(Meta.quot(alg)))' \ dy + getproperty(cache.cacheval, $(Meta.quot(alg)))' \ dy end - elseif alg in Symbol.((DefaultAlgorithmChoice.KrylovJL_GMRES,DefaultAlgorithmChoice.KrylovJL_LSMR, DefaultAlgorithmChoice.KrylovJL_CRAIGMR)) + elseif alg in Symbol.(( + DefaultAlgorithmChoice.KrylovJL_GMRES, DefaultAlgorithmChoice.KrylovJL_LSMR, + DefaultAlgorithmChoice.KrylovJL_CRAIGMR)) quote invprob = LinearSolve.LinearProblem(transpose(cache.A), dy) solve(invprob, cache.alg; @@ -422,10 +427,15 @@ end end ex = if ex == :() - Expr(:elseif, :(getproperty(DefaultAlgorithmChoice, $(Meta.quot(alg))) === cache.alg.alg), newex, + Expr(:elseif, + :(getproperty(DefaultAlgorithmChoice, $(Meta.quot(alg))) === cache.alg.alg), + newex, :(error("Algorithm Choice not Allowed"))) else - Expr(:elseif, :(getproperty(DefaultAlgorithmChoice, $(Meta.quot(alg))) === cache.alg.alg), newex, ex) + Expr(:elseif, + :(getproperty(DefaultAlgorithmChoice, $(Meta.quot(alg))) === cache.alg.alg), + newex, + ex) end end ex = Expr(:if, ex.args...) diff --git a/src/deprecated.jl b/src/deprecated.jl index a9bd5a846..19623c4c0 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -79,6 +79,7 @@ function set_cacheval(cache::LinearCache, alg_cache) end end -@deprecate SciMLBase.solve(cache::LinearCache, args...; kwargs...) SciMLBase.solve!(cache::LinearCache, +@deprecate SciMLBase.solve(cache::LinearCache, args...; kwargs...) SciMLBase.solve!( + cache::LinearCache, args...; kwargs...) false diff --git a/src/extension_algs.jl b/src/extension_algs.jl index d2e968911..4c2d4fecb 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -15,7 +15,7 @@ alternatively pass an already created solver to `HYPREAlgorithm` (and to the `Pl argument). See HYPRE.jl docs for how to set up solvers with specific options. !!! note - + Using HYPRE solvers requires Julia version 1.9 or higher, and that the package HYPRE.jl is installed. @@ -23,20 +23,20 @@ argument). See HYPRE.jl docs for how to set up solvers with specific options. The single positional argument `solver` has the following choices: -- `HYPRE.BiCGSTAB` -- `HYPRE.BoomerAMG` -- `HYPRE.FlexGMRES` -- `HYPRE.GMRES` -- `HYPRE.Hybrid` -- `HYPRE.ILU` -- `HYPRE.ParaSails` (as preconditioner only) -- `HYPRE.PCG` + - `HYPRE.BiCGSTAB` + - `HYPRE.BoomerAMG` + - `HYPRE.FlexGMRES` + - `HYPRE.GMRES` + - `HYPRE.Hybrid` + - `HYPRE.ILU` + - `HYPRE.ParaSails` (as preconditioner only) + - `HYPRE.PCG` ## Keyword Arguments -* `Pl`: A choice of left preconditioner. + - `Pl`: A choice of left preconditioner. -## Example +## Example For example, to use `HYPRE.PCG` as the solver, with `HYPRE.BoomerAMG` as the preconditioner, the algorithm should be defined as follows: @@ -64,11 +64,11 @@ end """ `CudaOffloadFactorization()` -An offloading technique used to GPU-accelerate CPU-based computations. +An offloading technique used to GPU-accelerate CPU-based computations. Requires a sufficiently large `A` to overcome the data transfer costs. !!! note - + Using this solver requires adding the package CUDA.jl, i.e. `using CUDA` """ struct CudaOffloadFactorization <: LinearSolve.AbstractFactorization @@ -85,15 +85,15 @@ end """ ```julia MKLPardisoFactorize(; nprocs::Union{Int, Nothing} = nothing, - matrix_type = nothing, - iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, - dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) + matrix_type = nothing, + iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, + dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) ``` A sparse factorization method using MKL Pardiso. !!! note - + Using this solver requires adding the package Pardiso.jl, i.e. `using Pardiso` ## Keyword Arguments @@ -108,15 +108,15 @@ MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...) """ ```julia MKLPardisoIterate(; nprocs::Union{Int, Nothing} = nothing, - matrix_type = nothing, - iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, - dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) + matrix_type = nothing, + iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing, + dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing) ``` A mixed factorization+iterative method using MKL Pardiso. !!! note - + Using this solver requires adding the package Pardiso.jl, i.e. `using Pardiso` ## Keyword Arguments @@ -131,10 +131,10 @@ MKLPardisoIterate(; kwargs...) = PardisoJL(; solver_type = 1, kwargs...) """ ```julia 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) ``` A generic method using MKL Pardiso. Specifying `solver_type` is required. @@ -158,10 +158,10 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm 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`") @@ -184,7 +184,7 @@ A generic iterative solver implementation allowing the choice of KrylovKit.jl solvers. !!! note - + Using this solver requires adding the package KrylovKit.jl, i.e. `using KrylovKit` """ struct KrylovKitJL{F, A, I, K} <: LinearSolve.AbstractKrylovSubspaceMethod @@ -202,8 +202,8 @@ KrylovKitJL_CG(args...; Pl = nothing, Pr = nothing, kwargs...) A generic CG implementation for Hermitian and positive definite linear systems !!! note - - Using this solver requires adding the package KrylovKit.jl, i.e. `using KrylovKit` + + Using this solver requires adding the package KrylovKit.jl, i.e. `using KrylovKit` """ function KrylovKitJL_CG end @@ -215,7 +215,7 @@ KrylovKitJL_GMRES(args...; Pl = nothing, Pr = nothing, gmres_restart = 0, kwargs A generic GMRES implementation. !!! note - + Using this solver requires adding the package KrylovKit.jl, i.e. `using KrylovKit` """ function KrylovKitJL_GMRES end @@ -223,18 +223,16 @@ function KrylovKitJL_GMRES end """ ```julia IterativeSolversJL(args...; - generate_iterator = IterativeSolvers.gmres_iterable!, - Pl = nothing, Pr = nothing, - gmres_restart = 0, kwargs...) + generate_iterator = IterativeSolvers.gmres_iterable!, + Pl = nothing, Pr = nothing, + gmres_restart = 0, kwargs...) ``` A generic wrapper over the IterativeSolvers.jl solvers. - !!! note - - Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` + Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` """ struct IterativeSolversJL{F, I, A, K} <: LinearSolve.AbstractKrylovSubspaceMethod generate_iterator::F @@ -251,24 +249,21 @@ IterativeSolversJL_CG(args...; Pl = nothing, Pr = nothing, kwargs...) A wrapper over the IterativeSolvers.jl CG. !!! note - - Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` + Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` """ function IterativeSolversJL_CG end """ ```julia -IterativeSolversJL_GMRES(args...; Pl = nothing, Pr = nothing, gmres_restart=0, kwargs...) +IterativeSolversJL_GMRES(args...; Pl = nothing, Pr = nothing, gmres_restart = 0, kwargs...) ``` A wrapper over the IterativeSolvers.jl GMRES. - !!! note - - Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` + Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` """ function IterativeSolversJL_GMRES end @@ -279,11 +274,9 @@ IterativeSolversJL_IDRS(args...; Pl = nothing, kwargs...) A wrapper over the IterativeSolvers.jl IDR(S). - !!! note - - Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` + Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` """ function IterativeSolversJL_IDRS end @@ -294,11 +287,9 @@ IterativeSolversJL_BICGSTAB(args...; Pl = nothing, Pr = nothing, kwargs...) A wrapper over the IterativeSolvers.jl BICGSTAB. - !!! note - - Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` + Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` """ function IterativeSolversJL_BICGSTAB end @@ -309,11 +300,9 @@ IterativeSolversJL_MINRES(args...; Pl = nothing, Pr = nothing, kwargs...) A wrapper over the IterativeSolvers.jl MINRES. - !!! note - - Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` + Using this solver requires adding the package IterativeSolvers.jl, i.e. `using IterativeSolvers` """ function IterativeSolversJL_MINRES end diff --git a/src/factorization.jl b/src/factorization.jl index 1f244b3c4..0ef7b754e 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -25,7 +25,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 @@ -36,18 +36,18 @@ end Julia's built in `lu`. Equivalent to calling `lu!(A)` -* On dense matrices, this uses the current BLAS implementation of the user's computer, -which by default is OpenBLAS but will use MKL if the user does `using MKL` in their -system. -* On sparse matrices, this will use UMFPACK from SparseArrays. Note that this will not -cache the symbolic factorization. -* On CuMatrix, it will use a CUDA-accelerated LU from CuSolver. -* On BandedMatrix and BlockBandedMatrix, it will use a banded LU. + - On dense matrices, this uses the current BLAS implementation of the user's computer, + which by default is OpenBLAS but will use MKL if the user does `using MKL` in their + system. + - On sparse matrices, this will use UMFPACK from SparseArrays. Note that this will not + cache the symbolic factorization. + - On CuMatrix, it will use a CUDA-accelerated LU from CuSolver. + - On BandedMatrix and BlockBandedMatrix, it will use a banded LU. ## Positional Arguments -* pivot: The choice of pivoting. Defaults to `LinearAlgebra.RowMaximum()`. The other choice is - `LinearAlgebra.NoPivot()`. + - pivot: The choice of pivoting. Defaults to `LinearAlgebra.RowMaximum()`. The other choice is + `LinearAlgebra.NoPivot()`. """ struct LUFactorization{P} <: AbstractFactorization pivot::P @@ -62,8 +62,8 @@ Has low overhead and is good for small matrices. ## Positional Arguments -* pivot: The choice of pivoting. Defaults to `LinearAlgebra.RowMaximum()`. The other choice is - `LinearAlgebra.NoPivot()`. + - pivot: The choice of pivoting. Defaults to `LinearAlgebra.RowMaximum()`. The other choice is + `LinearAlgebra.NoPivot()`. """ struct GenericLUFactorization{P} <: AbstractFactorization pivot::P @@ -94,36 +94,37 @@ function do_factorization(alg::GenericLUFactorization, A, b, u) return fact end -function init_cacheval(alg::Union{LUFactorization, GenericLUFactorization}, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) +function init_cacheval( + alg::Union{LUFactorization, GenericLUFactorization}, A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.lu_instance(convert(AbstractMatrix, A)) end function init_cacheval(alg::Union{LUFactorization, GenericLUFactorization}, - A::Union{<:Adjoint, <:Transpose}, b, u, Pl, Pr, maxiters::Int, abstol, reltol, - verbose::Bool, assumptions::OperatorAssumptions) + A::Union{<:Adjoint, <:Transpose}, b, u, Pl, Pr, maxiters::Int, abstol, reltol, + verbose::Bool, assumptions::OperatorAssumptions) if alg isa LUFactorization - return lu(A; check=false) + return lu(A; check = false) else A isa GPUArraysCore.AnyGPUArray && return nothing - return LinearAlgebra.generic_lufact!(copy(A), alg.pivot; check=false) + return LinearAlgebra.generic_lufact!(copy(A), alg.pivot; check = false) end 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 @@ -134,12 +135,12 @@ end Julia's built in `qr`. Equivalent to calling `qr!(A)`. -* On dense matrices, this uses the current BLAS implementation of the user's computer -which by default is OpenBLAS but will use MKL if the user does `using MKL` in their -system. -* On sparse matrices, this will use SPQR from SparseArrays -* On CuMatrix, it will use a CUDA-accelerated QR from CuSolver. -* On BandedMatrix and BlockBandedMatrix, it will use a banded QR. + - On dense matrices, this uses the current BLAS implementation of the user's computer + which by default is OpenBLAS but will use MKL if the user does `using MKL` in their + system. + - On sparse matrices, this will use SPQR from SparseArrays + - On CuMatrix, it will use a CUDA-accelerated QR from CuSolver. + - On BandedMatrix and BlockBandedMatrix, it will use a banded QR. """ struct QRFactorization{P} <: AbstractFactorization pivot::P @@ -168,20 +169,21 @@ 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), alg.pivot) end const PREALLOCATED_QR_ColumnNorm = ArrayInterface.qr_instance(rand(1, 1), ColumnNorm()) function init_cacheval(alg::QRFactorization{ColumnNorm}, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) return PREALLOCATED_QR_ColumnNorm end -function init_cacheval(alg::QRFactorization, A::Union{<:Adjoint, <:Transpose}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::QRFactorization, A::Union{<:Adjoint, <:Transpose}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) A isa GPUArraysCore.AnyGPUArray && return qr(A) return qr(A, alg.pivot) end @@ -189,13 +191,13 @@ end const PREALLOCATED_QR_NoPivot = ArrayInterface.qr_instance(rand(1, 1)) function init_cacheval(alg::QRFactorization{NoPivot}, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) return PREALLOCATED_QR_NoPivot 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 @@ -208,10 +210,10 @@ Julia's built in `cholesky`. Equivalent to calling `cholesky!(A)`. ## Keyword Arguments -* pivot: defaluts to NoPivot, can also be RowMaximum. -* tol: the tol argument in CHOLMOD. Only used for sparse matrices. -* shift: the shift argument in CHOLMOD. Only used for sparse matrices. -* perm: the perm argument in CHOLMOD. Only used for sparse matrices. + - pivot: defaluts to NoPivot, can also be RowMaximum. + - tol: the tol argument in CHOLMOD. Only used for sparse matrices. + - shift: the shift argument in CHOLMOD. Only used for sparse matrices. + - perm: the perm argument in CHOLMOD. Only used for sparse matrices. """ struct CholeskyFactorization{P, P2} <: AbstractFactorization pivot::P @@ -240,33 +242,33 @@ function do_factorization(alg::CholeskyFactorization, A, b, u) end function init_cacheval(alg::CholeskyFactorization, A::SMatrix{S1, S2}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {S1, S2} + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {S1, S2} cholesky(A) end function init_cacheval(alg::CholeskyFactorization, A::GPUArraysCore.AnyGPUArray, b, u, Pl, - Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) - cholesky(A; check=false) + Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + cholesky(A; check = false) 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 const PREALLOCATED_CHOLESKY = ArrayInterface.cholesky_instance(rand(1, 1), NoPivot()) 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) + A::Union{Diagonal, AbstractSciMLOperator}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) nothing end @@ -292,14 +294,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 @@ -310,9 +312,9 @@ end Julia's built in `svd`. Equivalent to `svd!(A)`. -* On dense matrices, this uses the current BLAS implementation of the user's computer -which by default is OpenBLAS but will use MKL if the user does `using MKL` in their -system. + - On dense matrices, this uses the current BLAS implementation of the user's computer + which by default is OpenBLAS but will use MKL if the user does `using MKL` in their + system. """ struct SVDFactorization{A} <: AbstractFactorization full::Bool @@ -332,22 +334,22 @@ function do_factorization(alg::SVDFactorization, A, b, u) end function init_cacheval(alg::SVDFactorization, A::Union{Matrix, SMatrix}, 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 @@ -361,7 +363,7 @@ Only for Symmetric matrices. ## Keyword Arguments -* rook: whether to perform rook pivoting. Defaults to false. + - rook: whether to perform rook pivoting. Defaults to false. """ Base.@kwdef struct BunchKaufmanFactorization <: AbstractFactorization rook::Bool = false @@ -374,9 +376,9 @@ 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 @@ -384,15 +386,15 @@ const PREALLOCATED_BUNCHKAUFMAN = ArrayInterface.bunchkaufman_instance(Symmetric 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 @@ -400,8 +402,8 @@ end """ `GenericFactorization(;fact_alg=LinearAlgebra.factorize)`: Constructs a linear solver from a generic - factorization algorithm `fact_alg` which complies with the Base.LinearAlgebra - factorization API. Quoting from Base: +factorization algorithm `fact_alg` which complies with the Base.LinearAlgebra +factorization API. Quoting from Base: * If `A` is upper or lower triangular (or diagonal), no factorization of `A` is required. The system is then solved with either forward or backward substitution. @@ -414,8 +416,8 @@ end ## Keyword Arguments -* fact_alg: the factorization algorithm to use. Defaults to `LinearAlgebra.factorize`, but can be - swapped to choices like `lu`, `qr` + - fact_alg: the factorization algorithm to use. Defaults to `LinearAlgebra.factorize`, but can be + swapped to choices like `lu`, `qr` """ struct GenericFactorization{F} <: AbstractFactorization fact_alg::F @@ -429,208 +431,225 @@ function do_factorization(alg::GenericFactorization, A, b, u) return fact end -function init_cacheval(alg::GenericFactorization{typeof(lu)}, A::AbstractMatrix, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(lu)}, A::AbstractMatrix, 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::AbstractMatrix, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(lu!)}, A::AbstractMatrix, 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::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) +function init_cacheval( + alg::GenericFactorization{typeof(lu!)}, A::Tridiagonal, 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::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(lu!)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(A) end -function init_cacheval(alg::GenericFactorization{typeof(lu)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(lu)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(A) end -function init_cacheval(alg::GenericFactorization{typeof(qr)}, A::AbstractMatrix, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(qr)}, A::AbstractMatrix, 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::AbstractMatrix, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(qr!)}, A::AbstractMatrix, 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::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(qr)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(A) end -function init_cacheval(alg::GenericFactorization{typeof(qr!)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(qr!)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(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) +function init_cacheval( + alg::GenericFactorization{typeof(qr!)}, A::Tridiagonal, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.qr_instance(A) end -function init_cacheval(alg::GenericFactorization{typeof(svd)}, A::AbstractMatrix, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(svd)}, A::AbstractMatrix, 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::AbstractMatrix, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(svd!)}, A::AbstractMatrix, 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::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) +function init_cacheval( + alg::GenericFactorization{typeof(svd)}, A::Tridiagonal, 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) + Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) ArrayInterface.svd_instance(A) end -function init_cacheval(alg::GenericFactorization{typeof(svd!)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(svd!)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(A) end -function init_cacheval(alg::GenericFactorization{typeof(svd)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(svd)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(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::GenericFactorization, A, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) end function init_cacheval(alg::GenericFactorization, A::AbstractMatrix, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) do_factorization(alg, A, b, u) 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) +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) 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) +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) if eltype(A) <: Complex return bunchkaufman!(Hermitian(A)) else @@ -642,53 +661,60 @@ end # Try to never use it. # Cholesky needs the posdef matrix, for GenericFactorization assume structure is needed -function init_cacheval(alg::GenericFactorization{typeof(cholesky)}, A::AbstractMatrix, b, u, Pl, Pr, +function init_cacheval( + alg::GenericFactorization{typeof(cholesky)}, A::AbstractMatrix, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + assumptions::OperatorAssumptions) newA = copy(convert(AbstractMatrix, A)) do_factorization(alg, newA, b, u) end -function init_cacheval(alg::GenericFactorization{typeof(cholesky!)}, A::AbstractMatrix, b, u, Pl, Pr, +function init_cacheval( + alg::GenericFactorization{typeof(cholesky!)}, A::AbstractMatrix, b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + assumptions::OperatorAssumptions) newA = copy(convert(AbstractMatrix, A)) do_factorization(alg, newA, b, u) end -function init_cacheval(alg::GenericFactorization{typeof(cholesky!)}, A::Diagonal, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval(alg::GenericFactorization{typeof(cholesky!)}, + A::Diagonal, b, u, Pl, Pr, maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end -function init_cacheval(alg::GenericFactorization{typeof(cholesky!)}, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(cholesky!)}, A::Tridiagonal, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.lu_instance(A) end -function init_cacheval(alg::GenericFactorization{typeof(cholesky!)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(cholesky!)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(A) end -function init_cacheval(alg::GenericFactorization{typeof(cholesky)}, A::Diagonal, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval(alg::GenericFactorization{typeof(cholesky)}, + A::Diagonal, b, u, Pl, Pr, maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) Diagonal(inv.(A.diag)) end -function init_cacheval(alg::GenericFactorization{typeof(cholesky)}, A::Tridiagonal, b, u, Pl, Pr, - maxiters::Int, - abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) +function init_cacheval( + alg::GenericFactorization{typeof(cholesky)}, A::Tridiagonal, b, u, Pl, Pr, + maxiters::Int, + abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) ArrayInterface.lu_instance(A) end -function init_cacheval(alg::GenericFactorization{typeof(cholesky)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) where {T, V} +function init_cacheval( + alg::GenericFactorization{typeof(cholesky)}, A::SymTridiagonal{T, V}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) where {T, V} LinearAlgebra.LDLt{T, SymTridiagonal{T, V}}(A) end - function init_cacheval(alg::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 @@ -719,23 +745,23 @@ const PREALLOCATED_UMFPACK = SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(0, 0 Int[], Float64[])) 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) + maxiters::Int, abstol, + reltol, + verbose::Bool, assumptions::OperatorAssumptions) A = convert(AbstractMatrix, A) zerobased = SparseArrays.getcolptr(A)[1] == 0 return SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A), @@ -792,23 +818,23 @@ const PREALLOCATED_KLU = KLU.KLUFactorization(SparseMatrixCSC(0, 0, [1], Int[], 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) + maxiters::Int, abstol, + reltol, + verbose::Bool, assumptions::OperatorAssumptions) A = convert(AbstractMatrix, A) return KLU.KLUFactorization(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))) @@ -865,8 +891,8 @@ Only supports sparse matrices. ## Keyword Arguments -* shift: the shift argument in CHOLMOD. -* perm: the perm argument in CHOLMOD + - shift: the shift argument in CHOLMOD. + - perm: the perm argument in CHOLMOD """ Base.@kwdef struct CHOLMODFactorization{T} <: AbstractFactorization shift::Float64 = 0.0 @@ -876,17 +902,18 @@ end const PREALLOCATED_CHOLMOD = cholesky(SparseMatrixCSC(0, 0, [1], Int[], 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}} + 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 @@ -928,34 +955,34 @@ 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 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 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) @@ -963,7 +990,7 @@ function SciMLBase.solve!(cache::LinearCache, alg::RFLUFactorization{P, T}; if length(ipiv) != min(size(A)...) ipiv = Vector{LinearAlgebra.BlasInt}(undef, min(size(A)...)) end - fact = RecursiveFactorization.lu!(A, ipiv, Val(P), Val(T), check=false) + fact = RecursiveFactorization.lu!(A, ipiv, Val(P), Val(T), check = false) cache.cacheval = (fact, ipiv) cache.isfresh = false end @@ -982,7 +1009,7 @@ be applied to well-conditioned matrices. ## Positional Arguments -* pivot: Defaults to RowMaximum(), but can be NoPivot() + - pivot: Defaults to RowMaximum(), but can be NoPivot() """ struct NormalCholeskyFactorization{P} <: AbstractFactorization pivot::P @@ -999,38 +1026,39 @@ default_alias_b(::NormalCholeskyFactorization, ::Any, ::Any) = true const PREALLOCATED_NORMALCHOLESKY = ArrayInterface.cholesky_instance(rand(1, 1), NoPivot()) function init_cacheval(alg::NormalCholeskyFactorization, - A::Union{AbstractSparseArray, GPUArraysCore.AnyGPUArray, - Symmetric{<:Number, <:AbstractSparseArray}}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + A::Union{AbstractSparseArray, GPUArraysCore.AnyGPUArray, + 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::SMatrix, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) return cholesky(Symmetric((A)' * 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) A_ = convert(AbstractMatrix, A) - return ArrayInterface.cholesky_instance(Symmetric(Matrix{eltype(A)}(undef,0,0)), alg.pivot) + return ArrayInterface.cholesky_instance( + Symmetric(Matrix{eltype(A)}(undef, 0, 0)), alg.pivot) end const PREALLOCATED_NORMALCHOLESKY_SYMMETRIC = ArrayInterface.cholesky_instance( Symmetric(rand(1, 1)), NoPivot()) function init_cacheval(alg::NormalCholeskyFactorization, A::Matrix{Float64}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) return PREALLOCATED_NORMALCHOLESKY_SYMMETRIC 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 @@ -1069,7 +1097,7 @@ be applied to well-conditioned matrices. ## Positional Arguments -* rook: whether to perform rook pivoting. Defaults to false. + - rook: whether to perform rook pivoting. Defaults to false. """ struct NormalBunchKaufmanFactorization <: AbstractFactorization rook::Bool @@ -1083,13 +1111,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 @@ -1111,12 +1139,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) @@ -1147,8 +1175,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 @@ -1184,30 +1212,30 @@ end FastQRFactorization() = FastQRFactorization(NoPivot(), 36) function init_cacheval(alg::FastQRFactorization{NoPivot}, A::AbstractMatrix, 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))) end function init_cacheval(::FastQRFactorization{ColumnNorm}, A::AbstractMatrix, 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))) end function init_cacheval(alg::FastQRFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) return init_cacheval(alg, convert(AbstractMatrix, A), b, u, Pl, Pr, maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) 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) @@ -1215,10 +1243,12 @@ function SciMLBase.solve!(cache::LinearCache, alg::FastQRFactorization{P}; # 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. if P === NoPivot - @set! ws_and_fact.factors = LinearAlgebra.QRCompactWY(LAPACK.geqrt!(ws_and_fact.workspace, + @set! ws_and_fact.factors = LinearAlgebra.QRCompactWY(LAPACK.geqrt!( + ws_and_fact.workspace, A)...) else - @set! ws_and_fact.factors = LinearAlgebra.QRPivoted(LAPACK.geqp3!(ws_and_fact.workspace, + @set! ws_and_fact.factors = LinearAlgebra.QRPivoted(LAPACK.geqp3!( + ws_and_fact.workspace, A)...) end cache.cacheval = ws_and_fact @@ -1254,25 +1284,26 @@ const PREALLOCATED_SPARSEPAK = sparspaklu(SparseMatrixCSC(0, 0, [1], Int[], Floa 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 A isa SparseArrays.AbstractSparseArray - return sparspaklu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), + return sparspaklu( + SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), factorize = false) else @@ -1282,7 +1313,7 @@ function init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int, end function init_cacheval(::SparspakFactorization, ::StaticArray, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, assumptions::OperatorAssumptions) nothing end @@ -1306,8 +1337,8 @@ end for alg in InteractiveUtils.subtypes(AbstractFactorization) @eval function init_cacheval(alg::$alg, A::MatrixOperator, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + 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) diff --git a/src/factorization_sparse.jl b/src/factorization_sparse.jl index 1bc0ac261..7c2ddb156 100644 --- a/src/factorization_sparse.jl +++ b/src/factorization_sparse.jl @@ -1,22 +1,29 @@ # 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, - SparseArrays.SPQR.QRSparse, - SparseArrays.CHOLMOD.Factor}, b::Vector) + A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY, + SparseArrays.SPQR.QRSparse, + SparseArrays.CHOLMOD.Factor}, b::Vector) x .= A \ b end function _ldiv!(x::AbstractVector, - A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY, - SparseArrays.SPQR.QRSparse, - SparseArrays.CHOLMOD.Factor}, b::AbstractVector) + A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY, + SparseArrays.SPQR.QRSparse, + SparseArrays.CHOLMOD.Factor}, b::AbstractVector) x .= A \ b end # Ambiguity removal -_ldiv!(::SVector, - A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, - b::AbstractVector) = (A \ b) -_ldiv!(::SVector, A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, - b::SVector) = (A \ b) +function _ldiv!(::SVector, + A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, + LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, + b::AbstractVector) + (A \ b) +end +function _ldiv!(::SVector, + A::Union{SparseArrays.CHOLMOD.Factor, LinearAlgebra.QR, + LinearAlgebra.QRCompactWY, SparseArrays.SPQR.QRSparse}, + b::SVector) + (A \ b) +end diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl index dbdef5b9b..3c03ea7d6 100644 --- a/src/iterative_wrappers.jl +++ b/src/iterative_wrappers.jl @@ -3,9 +3,9 @@ """ ```julia KrylovJL(args...; KrylovAlg = Krylov.gmres!, - Pl = nothing, Pr = nothing, - gmres_restart = 0, window = 0, - kwargs...) + Pl = nothing, Pr = nothing, + gmres_restart = 0, window = 0, + kwargs...) ``` A generic wrapper over the Krylov.jl krylov-subspace iterative solvers. @@ -19,8 +19,8 @@ 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) end @@ -30,7 +30,7 @@ default_alias_b(::KrylovJL, ::Any, ::Any) = true """ ```julia -KrylovJL_CG(args...; kwargs...) +KrylovJL_CG(args...; kwargs...) ``` A generic CG implementation for Hermitian and positive definite linear systems @@ -41,7 +41,7 @@ end """ ```julia -KrylovJL_MINRES(args...; kwargs...) +KrylovJL_MINRES(args...; kwargs...) ``` A generic MINRES implementation for Hermitian linear systems @@ -52,7 +52,7 @@ end """ ```julia -KrylovJL_GMRES(args...; gmres_restart = 0, window = 0, kwargs...) +KrylovJL_GMRES(args...; gmres_restart = 0, window = 0, kwargs...) ``` A generic GMRES implementation for square non-Hermitian linear systems @@ -63,7 +63,7 @@ end """ ```julia -KrylovJL_BICGSTAB(args...; kwargs...) +KrylovJL_BICGSTAB(args...; kwargs...) ``` A generic BICGSTAB implementation for square non-Hermitian linear systems @@ -74,7 +74,7 @@ end """ ```julia -KrylovJL_LSMR(args...; kwargs...) +KrylovJL_LSMR(args...; kwargs...) ``` A generic LSMR implementation for least-squares problems @@ -85,7 +85,7 @@ end """ ```julia -KrylovJL_CRAIGMR(args...; kwargs...) +KrylovJL_CRAIGMR(args...; kwargs...) ``` A generic CRAIGMR implementation for least-norm problems @@ -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 @@ -282,7 +282,8 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...) end stats = @get_cacheval(cache, :KrylovJL_GMRES).stats - resid = !isempty(stats.residuals) ? last(stats.residuals) : zero(eltype(stats.residuals)) + resid = !isempty(stats.residuals) ? last(stats.residuals) : + zero(eltype(stats.residuals)) retcode = if !stats.solved if stats.status == "maximum number of iterations exceeded" diff --git a/src/mkl.jl b/src/mkl.jl index db3106c08..a00882e1d 100644 --- a/src/mkl.jl +++ b/src/mkl.jl @@ -9,9 +9,9 @@ to avoid allocations and does not require libblastrampoline. struct MKLLUFactorization <: AbstractFactorization end function getrf!(A::AbstractMatrix{<:ComplexF64}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), - check = false) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -29,9 +29,9 @@ function getrf!(A::AbstractMatrix{<:ComplexF64}; end function getrf!(A::AbstractMatrix{<:ComplexF32}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), - check = false) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -49,9 +49,9 @@ function getrf!(A::AbstractMatrix{<:ComplexF32}; end function getrf!(A::AbstractMatrix{<:Float64}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), - check = false) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -69,9 +69,9 @@ function getrf!(A::AbstractMatrix{<:Float64}; end function getrf!(A::AbstractMatrix{<:Float32}; - ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), - info = Ref{BlasInt}(), - check = false) + ipiv = similar(A, BlasInt, min(size(A, 1), size(A, 2))), + info = Ref{BlasInt}(), + check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) @@ -89,10 +89,10 @@ function getrf!(A::AbstractMatrix{<:Float32}; end function getrs!(trans::AbstractChar, - A::AbstractMatrix{<:ComplexF64}, - ipiv::AbstractVector{BlasInt}, - B::AbstractVecOrMat{<:ComplexF64}; - info = Ref{BlasInt}()) + A::AbstractMatrix{<:ComplexF64}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:ComplexF64}; + info = Ref{BlasInt}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -114,10 +114,10 @@ function getrs!(trans::AbstractChar, end function getrs!(trans::AbstractChar, - A::AbstractMatrix{<:ComplexF32}, - ipiv::AbstractVector{BlasInt}, - B::AbstractVecOrMat{<:ComplexF32}; - info = Ref{BlasInt}()) + A::AbstractMatrix{<:ComplexF32}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:ComplexF32}; + info = Ref{BlasInt}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -139,10 +139,10 @@ function getrs!(trans::AbstractChar, end function getrs!(trans::AbstractChar, - A::AbstractMatrix{<:Float64}, - ipiv::AbstractVector{BlasInt}, - B::AbstractVecOrMat{<:Float64}; - info = Ref{BlasInt}()) + A::AbstractMatrix{<:Float64}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:Float64}; + info = Ref{BlasInt}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -164,10 +164,10 @@ function getrs!(trans::AbstractChar, end function getrs!(trans::AbstractChar, - A::AbstractMatrix{<:Float32}, - ipiv::AbstractVector{BlasInt}, - B::AbstractVecOrMat{<:Float32}; - info = Ref{BlasInt}()) + A::AbstractMatrix{<:Float32}, + ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{<:Float32}; + info = Ref{BlasInt}()) require_one_based_indexing(A, ipiv, B) LinearAlgebra.LAPACK.chktrans(trans) chkstride1(A, B, ipiv) @@ -197,20 +197,21 @@ const PREALLOCATED_MKL_LU = begin end function LinearSolve.init_cacheval(alg::MKLLUFactorization, A, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) PREALLOCATED_MKL_LU end -function LinearSolve.init_cacheval(alg::MKLLUFactorization, A::AbstractMatrix{<:Union{Float32,ComplexF32,ComplexF64}}, b, u, Pl, Pr, - maxiters::Int, abstol, reltol, verbose::Bool, - assumptions::OperatorAssumptions) +function LinearSolve.init_cacheval(alg::MKLLUFactorization, + A::AbstractMatrix{<:Union{Float32, ComplexF32, ComplexF64}}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, verbose::Bool, + assumptions::OperatorAssumptions) A = rand(eltype(A), 0, 0) ArrayInterface.lu_instance(A), Ref{BlasInt}() end function SciMLBase.solve!(cache::LinearCache, alg::MKLLUFactorization; - kwargs...) + kwargs...) A = cache.A A = convert(AbstractMatrix, A) if cache.isfresh @@ -239,4 +240,4 @@ function SciMLBase.solve!(cache::LinearCache, alg::MKLLUFactorization; SciMLBase.build_linear_solution(alg, cache.u, nothing, cache) =# -end \ No newline at end of file +end diff --git a/src/simplegmres.jl b/src/simplegmres.jl index 76682e547..c2596efc2 100644 --- a/src/simplegmres.jl +++ b/src/simplegmres.jl @@ -9,14 +9,16 @@ specialized dispatches. ## Arguments -* `restart::Bool`: If `true`, then the solver will restart after `memory` iterations. -* `memory::Int = 20`: The number of iterations before restarting. If restart is false, this - value is used to allocate memory and later expanded if more memory is required. -* `blocksize::Int = 0`: If blocksize is `> 0`, the solver assumes that the matrix has a - uniformly sized block diagonal structure with square blocks of size `blocksize`. Misusing - this option will lead to incorrect results. - * If this is set `≤ 0` and during runtime we get a Block Diagonal Matrix, then we will - check if the specialized dispatch can be used. + - `restart::Bool`: If `true`, then the solver will restart after `memory` iterations. + + - `memory::Int = 20`: The number of iterations before restarting. If restart is false, this + value is used to allocate memory and later expanded if more memory is required. + - `blocksize::Int = 0`: If blocksize is `> 0`, the solver assumes that the matrix has a + uniformly sized block diagonal structure with square blocks of size `blocksize`. Misusing + this option will lead to incorrect results. + + + If this is set `≤ 0` and during runtime we get a Block Diagonal Matrix, then we will + check if the specialized dispatch can be used. !!! warning @@ -38,13 +40,13 @@ struct SimpleGMRES{UBD} <: AbstractKrylovSubspaceMethod warm_start::Bool function SimpleGMRES{UBD}(; restart::Bool = true, blocksize::Int = 0, - warm_start::Bool = false, memory::Int = 20) where {UBD} + warm_start::Bool = false, memory::Int = 20) where {UBD} UBD && @assert blocksize > 0 return new{UBD}(restart, memory, blocksize, warm_start) end function SimpleGMRES(; restart::Bool = true, blocksize::Int = 0, - warm_start::Bool = false, memory::Int = 20) + warm_start::Bool = false, memory::Int = 20) return SimpleGMRES{blocksize > 0}(; restart, memory, blocksize, warm_start) end @@ -159,7 +161,7 @@ function init_cacheval(alg::SimpleGMRES{UDB}, args...; kwargs...) where {UDB} end function _init_cacheval(::Val{false}, alg::SimpleGMRES, A, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, ::Bool, ::OperatorAssumptions; zeroinit = true, kwargs...) + abstol, reltol, ::Bool, ::OperatorAssumptions; zeroinit = true, kwargs...) @unpack memory, restart, blocksize, warm_start = alg if zeroinit @@ -210,7 +212,8 @@ function _init_cacheval(::Val{false}, alg::SimpleGMRES, A, b, u, Pl, Pr, maxiter rNorm = β ε = abstol + reltol * rNorm - return SimpleGMRESCache{false}(memory, n, restart, maxiters, blocksize, ε, PlisI, PrisI, + return SimpleGMRESCache{false}( + memory, n, restart, maxiters, blocksize, ε, PlisI, PrisI, Pl, Pr, Δx, q, p, x, A, b, abstol, reltol, w, V, s, c, z, R, β, warm_start) end @@ -304,7 +307,8 @@ function SciMLBase.solve!(cache::SimpleGMRESCache{false}, lincache::LinearCache) # Compute and apply current Givens reflection Ωₖ. # [cₖ sₖ] [ r̄ₖ.ₖ ] = [rₖ.ₖ] # [s̄ₖ -cₖ] [hₖ₊₁.ₖ] [ 0 ] - (c[inner_iter], s[inner_iter], R[nr + inner_iter]) = _sym_givens(R[nr + inner_iter], + (c[inner_iter], s[inner_iter], R[nr + inner_iter]) = _sym_givens( + R[nr + inner_iter], Hbis) # Update zₖ = (Qₖ)ᴴβe₁ @@ -387,8 +391,8 @@ function SciMLBase.solve!(cache::SimpleGMRESCache{false}, lincache::LinearCache) end function _init_cacheval(::Val{true}, alg::SimpleGMRES, A, b, u, Pl, Pr, maxiters::Int, - abstol, reltol, ::Bool, ::OperatorAssumptions; zeroinit = true, - blocksize = alg.blocksize) + abstol, reltol, ::Bool, ::OperatorAssumptions; zeroinit = true, + blocksize = alg.blocksize) @unpack memory, restart, warm_start = alg if zeroinit diff --git a/src/simplelu.jl b/src/simplelu.jl index 8049b243d..9c1ad0bf7 100644 --- a/src/simplelu.jl +++ b/src/simplelu.jl @@ -120,7 +120,7 @@ A simple LU-factorization implementation without BLAS. Fast for small matrices. ## Positional Arguments -* pivot::Bool: whether to perform pivoting. Defaults to `true` + - pivot::Bool: whether to perform pivoting. Defaults to `true` """ struct SimpleLUFactorization <: AbstractFactorization pivot::Bool @@ -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 e26c7c621..5c74199cb 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 c28113654..cb64a1246 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -218,13 +218,12 @@ end end end - test_algs = [ LUFactorization(), QRFactorization(), SVDFactorization(), RFLUFactorization(), - LinearSolve.defaultalg(prob1.A, prob1.b), + LinearSolve.defaultalg(prob1.A, prob1.b) ] if VERSION >= v"1.9" && LinearSolve.usemkl @@ -289,9 +288,7 @@ end for alg in (("Default", IterativeSolversJL(kwargs...)), ("CG", IterativeSolversJL_CG(kwargs...)), ("GMRES", IterativeSolversJL_GMRES(kwargs...)), - ("IDRS", IterativeSolversJL_IDRS(kwargs...)) - # ("BICGSTAB",IterativeSolversJL_BICGSTAB(kwargs...)), - # ("MINRES",IterativeSolversJL_MINRES(kwargs...)), + ("IDRS", IterativeSolversJL_IDRS(kwargs...)) # ("BICGSTAB",IterativeSolversJL_BICGSTAB(kwargs...)), # ("MINRES",IterativeSolversJL_MINRES(kwargs...)), ) @testset "$(alg[1])" begin test_interface(alg[2], prob1, prob2) @@ -424,7 +421,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 @@ -432,7 +429,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 diff --git a/test/default_algs.jl b/test/default_algs.jl index 2ffc03b65..d2aafba9d 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -5,24 +5,24 @@ prob = LinearProblem(rand(3, 3), rand(3)) solve(prob) if LinearSolve.appleaccelerate_isavailable() - @test LinearSolve.defaultalg(nothing, zeros(50)).alg === - LinearSolve.DefaultAlgorithmChoice.AppleAccelerateLUFactorization + @test LinearSolve.defaultalg(nothing, zeros(50)).alg === + LinearSolve.DefaultAlgorithmChoice.AppleAccelerateLUFactorization else - @test LinearSolve.defaultalg(nothing, zeros(50)).alg === - LinearSolve.DefaultAlgorithmChoice.RFLUFactorization + @test LinearSolve.defaultalg(nothing, zeros(50)).alg === + LinearSolve.DefaultAlgorithmChoice.RFLUFactorization end prob = LinearProblem(rand(50, 50), rand(50)) solve(prob) if LinearSolve.usemkl - @test LinearSolve.defaultalg(nothing, zeros(600)).alg === - LinearSolve.DefaultAlgorithmChoice.MKLLUFactorization + @test LinearSolve.defaultalg(nothing, zeros(600)).alg === + LinearSolve.DefaultAlgorithmChoice.MKLLUFactorization elseif LinearSolve.appleaccelerate_isavailable() - @test LinearSolve.defaultalg(nothing, zeros(600)).alg === - LinearSolve.DefaultAlgorithmChoice.AppleAccelerateLUFactorization + @test LinearSolve.defaultalg(nothing, zeros(600)).alg === + LinearSolve.DefaultAlgorithmChoice.AppleAccelerateLUFactorization else - @test LinearSolve.defaultalg(nothing, zeros(600)).alg === - LinearSolve.DefaultAlgorithmChoice.LUFactorization + @test LinearSolve.defaultalg(nothing, zeros(600)).alg === + LinearSolve.DefaultAlgorithmChoice.LUFactorization end prob = LinearProblem(rand(600, 600), rand(600)) diff --git a/test/enzyme.jl b/test/enzyme.jl index 89903a858..ac552a45a 100644 --- a/test/enzyme.jl +++ b/test/enzyme.jl @@ -22,8 +22,8 @@ f(A, b1) # Uses BLAS Enzyme.autodiff(Reverse, f, Duplicated(copy(A), dA), Duplicated(copy(b1), db1)) -dA2 = ForwardDiff.gradient(x->f(x,eltype(x).(b1)), copy(A)) -db12 = ForwardDiff.gradient(x->f(eltype(x).(A),x), copy(b1)) +dA2 = ForwardDiff.gradient(x -> f(x, eltype(x).(b1)), copy(A)) +db12 = ForwardDiff.gradient(x -> f(eltype(x).(A), x), copy(b1)) @test dA ≈ dA2 @test db1 ≈ db12 @@ -33,13 +33,20 @@ dA = zeros(n, n); b1 = rand(n); db1 = zeros(n); -_ff = (x,y) -> f(x,y; alg = LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization)) +_ff = (x, y) -> f(x, + y; + alg = LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization)) _ff(copy(A), copy(b1)) -Enzyme.autodiff(Reverse, (x,y) -> f(x,y; alg = LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization)), Duplicated(copy(A), dA), Duplicated(copy(b1), db1)) +Enzyme.autodiff(Reverse, + (x, y) -> f(x, + y; + alg = LinearSolve.DefaultLinearSolver(LinearSolve.DefaultAlgorithmChoice.LUFactorization)), + Duplicated(copy(A), dA), + Duplicated(copy(b1), db1)) -dA2 = ForwardDiff.gradient(x->f(x,eltype(x).(b1)), copy(A)) -db12 = ForwardDiff.gradient(x->f(eltype(x).(A),x), copy(b1)) +dA2 = ForwardDiff.gradient(x -> f(x, eltype(x).(b1)), copy(A)) +db12 = ForwardDiff.gradient(x -> f(eltype(x).(A), x), copy(b1)) @test dA ≈ dA2 @test db1 ≈ db12 @@ -51,7 +58,6 @@ b1 = rand(n); db1 = zeros(n); db12 = zeros(n); - # Batch test n = 4 A = rand(n, n); @@ -79,11 +85,12 @@ end y = [0.0] dy1 = [1.0] dy2 = [1.0] -Enzyme.autodiff(Reverse, fbatch, Duplicated(y, dy1), Duplicated(copy(A), dA), Duplicated(copy(b1), db1)) +Enzyme.autodiff( + Reverse, fbatch, Duplicated(y, dy1), Duplicated(copy(A), dA), Duplicated(copy(b1), db1)) -@test y[1] ≈ f(copy(A),b1) -dA_2 = ForwardDiff.gradient(x->f(x,eltype(x).(b1)), copy(A)) -db1_2 = ForwardDiff.gradient(x->f(eltype(x).(A),x), copy(b1)) +@test y[1] ≈ f(copy(A), b1) +dA_2 = ForwardDiff.gradient(x -> f(x, eltype(x).(b1)), copy(A)) +db1_2 = ForwardDiff.gradient(x -> f(eltype(x).(A), x), copy(b1)) @test dA ≈ dA_2 @test db1 ≈ db1_2 @@ -95,7 +102,8 @@ dA .= 0 dA2 .= 0 db1 .= 0 db12 .= 0 -Enzyme.autodiff(Reverse, fbatch, BatchDuplicated(y, (dy1, dy2)), BatchDuplicated(copy(A), (dA, dA2)), BatchDuplicated(copy(b1), (db1, db12))) +Enzyme.autodiff(Reverse, fbatch, BatchDuplicated(y, (dy1, dy2)), + BatchDuplicated(copy(A), (dA, dA2)), BatchDuplicated(copy(b1), (db1, db12))) @test dA ≈ dA_2 @test db1 ≈ db1_2 @@ -119,11 +127,12 @@ b2 = rand(n); db2 = zeros(n); f(A, b1, b2) -Enzyme.autodiff(Reverse, f, Duplicated(copy(A), dA), Duplicated(copy(b1), db1), Duplicated(copy(b2), db2)) +Enzyme.autodiff(Reverse, f, Duplicated(copy(A), dA), + Duplicated(copy(b1), db1), Duplicated(copy(b2), db2)) -dA2 = ForwardDiff.gradient(x->f(x,eltype(x).(b1),eltype(x).(b2)), copy(A)) -db12 = ForwardDiff.gradient(x->f(eltype(x).(A),x,eltype(x).(b2)), copy(b1)) -db22 = ForwardDiff.gradient(x->f(eltype(x).(A),eltype(x).(b1),x), copy(b2)) +dA2 = ForwardDiff.gradient(x -> f(x, eltype(x).(b1), eltype(x).(b2)), copy(A)) +db12 = ForwardDiff.gradient(x -> f(eltype(x).(A), x, eltype(x).(b2)), copy(b1)) +db22 = ForwardDiff.gradient(x -> f(eltype(x).(A), eltype(x).(b1), x), copy(b2)) @test dA ≈ dA2 @test db1 ≈ db12 @@ -142,7 +151,8 @@ f2(A, b1, b2) dA = zeros(n, n); db1 = zeros(n); db2 = zeros(n); -Enzyme.autodiff(Reverse, f2, Duplicated(copy(A), dA), Duplicated(copy(b1), db1), Duplicated(copy(b2), db2)) +Enzyme.autodiff(Reverse, f2, Duplicated(copy(A), dA), + Duplicated(copy(b1), db1), Duplicated(copy(b2), db2)) @test dA ≈ dA2 @test db1 ≈ db12 @@ -169,10 +179,9 @@ A = rand(n, n); dA = zeros(n, n); b1 = rand(n); for alg in ( - LUFactorization(), - RFLUFactorization(), - # KrylovJL_GMRES(), fails - ) + LUFactorization(), + RFLUFactorization() # KrylovJL_GMRES(), fails +) @show alg function fb(b) prob = LinearProblem(A, b) @@ -192,7 +201,7 @@ for alg in ( end |> collect @show en_jac - @test en_jac ≈ fd_jac rtol=1e-4 + @test en_jac≈fd_jac rtol=1e-4 function fA(A) prob = LinearProblem(A, b1) @@ -212,5 +221,5 @@ for alg in ( end |> collect @show en_jac - @test en_jac ≈ fd_jac rtol=1e-4 -end \ No newline at end of file + @test en_jac≈fd_jac rtol=1e-4 +end diff --git a/test/gpu/cuda.jl b/test/gpu/cuda.jl index 6ce8d2916..d8ea3106b 100644 --- a/test/gpu/cuda.jl +++ b/test/gpu/cuda.jl @@ -27,19 +27,19 @@ function test_interface(alg, prob1, prob2) x2 = prob2.u0 y = solve(prob1, alg; cache_kwargs...) - @test CUDA.@allowscalar(Array(A1 * y) ≈ Array(b1)) + @test CUDA.@allowscalar(Array(A1 * y)≈Array(b1)) cache = SciMLBase.init(prob1, alg; cache_kwargs...) # initialize cache solve!(cache) - @test CUDA.@allowscalar(Array(A1 * cache.u) ≈ Array(b1)) + @test CUDA.@allowscalar(Array(A1 * cache.u)≈Array(b1)) cache.A = copy(A2) solve!(cache) - @test CUDA.@allowscalar(Array(A2 * cache.u) ≈ Array(b1)) + @test CUDA.@allowscalar(Array(A2 * cache.u)≈Array(b1)) cache.b = copy(b2) solve!(cache) - @test CUDA.@allowscalar(Array(A2 * cache.u) ≈ Array(b2)) + @test CUDA.@allowscalar(Array(A2 * cache.u)≈Array(b2)) return end @@ -83,7 +83,7 @@ prob1 = LinearProblem(A', b) prob2 = LinearProblem(transpose(A), b) @testset "Adjoint/Transpose Type: $(alg)" for alg in (NormalCholeskyFactorization(), - CholeskyFactorization(), LUFactorization(), QRFactorization(), nothing) + CholeskyFactorization(), LUFactorization(), QRFactorization(), nothing) sol = solve(prob1, alg; alias_A = false) @test norm(A' * sol.u .- b) < 1e-5 diff --git a/test/hypretests.jl b/test/hypretests.jl index 3cceb6765..c41079c1c 100644 --- a/test/hypretests.jl +++ b/test/hypretests.jl @@ -1,7 +1,7 @@ using HYPRE using HYPRE.LibHYPRE: HYPRE_BigInt, - HYPRE_Complex, HYPRE_IJMatrixGetValues, - HYPRE_IJVectorGetValues, HYPRE_Int + HYPRE_Complex, HYPRE_IJMatrixGetValues, + HYPRE_IJVectorGetValues, HYPRE_Int using LinearAlgebra using LinearSolve using MPI diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl index 36c24ece5..6753be21b 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) @@ -119,7 +119,7 @@ iparm = [ (61, 0), (62, 0), (63, 0), - (64, 0), + (64, 0) ] for i in iparm diff --git a/test/resolve.jl b/test/resolve.jl index 9be8b663d..bd2922f54 100644 --- a/test/resolve.jl +++ b/test/resolve.jl @@ -3,14 +3,14 @@ using LinearSolve, LinearAlgebra, SparseArrays, InteractiveUtils, Test for alg in subtypes(LinearSolve.AbstractFactorization) @show alg if !(alg in [ - DiagonalFactorization, - CudaOffloadFactorization, - AppleAccelerateLUFactorization, - MetalLUFactorization, - ]) && + DiagonalFactorization, + CudaOffloadFactorization, + AppleAccelerateLUFactorization, + MetalLUFactorization + ]) && (!(alg == AppleAccelerateLUFactorization) || LinearSolve.appleaccelerate_isavailable()) && - (!(alg == MKLLUFactorization) || LinearSolve.usemkl) + (!(alg == MKLLUFactorization) || LinearSolve.usemkl) A = [1.0 2.0; 3.0 4.0] alg in [KLUFactorization, UMFPACKFactorization, SparspakFactorization] && (A = sparse(A)) @@ -50,19 +50,19 @@ 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 = [1.0 2.0 - 2.0 1.0] + 2.0 1.0] A = Symmetric(A * A') b = [1.0, 2.0] prob = LinearProblem(A, b) @@ -70,7 +70,7 @@ linsolve = init(prob, CholeskyFactorization(), alias_A = false, alias_b = false) @test solve!(linsolve).u ≈ [-1 / 3, 2 / 3] @test solve!(linsolve).u ≈ [-1 / 3, 2 / 3] A = [1.0 2.0 - 2.0 1.0] + 2.0 1.0] A = Symmetric(A * A') b = [1.0, 2.0] @test solve!(linsolve).u ≈ [-1 / 3, 2 / 3] diff --git a/test/sparse_vector.jl b/test/sparse_vector.jl index 7f80ae1a5..a7ace0202 100644 --- a/test/sparse_vector.jl +++ b/test/sparse_vector.jl @@ -12,7 +12,7 @@ function hess_sparse(x::Vector{T}) where {T} 1.0, 1.0, 12.0 * x[5]^2 + 1.0, - 1.0, + 1.0 ] end rowval = [1, 1, 2, 2, 3, 4, 5, 6] @@ -31,7 +31,7 @@ x0 = [ 1.023999999999997e-7, -1.0, 0.33141395338218227, - -1.0, + -1.0 ] n = length(x0) hess_mat = sparse(rowval, colval, hess_sparse(x0), n, n) From 0f63778449e33e79b4a827e3feb095866ff4e452 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Sun, 7 Jan 2024 15:47:16 -0500 Subject: [PATCH 52/55] Fix UMFPACK to use check=false --- src/factorization.jl | 17 +++++++++++------ test/default_algs.jl | 7 +++++++ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/src/factorization.jl b/src/factorization.jl index 0ef7b754e..efcaaaf00 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -780,21 +780,26 @@ function SciMLBase.solve!(cache::LinearCache, alg::UMFPACKFactorization; kwargs. SparseArrays.decrement(SparseArrays.getrowval(A)) == cacheval.rowval) fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + nonzeros(A)), check=false) else fact = lu!(cacheval, SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), - nonzeros(A))) + nonzeros(A)), check=false) end else - fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A))) + fact = lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)), check=false) end cache.cacheval = fact cache.isfresh = false end - y = ldiv!(cache.u, @get_cacheval(cache, :UMFPACKFactorization), cache.b) - SciMLBase.build_linear_solution(alg, y, nothing, cache) + F = @get_cacheval(cache, :UMFPACKFactorization) + if F.status == SparseArrays.UMFPACK.UMFPACK_OK + y = ldiv!(cache.u, F, cache.b) + SciMLBase.build_linear_solution(alg, y, nothing, cache) + else + SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode=ReturnCode.Infeasible) + end end """ @@ -840,10 +845,10 @@ function init_cacheval(alg::KLUFactorization, A::AbstractSparseArray, b, u, Pl, nonzeros(A))) end +# TODO: guard this against errors function SciMLBase.solve!(cache::LinearCache, alg::KLUFactorization; kwargs...) A = cache.A A = convert(AbstractMatrix, A) - if cache.isfresh cacheval = @get_cacheval(cache, :KLUFactorization) if cacheval !== nothing && alg.reuse_symbolic diff --git a/test/default_algs.jl b/test/default_algs.jl index d2aafba9d..20fba16b9 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -35,6 +35,13 @@ solve(prob) LinearSolve.OperatorAssumptions(false)).alg === LinearSolve.DefaultAlgorithmChoice.QRFactorization + +A = spzeros(100, 100) +A[1,1]=1 +# test that solving a singluar problem doesn't error +prob = LinearProblem(A, ones(100)) +solve(prob) + @test LinearSolve.defaultalg(sprand(10^4, 10^4, 1e-5) + I, zeros(1000)).alg === LinearSolve.DefaultAlgorithmChoice.KLUFactorization prob = LinearProblem(sprand(1000, 1000, 0.5), zeros(1000)) From 7452db86d707d9a7dc6fdd10f89d6e1ba2281fbf Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 23 Feb 2024 16:41:28 -0500 Subject: [PATCH 53/55] typo --- test/default_algs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/default_algs.jl b/test/default_algs.jl index 20fba16b9..da998fad1 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -38,7 +38,7 @@ solve(prob) A = spzeros(100, 100) A[1,1]=1 -# test that solving a singluar problem doesn't error +# test that solving a singular problem doesn't error prob = LinearProblem(A, ones(100)) solve(prob) From 86e7191bc84891936a0d84b5983c4e0aef1cc94d Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 23 Feb 2024 17:07:07 -0500 Subject: [PATCH 54/55] update tests --- test/default_algs.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/default_algs.jl b/test/default_algs.jl index da998fad1..c604129c9 100644 --- a/test/default_algs.jl +++ b/test/default_algs.jl @@ -36,11 +36,12 @@ solve(prob) LinearSolve.DefaultAlgorithmChoice.QRFactorization -A = spzeros(100, 100) -A[1,1]=1 +A = spzeros(2, 2) # test that solving a singular problem doesn't error -prob = LinearProblem(A, ones(100)) -solve(prob) +prob = LinearProblem(A, ones(2)) +@test solve(prob, UMFPACKFactorization()).retcode == ReturnCode.Infeasible +@test_broken solve(prob, KLUFactorization()).retcode == ReturnCode.Infeasible + @test LinearSolve.defaultalg(sprand(10^4, 10^4, 1e-5) + I, zeros(1000)).alg === LinearSolve.DefaultAlgorithmChoice.KLUFactorization From a206054e220801e4accaaa9a86aaed95c2aa2010 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 24 Feb 2024 02:41:08 -0500 Subject: [PATCH 55/55] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2b3422754..4f9505afc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LinearSolve" uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" authors = ["SciML"] -version = "2.23.4" +version = "2.24.0" [deps] ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"