diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index c99162a..fb88e6a 100644
--- a/.github/workflows/ci.yml
+++ b/.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: |
diff --git a/CHANGELOG.md b/CHANGELOG.md
index fc41761..16a6b8a 100644
--- a/CHANGELOG.md
+++ b/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
 
diff --git a/Project.toml b/Project.toml
index b152216..eabf350 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "1.5.3"
+version = "1.6.0"
 
 [deps]
 AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288"
@@ -19,6 +19,7 @@ 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]
@@ -26,13 +27,15 @@ 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"
diff --git a/docs/Project.toml b/docs/Project.toml
index 5106403..5a3ae26 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/docs/make.jl b/docs/make.jl
index 2eb1e8d..c9d39b1 100644
--- a/docs/make.jl
+++ b/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
diff --git a/docs/src/extsparse.md b/docs/src/extsparse.md
index c4d2cd6..e756042 100644
--- a/docs/src/extsparse.md
+++ b/docs/src/extsparse.md
@@ -2,6 +2,10 @@
 
 ## Matrix creation and update API
 
+```@docs
+ExtendableSparseMatrix
+```
+
 ```@autodocs
 Modules = [ExtendableSparse]
 Pages = ["extendable.jl"]
diff --git a/docs/src/internal.md b/docs/src/internal.md
index 71a10f3..477058c 100644
--- a/docs/src/internal.md
+++ b/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",
diff --git a/docs/src/iter.md b/docs/src/iter.md
index 22b68d0..31030c0 100644
--- a/docs/src/iter.md
+++ b/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
 
diff --git a/docs/src/linearsolve.md b/docs/src/linearsolve.md
index 22d87f6..ffaa9de 100644
--- a/docs/src/linearsolve.md
+++ b/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)
 ```
+
diff --git a/ext/ExtendableSparseAMGCLWrapExt.jl b/ext/ExtendableSparseAMGCLWrapExt.jl
index 3389ce8..f078a0e 100644
--- a/ext/ExtendableSparseAMGCLWrapExt.jl
+++ b/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
diff --git a/ext/ExtendableSparseAlgebraicMultigridExt.jl b/ext/ExtendableSparseAlgebraicMultigridExt.jl
index a21b935..7abdbc2 100644
--- a/ext/ExtendableSparseAlgebraicMultigridExt.jl
+++ b/ext/ExtendableSparseAlgebraicMultigridExt.jl
@@ -1,16 +1,37 @@
 module ExtendableSparseAlgebraicMultigridExt
 using ExtendableSparse
-using AlgebraicMultigrid
+using AlgebraicMultigrid: AlgebraicMultigrid, ruge_stuben, smoothed_aggregation, aspreconditioner
+using SparseArrays: SparseMatrixCSC, AbstractSparseMatrixCSC
+using LinearAlgebra: I
+
+
+import ExtendableSparse: SmoothedAggregationPreconBuilder
+import ExtendableSparse: RugeStubenPreconBuilder
+
+(b::SmoothedAggregationPreconBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I)
+(b::RugeStubenPreconBuilder)(A::AbstractSparseMatrixCSC,p)= (aspreconditioner(ruge_stuben(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)),I)
+
+
+####
+# Deprecated from here on
+# TODO remove in v2.0
 
 import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update!
 
 ######################################################################################
+rswarned=false
+
 mutable struct RS_AMGPreconditioner <: AbstractPreconditioner
     A::ExtendableSparseMatrix
     factorization::AlgebraicMultigrid.Preconditioner
     kwargs
     blocksize
     function ExtendableSparse.RS_AMGPreconditioner(blocksize=1; kwargs...)
+        global rswarned
+        if !rswarned
+            @warn "RS_AMGPreconditioner is deprecated. Use LinearSolve with `precs=RugeStubenPreconBuilder()` instead"
+            rswarned=true
+        end
         precon = new()
         precon.kwargs = kwargs
         precon.blocksize=blocksize
