Finish and test well-conditioned form

Author: ChrisRackauckas

src/LinearSolve.jl

@@ -106,6 +106,7 @@ end
 
 export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
        GenericLUFactorization, SimpleLUFactorization, RFLUFactorization,
+       NormalCholeskyFactorization, NormalBunchKaufmanFactorization,
        UMFPACKFactorization, KLUFactorization, FastLUFactorization, FastQRFactorization,
        SparspakFactorization, DiagonalFactorization
 
@@ -119,4 +120,6 @@ export KrylovJL, KrylovJL_CG, KrylovJL_MINRES, KrylovJL_GMRES,
 
 export HYPREAlgorithm
 
+export OperatorAssumptions, OperatorCondition
+
 end

src/common.jl

@@ -63,7 +63,7 @@ __issquare(::OperatorAssumptions{issq,condition}) where {issq,condition} = issq
 __conditioning(::OperatorAssumptions{issq,condition}) where {issq,condition} = condition
 
 
-struct LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq}
+struct LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq, condition}
     A::TA
     b::Tb
     u::Tu
@@ -77,7 +77,7 @@ struct LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq}
     reltol::Ttol
     maxiters::Int
     verbose::Bool
-    assumptions::OperatorAssumptions{issq}
+    assumptions::OperatorAssumptions{issq,condition}
 end
 
 """
@@ -143,9 +143,12 @@ default_tol(::Type{<:Rational}) = 0
 default_tol(::Type{<:Integer}) = 0
 default_tol(::Type{Any}) = 0
 
-function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorithm, Nothing},
+default_alias_A(::Any,::Any,::Any) = false
+default_alias_b(::Any,::Any,::Any) = false
+
+function SciMLBase.init(prob::LinearProblem, alg::SciMLLinearSolveAlgorithm,
                         args...;
-                        alias_A = false, alias_b = false,
+                        alias_A = default_alias_A(alg, prob.A, prob.b), alias_b = default_alias_b(alg, prob.A, prob.b),
                         abstol = default_tol(eltype(prob.A)),
                         reltol = default_tol(eltype(prob.A)),
                         maxiters::Int = length(prob.b),
@@ -187,7 +190,8 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
                         typeof(Pl),
                         typeof(Pr),
                         typeof(reltol),
-                        __issquare(assumptions)
+                        __issquare(assumptions),
+                        __conditioning(assumptions)
                         }(A,
                           b,
                           u0,

src/default.jl

@@ -220,6 +220,10 @@ function defaultalg(A, b, assump::OperatorAssumptions{true})
     alg
 end
 
+function defaultalg(A, b, ::OperatorAssumptions{false,OperatorCondition.WellConditioned})
+    NormalCholeskyFactorization()
+end
+
 function defaultalg(A, b, ::OperatorAssumptions{false,OperatorCondition.IllConditioned})
     QRFactorization()
 end
@@ -234,6 +238,14 @@ end
 
 ## Catch high level interface
 
+function SciMLBase.init(prob::LinearProblem, alg::Nothing,
+                        args...;
+                        assumptions = OperatorAssumptions(Val(issquare(prob.A))),
+                        kwargs...)
+    alg = defaultalg(prob.A, prob.b, assumptions)
+    SciMLBase.init(prob, alg, args...; assumptions, kwargs...)
+end
+
 function SciMLBase.solve(cache::LinearCache, alg::Nothing,
                          args...; assumptions::OperatorAssumptions = OperatorAssumptions(),
                          kwargs...)

src/factorization.jl

@@ -323,8 +323,7 @@ function init_cacheval(alg::GenericFactorization,
 end
 
 # Cholesky needs the posdef matrix, for GenericFactorization assume structure is needed
-function init_cacheval(alg::Union{GenericFactorization,
-                                  GenericFactorization{typeof(cholesky)},
+function init_cacheval(alg::Union{GenericFactorization{typeof(cholesky)},
                                   GenericFactorization{typeof(cholesky!)}}, A, b, u, Pl, Pr,
                        maxiters::Int, abstol, reltol, verbose::Bool,
                        assumptions::OperatorAssumptions)
@@ -476,20 +475,22 @@ end
 
 struct NormalCholeskyFactorization{P} <: AbstractFactorization
     pivot::P
-    perm::Bool
 end
 
-function NormalCholeskyFactorization(; pivot = nothing, perm = nothing)
+function NormalCholeskyFactorization(; pivot = nothing)
     if pivot === nothing
-        @static if VERSION < v"1.7beta"
+        pivot = @static if VERSION < v"1.7beta"
             Val(true)
         else
             RowMaximum()
         end
     end
-    NormalCholeskyFactorization(pivot, perm)
+    NormalCholeskyFactorization(pivot)
 end
 
+default_alias_A(::NormalCholeskyFactorization,::Any,::Any) = true
+default_alias_b(::NormalCholeskyFactorization,::Any,::Any) = true
+
 function init_cacheval(alg::NormalCholeskyFactorization, A, b, u, Pl, Pr,
                        maxiters::Int, abstol, reltol, verbose::Bool,
                        assumptions::OperatorAssumptions)
@@ -500,9 +501,48 @@ function SciMLBase.solve(cache::LinearCache, alg::NormalCholeskyFactorization;
                          kwargs...)
     A = cache.A
     A = convert(AbstractMatrix, A)
-    fact, ipiv = cache.cacheval
     if cache.isfresh
-        fact = cholesky(Symmetric((A)' * A), alg.pivot)
+        if A isa SparseMatrixCSC
+            fact = cholesky(Symmetric((A)' * A))
+        else
+            fact = cholesky(Symmetric((A)' * A), alg.pivot)
+        end
+        cache = set_cacheval(cache, fact)
+    end
+    if A isa SparseMatrixCSC
+        cache.u .= cache.cacheval \ (A' * cache.b)
+        y = cache.u
+    else
+        y = ldiv!(cache.u, cache.cacheval, A' * cache.b)
+    end
+    SciMLBase.build_linear_solution(alg, y, nothing, cache)
+end
+
+## NormalBunchKaufmanFactorization
+
+struct NormalBunchKaufmanFactorization <: AbstractFactorization
+    rook::Bool
+end
+
+function NormalBunchKaufmanFactorization(; rook = false)
+    NormalBunchKaufmanFactorization(rook)
+end
+
+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)
+    ArrayInterface.bunchkaufman_instance(convert(AbstractMatrix, A))
+end
+
+function SciMLBase.solve(cache::LinearCache, alg::NormalBunchKaufmanFactorization;
+                         kwargs...)
+    A = cache.A
+    A = convert(AbstractMatrix, A)
+    if cache.isfresh
+        fact = bunchkaufman(Symmetric((A)' * A), alg.rook)
         cache = set_cacheval(cache, fact)
     end
     y = ldiv!(cache.u, cache.cacheval, A' * cache.b)

src/iterative_wrappers.jl

@@ -15,6 +15,9 @@ function KrylovJL(args...; KrylovAlg = Krylov.gmres!,
                     args, kwargs)
 end
 
+default_alias_A(::KrylovJL,::Any,::Any) = true
+default_alias_b(::KrylovJL,::Any,::Any) = true
+
 function KrylovJL_CG(args...; kwargs...)
     KrylovJL(args...; KrylovAlg = Krylov.cg!, kwargs...)
 end
@@ -205,6 +208,9 @@ function IterativeSolversJL(args...;
                               args, kwargs)
 end
 
+default_alias_A(::IterativeSolversJL,::Any,::Any) = true
+default_alias_b(::IterativeSolversJL,::Any,::Any) = true
+
 function IterativeSolversJL_CG(args...; kwargs...)
     IterativeSolversJL(args...;
                        generate_iterator = IterativeSolvers.cg_iterator!,
@@ -312,6 +318,9 @@ function KrylovKitJL(args...;
     return KrylovKitJL(KrylovAlg, gmres_restart, args, kwargs)
 end
 
+default_alias_A(::KrylovKitJL,::Any,::Any) = true
+default_alias_b(::KrylovKitJL,::Any,::Any) = true
+
 function KrylovKitJL_CG(args...; kwargs...)
     KrylovKitJL(args...; KrylovAlg = KrylovKit.CG, kwargs..., isposdef = true)
 end

test/nonsquare.jl

@@ -29,4 +29,19 @@ A = sprandn(1000, 100, 0.1)
 b = randn(1001)
 prob = LinearProblem(A, view(b, 1:1000))
 linsolve = init(prob, QRFactorization())
-solve(linsolve)
+solve(linsolve)
+
+A = randn(1000, 100)
+b = randn(1000)
+@test isapprox(solve(LinearProblem(A, b)).u, Symmetric(A' * A) \ (A' * b))
+solve(LinearProblem(A, b)).u;
+solve(LinearProblem(A, b), (LinearSolve.NormalCholeskyFactorization())).u;
+solve(LinearProblem(A, b), (LinearSolve.NormalBunchKaufmanFactorization())).u;
+solve(LinearProblem(A, b), assumptions = (OperatorAssumptions(false; condition = OperatorCondition.WellConditioned))).u;
+
+A = sprandn(5000, 100, 0.1)
+b = randn(5000)
+@test isapprox(solve(LinearProblem(A, b)).u, ldlt(A' * A) \ (A' * b))
+solve(LinearProblem(A, b)).u;
+solve(LinearProblem(A, b), (LinearSolve.NormalCholeskyFactorization())).u;
+solve(LinearProblem(A, b), assumptions = (OperatorAssumptions(false; condition = OperatorCondition.WellConditioned))).u;