Support new precs API of LinearSolve

.github/workflows/ci.yml

@@ -15,9 +15,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - 'lts' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1'   # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
-#          - 'nightly'
+          - 'nightly'
         os:
           - ubuntu-latest
           - windows-latest
@@ -34,7 +34,7 @@ jobs:
             arch: x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
@@ -59,7 +59,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
       - run: |

CHANGELOG.md

@@ -1,5 +1,16 @@
 # Changelog
 
+## [2.0.0] - Planned
+
+### Breaking
+- remove solver + precon API which is not based on precs or directly overloading `\`.
+  Fully rely on LinearSolve (besides `\`)
+- Move AMGBuilder, ILUZeroBuilder etc. to the correspondig packages (depending on the PRs)
+- remove "old" SparseMatrixLNK (need to benchmark before)
+
+## [1.6.0] - WIP
+- Support precs API of LinearSolve.jl
+
 ## [1.5.3] - 2024-10-07
 - Moved repo to WIAS-PDELib organization
 

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <[email protected]>"]
-version = "1.5.3"
+version = "1.6.0"
 
 [deps]
 AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288"
@@ -19,20 +19,23 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288"
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 
 [extensions]
 ExtendableSparseAMGCLWrapExt = "AMGCLWrap"
 ExtendableSparseAlgebraicMultigridExt = "AlgebraicMultigrid"
 ExtendableSparseIncompleteLUExt = "IncompleteLU"
 ExtendableSparsePardisoExt = "Pardiso"
+ExtendableSparseLinearSolveExt = "LinearSolve"
 
 [compat]
-AMGCLWrap = "0.4,1"
+AMGCLWrap = "2"
 AlgebraicMultigrid = "0.4,0.5,0.6"
 DocStringExtensions = "0.8, 0.9"
 ILUZero = "0.2"
 IncompleteLU = "^0.2.1"
+LinearSolve = "2.36.0"
 Pardiso = "0.5.1"
 Sparspak = "0.3.6"
 StaticArrays = "1.5.24"

docs/Project.toml

@@ -9,8 +9,10 @@ ILUZero = "88f59080-6952-5380-9ea5-54057fb9a43f"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
 
 [compat]
 Documenter = "1.0"
 IterativeSolvers = "0.9"
+LinearSolve = "2.36.0"

docs/make.jl

@@ -1,4 +1,3 @@
-push!(LOAD_PATH, "../src/")
 using Documenter, ExtendableSparse, AlgebraicMultigrid, IncompleteLU, Sparspak, LinearAlgebra
 
 function mkdocs()
@@ -12,10 +11,10 @@ function mkdocs()
              pages = [
                  "Home" => "index.md",
                  "example.md",
-                 "linearsolve.md",
                  "extsparse.md",
-                 "iter.md",
+                 "linearsolve.md",
                  "internal.md",
+                 "iter.md",
                  "changes.md",
              ])
 end

docs/src/extsparse.md

@@ -2,6 +2,10 @@
 
 ## Matrix creation and update API
 
+```@docs
+ExtendableSparseMatrix
+```
+
 ```@autodocs
 Modules = [ExtendableSparse]
 Pages = ["extendable.jl"]

docs/src/internal.md

@@ -15,7 +15,7 @@ Pages = ["sparsematrixcsc.jl"]
 ```
 ## New API 
 Under development - aimed at multithreading
-```@autodocs
+```@autodocs; canonical = false
 Modules = [ExtendableSparse]
 Pages = ["abstractsparsematrixextension.jl",
     "abstractextendablesparsematrixcsc.jl",

docs/src/iter.md

@@ -1,6 +1,8 @@
 # Factorizations & Preconditioners
 
-This functionality probably will be reduced in favor of LinearSolve.jl.
+This functionality is deprecated as of version 1.6 an will be removed
+with  version 2.0. Use the  [integration with LinearSolve.jl](/linearsolve/#Integration-with-LinearSolve.jl).
+
 
 ## Factorizations
 

docs/src/linearsolve.md

@@ -1,63 +1,152 @@
-# Integration with LinearSolve.jl
+# Linear System solution
+
+## The `\` operator
+The packages overloads `\` for the ExtendableSparseMatrix type.
+The following example uses [`fdrand`](@ref) to create a test matrix and solve
+the corresponding linear system of equations.
+
+```@example
+using ExtendableSparse
+A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
+x = ones(1000)
+b = A * x
+y = A \ b
+sum(y)
+```
+
+This works as well for number types besides `Float64` and related, in this case,
+by default a LU factorization based on Sparspak ist used.
+
+```@example
+using ExtendableSparse
+using MultiFloats
+A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix)
+x = ones(Float64x2,1000)
+b = A * x
+y = A \ b
+sum(y)
+```
+
+## Solving with LinearSolve.jl
 
 Starting with version 0.9.6, ExtendableSparse is compatible
 with [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl).
 Since version 0.9.7, this is facilitated via the
-AbstractSparseMatrixCSC interface.
+AbstractSparseMatrixCSC interface. 
 
-```@autodocs
-Modules = [ExtendableSparse]
-Pages = ["linearsolve.jl"]
-```
 
-We can create a test problem and solve it with the `\` operator.
+The same problem can be solved via `LinearSolve.jl`:
 
 ```@example
-using ExtendableSparse # hide
+using ExtendableSparse
+using LinearSolve
 A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
 x = ones(1000)
 b = A * x