@@ -31,13 +52,20 @@ end
 allow_views(::RS_AMGPreconditioner)=true
 allow_views(::Type{RS_AMGPreconditioner})=true
 
+
 ######################################################################################
+sawarned=false
 mutable struct SA_AMGPreconditioner <: AbstractPreconditioner
     A::ExtendableSparseMatrix
     factorization::AlgebraicMultigrid.Preconditioner
     kwargs
     blocksize
     function ExtendableSparse.SA_AMGPreconditioner(blocksize=1; kwargs...)
+        global sawarned
+        if !sawarned
+            @warn "SA_AMGPreconditioner is deprecated. Use LinearSolve with `precs=SmoothedAggregationPreconBuilder()` instead"
+            sawarned=true
+        end
         precon = new()
         precon.kwargs = kwargs
         precon.blocksize=blocksize
diff --git a/ext/ExtendableSparseIncompleteLUExt.jl b/ext/ExtendableSparseIncompleteLUExt.jl
index dc6ff25..bfeb7ed 100644
--- a/ext/ExtendableSparseIncompleteLUExt.jl
+++ b/ext/ExtendableSparseIncompleteLUExt.jl
@@ -1,14 +1,29 @@
 module ExtendableSparseIncompleteLUExt
 using ExtendableSparse
 using IncompleteLU 