-y = A \ b
+y = solve(LinearProblem(A, b)).u
 sum(y)
 ```
 
-The same problem can be solved by the tools available via `LinearSolve.jl`:
+```@example
+using ExtendableSparse
+using LinearSolve
+using MultiFloats
+A = fdrand(Float64x2, 10, 10, 10; matrixtype = ExtendableSparseMatrix)
+x = ones(Float64x2,1000)
+b = A * x
+y = solve(LinearProblem(A, b), SparspakFactorization()).u
+sum(y)
+```
+
+## Preconditioned Krylov solvers  with LinearSolve.jl
+
+Since version 1.6, preconditioner constructors can be passed to iterative solvers via the [`precs`
+keyword argument](https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners/#prec)
+to the iterative solver specification.
 
 ```@example
-using ExtendableSparse # hide
-using LinearSolve # hide
+using ExtendableSparse
+using LinearSolve 
+using ExtendableSparse: ILUZeroPreconBuilder
 A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
 x = ones(1000)
 b = A * x
-y = solve(LinearProblem(A, b), SparspakFactorization()).u
+y = LinearSolve.solve(LinearProblem(A, b), 
+                      KrylovJL_CG(; precs=ILUZeroPreconBuilder())).u
 sum(y)
 ```
 
-Also, the iterative method interface works with ExtendableSparse.
+## Available preconditioners
+ExtendableSparse provides constructors for preconditioners wich can be used as `precs`.
+These generally return a tuple `(Pl,I)` of a left preconditioner and a trivial right preconditioner.
+
+ExtendableSparse has a number of package extensions which construct preconditioners
+from some other packages. In the future, these packages may provide this functionality on their own.
+
+```@docs
+ExtendableSparse.ILUZeroPreconBuilder
+ExtendableSparse.ILUTPreconBuilder
+ExtendableSparse.SmoothedAggregationPreconBuilder
+ExtendableSparse.RugeStubenPreconBuilder
+```
+
+In addition, ExtendableSparse implements some preconditioners:
+
+```@docs
+ExtendableSparse.JacobiPreconBuilder
+```
+
+LU factorizations of  matrices from previous iteration steps may be good
+preconditioners for Krylov solvers called during a nonlinear solve via
+Newton's method. For this purpose, ExtendableSparse provides a preconditioner constructor
+which wraps sparse LU factorizations  supported by LinearSolve.jl
+```@docs
+ExtendableSparse.LinearSolvePreconBuilder
+```
+
+
+Block preconditioner constructors are provided as well
+```@docs;  canonical=false
+ExtendableSparse.BlockPreconBuilder
+```
+
 
+The example beloww shows how to create a  block Jacobi preconditioner where the blocks are defined by even and odd
+degrees of freedom, and the diagonal blocks are solved using UMFPACK.
 ```@example
-using ExtendableSparse # hide
-using LinearSolve # hide
-using SparseArrays # hide
-using ILUZero # hide
+using ExtendableSparse
+using LinearSolve 
+using ExtendableSparse: LinearSolvePreconBuilder, BlockPreconBuilder
 A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
 x = ones(1000)
 b = A * x
-y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG();
-                      Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u
+partitioning=A->[1:2:size(A,1), 2:2:size(A,1)]
+umfpackprecs=LinearSolvePreconBuilder(UMFPACKFactorization())
+blockprecs=BlockPreconBuilder(;precs=umfpackprecs, partitioning)
+y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG(; precs=blockprecs)).u
 sum(y)
 ```
+`umpfackpreconbuilder` e.g. could be replaced by `SmoothedAggregationPreconBuilder()`. Moreover, this approach
+works for any `AbstractSparseMatrixCSC`.
+
+
+## Deprecated API
+Passing a preconditioner via the `Pl` or `Pr` keyword arguments
+will be deprecated in LinearSolve. ExtendableSparse used to
+export a number of wrappers for preconditioners from other packages
+for this purpose. This approach is deprecated as of v1.6 and will be removed
+with v2.0.
 
-However, ExtendableSparse provides a number of wrappers around preconditioners
-from various Julia packages.
 ```@example
-using ExtendableSparse # hide
-using LinearSolve # hide
-using ILUZero # hide
+using ExtendableSparse
+using LinearSolve
+using SparseArray
+using ILUZero
 A = fdrand(10, 10, 10; matrixtype = ExtendableSparseMatrix)
 x = ones(1000)
 b = A * x
 y = LinearSolve.solve(LinearProblem(A, b), KrylovJL_CG();
-                      Pl = ILUZeroPreconditioner(A)).u
+                      Pl = ILUZero.ilu0(SparseMatrixCSC(A))).u
 sum(y)
 ```
+

ext/ExtendableSparseAMGCLWrapExt.jl

@@ -1,3 +1,5 @@
+# The whole extension is deprecated
+# TODO remove in v2.0
 module ExtendableSparseAMGCLWrapExt
 using ExtendableSparse
 using AMGCLWrap
@@ -10,6 +12,8 @@ mutable struct AMGCL_AMGPreconditioner <: AbstractPreconditioner
     factorization::AMGCLWrap.AMGPrecon
     kwargs
     function ExtendableSparse.AMGCL_AMGPreconditioner(; kwargs...)
+        Base.depwarn("AMGCL_AMGPreconditioner() is deprecated. Use LinearSolve with `precs=AMGCLWrap.AMGPreconBuilder()`  instead.",
+                     :AMGCL_AMGPreconditioner)
         precon = new()
         precon.kwargs = kwargs
         precon
@@ -35,6 +39,8 @@ mutable struct AMGCL_RLXPreconditioner <: AbstractPreconditioner
     factorization::AMGCLWrap.RLXPrecon
     kwargs
     function ExtendableSparse.AMGCL_RLXPreconditioner(; kwargs...)
+        Base.depwarn("AMGCL_RLXPreconditioner() is deprecated. Use LinearSolve with  `precs=AMGCLWrap.RLXPreconBuilder()`  instead.",
+                     :AMGCL_RLXPreconditioner)
         precon = new()
         precon.kwargs = kwargs
         precon