+using LinearAlgebra: I
+using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC, getcolptr, rowvals, nonzeros
+
+import ExtendableSparse: ILUTPreconBuilder
+
+(b::ILUTPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(IncompleteLU.ilu(SparseMatrixCSC(A); τ = b.droptol),I)
+
 
 import ExtendableSparse: @makefrommatrix, AbstractPreconditioner, update!
 
+
+# Deprecated from here
+warned=false
 mutable struct ILUTPreconditioner <: AbstractPreconditioner
     A::ExtendableSparseMatrix
     factorization::IncompleteLU.ILUFactorization
     droptol::Float64
     function ExtendableSparse.ILUTPreconditioner(; droptol = 1.0e-3)
+        global warned
+        if !warned
+            @warn "ILUTPreconditioner is deprecated. Use LinearSolve with `precs=ILUTPreconBuilder()` instead"
+            warned=true
+        end
         p = new()
         p.droptol = droptol
         p
@@ -25,5 +40,6 @@ function update!(precon::ILUTPreconditioner)
     @inbounds flush!(A)
     precon.factorization = IncompleteLU.ilu(A.cscmatrix; τ = precon.droptol)
 end
+
 end
 
diff --git a/ext/ExtendableSparseLinearSolveExt.jl b/ext/ExtendableSparseLinearSolveExt.jl
new file mode 100644
index 0000000..2c5335a
--- /dev/null
+++ b/ext/ExtendableSparseLinearSolveExt.jl
@@ -0,0 +1,26 @@
+module ExtendableSparseLinearSolveExt
+using LinearSolve
+import ExtendableSparse: LinearSolvePreconBuilder
+import LinearAlgebra
+using SparseArrays: AbstractSparseMatrixCSC
+
+
+struct LinearSolvePrecon{T}
+    cache::T
+end
+
+function LinearSolvePrecon(A,method::LinearSolve.AbstractFactorization)
+    pr = LinearProblem(A, zeros(eltype(A), size(A, 1)))
+    LinearSolvePrecon(init(pr, method))
+end
+
+function LinearAlgebra.ldiv!(u, P::LinearSolvePrecon, b)
+    P.cache.b = b
+    sol = solve!(P.cache)
+    copyto!(u, sol.u)
+end
+
+(b::LinearSolvePreconBuilder)(A::AbstractSparseMatrixCSC,p) = (LinearSolvePrecon(A,b.method), LinearAlgebra.I)
+
+end
+
diff --git a/ext/ExtendableSparsePardisoExt.jl b/ext/ExtendableSparsePardisoExt.jl
index 4faaa2f..69b8e29 100644
--- a/ext/ExtendableSparsePardisoExt.jl
+++ b/ext/ExtendableSparsePardisoExt.jl
@@ -2,7 +2,7 @@ module ExtendableSparsePardisoExt
 using ExtendableSparse
 using LinearAlgebra
 using Pardiso
-
+# Deprecated
 import ExtendableSparse: @makefrommatrix, update!, AbstractLUFactorization
 
 abstract type AbstractPardisoLU <: AbstractLUFactorization end
diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl
index 9beee11..e4a67c8 100644
--- a/src/ExtendableSparse.jl
+++ b/src/ExtendableSparse.jl
@@ -21,6 +21,7 @@ if USE_GPL_LIBS
 end
 
 
+include("compat.jl") # @public
 
 include("matrix/sparsematrixcsc.jl")
 include("matrix/abstractsparsematrixextension.jl")
@@ -31,10 +32,33 @@ include("matrix/extendable.jl")
 include("matrix/genericmtextendablesparsematrixcsc.jl")
 include("matrix/genericextendablesparsematrixcsc.jl")
 
+"""
+    ExtendableSparseMatrix
+
+Aliased to [`ExtendableSparseMatrixCSC`](@ref)
+"""
 const ExtendableSparseMatrix=ExtendableSparseMatrixCSC
+
+
+"""
+    MTExtendableSparseMatrixCSC
+
+Multithreaded extendable sparse matrix  (Experimental).
+
+Aliased to [`GenericMTExtendableSparseMatricCSC`](@ref) with [`SparseMatrixDILNKC`](@ref) 
+scalar matrix parameter.
+"""
 const MTExtendableSparseMatrixCSC{Tv,Ti}=GenericMTExtendableSparseMatrixCSC{SparseMatrixDILNKC{Tv,Ti},Tv,Ti}
 MTExtendableSparseMatrixCSC(m,n,args...)=MTExtendableSparseMatrixCSC{Float64,Int64}(m,n,args...)
 
+"""
+    STExtendableSparseMatrixCSC
+
+Single threaded extendable sparse matrix (Experimental).
+
+Aliased to [`GenericExtendableSparseMatricCSC`](@ref) with [`SparseMatrixDILNKC`](@ref) 
+scalar matrix parameter.
+"""
 const STExtendableSparseMatrixCSC{Tv,Ti}=GenericExtendableSparseMatrixCSC{SparseMatrixDILNKC{Tv,Ti},Tv,Ti}
 STExtendableSparseMatrixCSC(m,n,args...)=STExtendableSparseMatrixCSC{Float64,Int64}(m,n,args...)
 
@@ -52,13 +76,21 @@ include("factorizations/factorizations.jl")
 include("factorizations/simple_iteration.jl")
 export simple, simple!
 
+include("preconbuilders.jl")
+export LinearSolvePreconBuilder, BlockPreconBuilder, JacobiPreconBuilder
+
+@public ILUZeroPreconBuilder, ILUTPreconBuilder, SmoothedAggregationPreconBuilder, RugeStubenPreconBuilder
+
+
 include("matrix/sprand.jl")
 export sprand!, sprand_sdd!, fdrand, fdrand!, fdrand_coo, solverbenchmark
 
 export rawupdateindex!, updateindex!
 
 
-
+#####
+# All of the following is deprecated in favor of the precs bases
+# API
 
 export JacobiPreconditioner,
     ILU0Preconditioner,
@@ -208,4 +240,5 @@ Pardiso.set_iparm!(plu.ps,5,13.0)
 function MKLPardisoLU end
 export MKLPardisoLU
 
+
 end # module
diff --git a/src/compat.jl b/src/compat.jl
new file mode 100644
index 0000000..570f712
--- /dev/null
+++ b/src/compat.jl
@@ -0,0 +1,18 @@
+# Copied from MIT licensed
+# https://github.com/JuliaDiff/DifferentiationInterface.jl/blob/main/DifferentiationInterface/src/compat.jl
+#
+macro public(ex)
+    if VERSION >= v"1.11.0-DEV.469"
+        args = if ex isa Symbol
+            (ex,)
+        elseif Base.isexpr(ex, :tuple)
+            ex.args
+        else
+            error("something informative")
+        end
+        esc(Expr(:public, args...))
+    else
+        nothing
+    end
+end
+
diff --git a/src/factorizations/blockpreconditioner.jl b/src/factorizations/blockpreconditioner.jl
index af13986..e9bc1f8 100644
--- a/src/factorizations/blockpreconditioner.jl
+++ b/src/factorizations/blockpreconditioner.jl
@@ -30,7 +30,7 @@ function BlockPreconditioner end
 
 Factorizations on matrix partitions within a block preconditioner may or may not work with array views.
 E.g. the umfpack factorization cannot work with views, while ILUZeroPreconditioner can.
-Implementing a method for `allow_views` returning `false` resp. `true` allows to dispatch to the proper case.
+ Implementing a method for `allow_views` returning `false` resp. `true` allows to dispatch to the proper case.
 """
 allow_views(::Any)=false
 
@@ -98,3 +98,32 @@ function LinearAlgebra.ldiv!(u,p::BlockPreconditioner,v)
 end
 
 Base.eltype(p::BlockPreconditioner)=eltype(p.facts[1])
+
+
+"""
+     BlockPreconBuilder(;precs=UMFPACKPreconBuilder(),  
+                         partitioning = A -> [1:size(A,1)]
+
+Return callable object constructing a left block Jacobi preconditioner 
+from partition of unknowns.
+
+- `partitioning(A)`  shall return a vector of AbstractVectors describing the indices of the partitions of the matrix. 
+  For a matrix of size `n x n`, e.g. partitioning could be `[ 1:n÷2, (n÷2+1):n]` or [ 1:2:n, 2:2:n].
+
+- `precs(A,p)` shall return a left precondioner for a matrix block.
+"""
+Base.@kwdef mutable struct BlockPreconBuilder
+    precs=UMFPACKPreconBuilder()
+    partitioning= A -> [1:size(A,1)]
+end
+
+function (blockprecs::BlockPreconBuilder)(A,p)
+    (;precs, partitioning)=blockprecs
+    factorization= A->precs(A,p)[1]
+    bp=BlockPreconditioner(A;partitioning=partitioning(A), factorization)
+    (bp,LinearAlgebra.I)
+end
+
+"""
+    Allow array for precs => different precoms
+"""
diff --git a/src/factorizations/cholmod_cholesky.jl b/src/factorizations/cholmod_cholesky.jl
index 1e3fbc8..66228b4 100644
--- a/src/factorizations/cholmod_cholesky.jl
+++ b/src/factorizations/cholmod_cholesky.jl
@@ -5,6 +5,7 @@ mutable struct CholeskyFactorization <: AbstractLUFactorization
     A64::Any
 end
 
+cholwarned=false
 """
 CholeskyFactorization(;valuetype=Float64, indextype=Int64)
 CholeskyFactorization(matrix)
@@ -12,6 +13,11 @@ CholeskyFactorization(matrix)
 Default Cholesky factorization via cholmod.
 """
 function CholeskyFactorization()
+    global cholwarned
+    if !cholwarned
+        @warn "ExtendableSparse.CholeskyFactorization is deprecated. Use LinearSolve.CholeskyFactorization` instead"
+        cholwarned=true
+    end
     CholeskyFactorization(nothing, nothing, 0, nothing)
 end
 
diff --git a/src/factorizations/iluzero.jl b/src/factorizations/iluzero.jl
index 2f068c5..efc7ef7 100644
--- a/src/factorizations/iluzero.jl
+++ b/src/factorizations/iluzero.jl
@@ -1,8 +1,14 @@
+iluzerowarned=false
 mutable struct ILUZeroPreconditioner <: AbstractPreconditioner
     A::ExtendableSparseMatrix
     factorization::ILUZero.ILU0Precon
     phash::UInt64
     function ILUZeroPreconditioner()
+        global iluzerowarned
+        if !iluzerowarned
+            @warn "ILUZeroPreconditioner is deprecated. Use LinearSolve with `precs=ILUZeroPreconBuilder()` instead"
+            iluzerowarned=true
+        end
         p = new()
         p.phash = 0
         p
@@ -35,12 +41,18 @@ allow_views(::Type{ILUZeroPreconditioner})=true
 
 
 
+biluzerowarned=false
 mutable struct PointBlockILUZeroPreconditioner <: AbstractPreconditioner
     A::ExtendableSparseMatrix
     factorization::ILUZero.ILU0Precon
     phash::UInt64
     blocksize::Int
     function PointBlockILUZeroPreconditioner(;blocksize=1)
+        global biluzerowarned
+        if !biluzerowarned
+            @warn "PointBlockILUZeroPreconditioner is deprecated. Use LinearSolve with `precs=ILUZeroPreconBuilder(; blocksize=$(blocksize))` instead"
+            biluzerowarned=true
+        end
         p = new()
         p.phash = 0
         p.blocksize=blocksize
diff --git a/src/factorizations/linearsolve.jl b/src/factorizations/linearsolve.jl
deleted file mode 100644
index 81bfe57..0000000
--- a/src/factorizations/linearsolve.jl
+++ /dev/null
@@ -1,17 +0,0 @@
-"""
-    $(SIGNATURES)
-
-Create LinearProblem from ExtendableSparseMatrix.
-"""
-function LinearSolve.LinearProblem(A::ExtendableSparseMatrix, b::AbstractArray)
-    LinearSolve.LinearProblem(SparseMatrixCSC(A), b)
-end
-
-"""
-    $(SIGNATURES)
-
-Update LinearSolve cache.
-"""
-function LinearSolve.set_A(cache::LinearSolve.LinearCache, A::ExtendableSparseMatrix)
-    LinearSolve.set_A(cache, SparseMatrixCSC(A))
-end
diff --git a/src/preconbuilders.jl b/src/preconbuilders.jl
new file mode 100644
index 0000000..9cd5782
--- /dev/null
+++ b/src/preconbuilders.jl
@@ -0,0 +1,105 @@
+"""
+        LinearSolvePreconBuilder(; method=UMFPACKFactorization())
+
+Return callable object constructing a formal left preconditioner from a sparse LU factorization using LinearSolve
+as the `precs` parameter for a  [`BlockPreconBuilder`](@ref)  or  iterative methods wrapped by LinearSolve.jl.
+"""
+Base.@kwdef struct LinearSolvePreconBuilder
+    method=UMFPACKFactorization()
+end
+(::LinearSolvePreconBuilder)(A,p)= error("import LinearSolve in order to use LinearSolvePreconBuilder")
+
+
+"""
+    JacobiPreconBuilder()
+
+Return callable object constructing a left Jacobi preconditioner
+to be passed as the `precs` parameter to iterative methods wrapped by LinearSolve.jl.
+"""
+struct JacobiPreconBuilder end
+(::JacobiPreconBuilder)(A::AbstractSparseMatrixCSC,p)=(JacobiPreconditioner(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),nonzeros(A))),LinearAlgebra.I)
+
+
+"""
+    ILUZeroPreconBuilder(;blocksize=1)
+
+Return callable object constructing a left zero fill-in ILU preconditioner 
+using [ILUZero.jl](https://github.com/mcovalt/ILUZero.jl)
+"""
+Base.@kwdef mutable struct ILUZeroPreconBuilder
+    blocksize::Int = 1
+end
+
+struct ILUBlockPrecon{N,NN,Tv,Ti}
+    ilu0::ILUZero.ILU0Precon{SMatrix{N, N, Tv, NN}, Ti, SVector{N, Tv}}
+end
+
+function LinearAlgebra.ldiv!(Y::Vector{Tv},
+                             A::ILUBlockPrecon{N,NN,Tv,Ti},
+                             B::Vector{Tv}) where {N,NN,Tv,Ti}
+    BY=reinterpret(SVector{N,Tv},Y)
+    BB=reinterpret(SVector{N,Tv},B)
+    ldiv!(BY,A.ilu0,BB)
+    Y
+end
+
+function (b::ILUZeroPreconBuilder)(A0,p)
+    A=SparseMatrixCSC(size(A0)..., getcolptr(A0), rowvals(A0),nonzeros(A0))
+    if b.blocksize==1
+        (ILUZero.ilu0(A),LinearAlgebra.I)
+    else
+        (ILUBlockPrecon(ILUZero.ilu0(pointblock(A,b.blocksize),SVector{b.blocksize,eltype(A)})),LinearAlgebra.I)
+    end
+end
+
+
+"""
+    ILUTPreconBuilder(; droptol=0.1)
+
+Return callable object constructing a left ILUT preconditioner 
+using [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl)
+"""
+Base.@kwdef struct ILUTPreconBuilder
+    droptol::Float64=0.1
+end
+(::ILUTPreconBuilder)(A,p)= error("import IncompleteLU.jl in order to use ILUTBuilder")
+
+
+
+"""
+    SmoothedAggregationPreconBuilder(;blocksize=1, kwargs...)
+
+Return callable object constructing a left smoothed aggregation algebraic multigrid preconditioner
+using [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl).
+
+Needs `import AlgebraicMultigrid` to trigger the corresponding extension.
+"""
+struct SmoothedAggregationPreconBuilder{Tk}
+    blocksize::Int
+    kwargs::Tk
+end
+
+function SmoothedAggregationPreconBuilder(;blocksize=1, kwargs...)
+    SmoothedAggregationPreconBuilder(blocksize,kwargs)
+end
+
+(::SmoothedAggregationPreconBuilder)(A,p)= error("import AlgebraicMultigrid in order to use SmoothedAggregationPreconBuilder")
+
+"""
+   RugeStubenPreconBuilder(;blocksize=1, kwargs...)
+
+Return callable object constructing a  left algebraic multigrid preconditioner after Ruge & Stüben
+using [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl).
+
+Needs `import AlgebraicMultigrid` to trigger the corresponding extension.
+"""
+struct RugeStubenPreconBuilder{Tk}
+    blocksize::Int
+    kwargs::Tk
+end
+
+function RugeStubenPreconBuilder(;blocksize=1, kwargs...)
+    SmoothedAggregationPreconBuilder(blocksize,kwargs)
+end
+
+(::RugeStubenPreconBuilder)(A,p)= error("import AlgebraicMultigrid in order to use RugeStubenAMGBuilder")
diff --git a/test/Project.toml b/test/Project.toml
index 76e3ae8..e6a6a68 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -15,6 +15,7 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
+Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
@@ -25,3 +26,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 ExtendableGrids = "1.9"
 IterativeSolvers = "0.9"
+LinearSolve = "2.36.0"
diff --git a/test/test_block.jl b/test/test_block.jl
index 1a055a1..22cabc1 100644
--- a/test/test_block.jl
+++ b/test/test_block.jl
@@ -5,7 +5,7 @@ using ILUZero, AlgebraicMultigrid
 using IterativeSolvers
 using LinearAlgebra
 using Sparspak
-
+using AMGCLWrap
 
 ExtendableSparse.allow_views(::typeof(ilu0))=true
 
@@ -48,4 +48,5 @@ function main(;n=100)
 end
 
 main(n=100)
+
 end
diff --git a/test/test_linearsolve.jl b/test/test_linearsolve.jl
index db161cd..57355be 100644
--- a/test/test_linearsolve.jl
+++ b/test/test_linearsolve.jl
@@ -8,6 +8,10 @@ using LinearSolve
 using ForwardDiff
 using MultiFloats
 
+using AMGCLWrap
+using ILUZero, IncompleteLU, AlgebraicMultigrid
+import Pardiso
+
 f64(x::ForwardDiff.Dual{T}) where {T} = Float64(ForwardDiff.value(x))
 f64(x::Number) = Float64(x)
 const Dual64 = ForwardDiff.Dual{Float64, Float64, 1}
@@ -17,59 +21,105 @@ function test_ls1(T, k, l, m; linsolver = SparspakFactorization())
     b = A*ones(k * l * m)
     x0 = A \ b
     p = LinearProblem(T.(A), T.(b))
-    x1 = solve(p, linsolver)
+    x1 = solve(p, linsolver, abstol=1.0e-12)
     x0 ≈ x1
 end
 
-@test test_ls1(Float64, 10, 10, 10, linsolver = KLUFactorization())
-@test test_ls1(Float64, 25, 40, 1, linsolver = KLUFactorization())
-@test test_ls1(Float64, 100, 1, 1, linsolver = KLUFactorization())
-
-for T in [Float32, Float64, Float64x1, Float64x2, Dual64]
-    println("$T:")
-    @test test_ls1(T, 10, 10, 10, linsolver = SparspakFactorization())
-    @test test_ls1(T, 25, 40, 1, linsolver = SparspakFactorization())
-    @test test_ls1(T, 100, 1, 1, linsolver = SparspakFactorization())
-end
-
 function test_ls2(T, k, l, m; linsolver = SparspakFactorization())
     A = fdrand(T, k, l, m; rand = () -> 1, matrixtype = ExtendableSparseMatrix)
     b = T.(rand(k * l * m))
     p = LinearProblem(A, b)
     x0 = solve(p, linsolver)
     cache = x0.cache
-    x0 = A \ b
-    for i = 4:(k * l * m - 3)
-        A[i, i + 3] -= 1.0e-4
-        A[i - 3, i] -= 1.0e-4
+    x0=copy(x0)
+    nonzeros(A).-=1.0e-4
+    for i = 1:k*l*m
+        A[i, i] += 1.0e-4
     end
 
-    LinearSolve.set_A(cache, A)
-    x1 = solve(p, linsolver)
-    x1 = A \ b
+    reinit!(cache; A, reuse_precs=true)
+    x1 = solve!(cache, linsolver)
     all(x0 .< x1)
 end
 
-@test test_ls1(Float64, 10, 10, 10, linsolver = KLUFactorization())
-@test test_ls1(Float64, 25, 40, 1, linsolver = KLUFactorization())
-@test test_ls1(Float64, 100, 1, 1, linsolver = KLUFactorization())
 
-for T in [Float32, Float64, Float64x1, Float64x2, Dual64]
-    println("$T:")
-    @test test_ls1(T, 10, 10, 10, linsolver = SparspakFactorization())
-    @test test_ls1(T, 25, 40, 1, linsolver = SparspakFactorization())
-    @test test_ls1(T, 100, 1, 1, linsolver = SparspakFactorization())
+@testset "Sparspak" begin
+    for T in [Float32, Float64, Float64x1, Float64x2, Dual64]
+        @test test_ls1(T, 10, 10, 10, linsolver = SparspakFactorization())
+        @test test_ls1(T, 25, 40, 1, linsolver = SparspakFactorization())
+        @test test_ls1(T, 100, 1, 1, linsolver = SparspakFactorization())
+        
+        @test test_ls2(T, 10, 10, 10, linsolver = SparspakFactorization())
+        @test test_ls2(T, 25, 40, 1, linsolver = SparspakFactorization())
+        @test test_ls2(T, 100, 1, 1, linsolver = SparspakFactorization())
+    end
+
+end
+
+factorizations=[UMFPACKFactorization(),
+                SparspakFactorization(),
+                KLUFactorization(reuse_symbolic=false)]
+
+if !Sys.isapple()
+    push!(factorizations,MKLPardisoFactorize())
 end
 
-@test test_ls2(Float64, 10, 10, 10, linsolver = KLUFactorization())
-@test test_ls2(Float64, 25, 40, 1, linsolver = KLUFactorization())
-@test test_ls2(Float64, 100, 1, 1, linsolver = KLUFactorization())
+@testset "Factorizations" begin
+    
+    for factorization in factorizations
 
-for T in [Float32, Float64, Float64x1, Float64x2, Dual64]
-    println("$T:")
-    @test test_ls2(T, 10, 10, 10, linsolver = SparspakFactorization())
-    @test test_ls2(T, 25, 40, 1, linsolver = SparspakFactorization())
-    @test test_ls2(T, 100, 1, 1, linsolver = SparspakFactorization())
+        @test test_ls1(Float64, 10, 10, 10, linsolver = factorization)
+        @test test_ls1(Float64, 25, 40, 1, linsolver = factorization)
+        @test test_ls1(Float64, 100, 1, 1, linsolver = factorization)
+        
+        @test test_ls2(Float64, 10, 10, 10, linsolver = factorization)
+        @test test_ls2(Float64, 25, 40, 1, linsolver = factorization)
+        @test test_ls2(Float64, 100, 1, 1, linsolver = factorization)
+    end
+end
+
+allprecs=[
+    AMGCLWrap.AMGPreconBuilder(),
+    AMGCLWrap.AMGPreconBuilder(),                   
+    AMGCLWrap.RLXPreconBuilder(),                   
+    ExtendableSparse.ILUZeroPreconBuilder(),              
+    ExtendableSparse.ILUZeroPreconBuilder(;blocksize=2),              
+    ExtendableSparse.ILUTPreconBuilder(),              
+    ExtendableSparse.JacobiPreconBuilder(),              
+    ExtendableSparse.SmoothedAggregationPreconBuilder(),
+    ExtendableSparse.RugeStubenPreconBuilder()
+]         
+
+@testset "iterations" begin
+    for precs in allprecs
+        iteration=KrylovJL_GMRES(precs;)
+        
+        @test test_ls1(Float64, 10, 10, 10, linsolver = iteration)
+        @test test_ls1(Float64, 25, 40, 1, linsolver = iteration)
+        @test test_ls1(Float64, 100, 1, 1, linsolver = iteration)
+        
+        @test test_ls2(Float64, 10, 10, 10, linsolver = iteration)
+        @test test_ls2(Float64, 25, 40, 1, linsolver = iteration)
+        @test test_ls2(Float64, 100, 1, 1, linsolver = iteration)
+    end
+end
+
+
+luprecs=[ExtendableSparse.LinearSolvePreconBuilder(factorization) for  factorization in factorizations]
+
+@testset "block preconditioning" begin
+    n=100
+    A=fdrand(n,n)
+    partitioning=A->[1:2:size(A,1), 2:2:size(A,1)]
+    sol0=ones(n^2)
+    b=A*ones(n^2);
+    
+    for precs in vcat(allprecs, luprecs)
+        iteration=KrylovJL_CG(precs=BlockPreconBuilder(;precs, partitioning))
+        p=LinearProblem(A,b)
+        sol=solve(p, KrylovJL_CG(;precs), abstol=1.0e-12)
+        @test isapprox(sol, sol0, atol=1e-6)
+    end
 end
 
 end
diff --git a/test/test_pardiso.jl b/test/test_pardiso.jl
index 24f6081..1bc07b7 100644
--- a/test/test_pardiso.jl
+++ b/test/test_pardiso.jl
@@ -14,7 +14,4 @@ include("test_lu.jl")
 @test test_lu2(Float64, 10, 10, 10, lufac = PardisoLU())
 @test test_lu2(Float64, 25, 40, 1, lufac = PardisoLU())
 @test test_lu2(Float64, 1000, 1, 1, lufac = PardisoLU())
-# @test test_lu2(10,10,10,lufac=PardisoLU(mtype=2))
-# @test test_lu2(25,40,1,lufac=PardisoLU(mtype=2))
-# @test test_lu2(1000,1,1,lufac=PardisoLU(mtype=2))
 end