Add support for LinearSolve (#24)

* Add set_A and LinearProblem methods if LinearSolve is loaded.
* format all files
* format & update docs

Author: j-fu

.JuliaFormatter.toml

@@ -0,0 +1,4 @@
+style = "sciml"
+always_for_in = false
+separate_kwargs_with_semicolon = true
+format_markdown = true

.gitignore

@@ -1,7 +1,7 @@
 QUARRY
-.project.jl
 docs/build
 *~
 Manifest.toml
-.repl_history.jl
+.repl_history
+quarry
 

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <[email protected]>"]
-version = "0.9.5"
+version = "0.9.6"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -16,9 +16,9 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DocStringExtensions = "0.8, 0.9"
-ILUZero= "0.2"
+ILUZero = "0.2"
 Requires = "1.1.3"
-Sparspak = "0.3.3"
+Sparspak = "0.3.6"
 julia = "1.6"
 
 [extras]

README.md

@@ -7,10 +7,48 @@
 
 Sparse matrix class with efficient successive insertion of entries and entry update, supporting general number types.
 
+## Summary
+
+The package allows for efficient assembly of a sparse matrix without
+the need to handle intermediate arrays:
+
+```
+using ExtendableSparse
+A=ExtendableSparseMatrix(10,10)
+A[1,1]=1
+for i = 1:9
+   A[i + 1, i] += -1
+   A[i, i + 1] += -1
+   A[i, i] += 1
+   A[i + 1, i + 1] += 1
+end
+b=ones(10)
+x=A\b
+```
+
+While one could replace here  `ExtendableSparseMatrix(10,10)` by `spzeros(10,10)`, the later is inefficient for large problems. Without this package, the general advise is  to [construct a sparse matrix via the COO format](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#Sparse-Vector-and-Matrix-Constructors).
+
+Instead of `\`, the methods from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) can be used:
+
+```
+using LinearSolve
+p=LinearProblem(A,b)
+LinearSolve.solve(p)
+```
+
+With the help of [Sparspak.jl](https://github.com/PetrKryslUCSD/Sparspak.jl), these examples work for general number types.
+
+`sparse(A)` creates a standard `SparseMatrixCSC` from a filled `ExtendableSparseMatrix` which can be used e.g. to create preconditioners. So one can instead perform e.g.
+
+```
+LinearSolve.solve(p, KrylovJL_CG(); Pl = ILUZero.ilu0(sparse(A)))
+```
+
 ## Rationale
+
 Without an intermediate data structure, efficient successive insertion/update of possibly duplicate entries in random order into a standard compressed column storage structure appears to be not possible. The package introduces `ExtendableSparseMatrix`, a delegating wrapper containing a Julia standard `SparseMatrixCSC` struct for performing linear algebra operations and a `SparseMatrixLNK` struct realising a linked list based (but realised in vectors) format collecting new entries.
 
-The later is modeled after the linked list sparse matrix format described in the [whitepaper](https://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps) by Y. Saad. See also exercise P.3-16  in his [book](https://www-users.cs.umn.edu/~saad/IterMethBook_2ndEd.pdf).
+The later is modeled after the linked list sparse matrix format described in the [whitepaper](https://www-users.cs.umn.edu/%7Esaad/software/SPARSKIT/paper.ps) by Y. Saad. See also exercise P.3-16  in his [book](https://www-users.cs.umn.edu/%7Esaad/IterMethBook_2ndEd.pdf).
 
 Any linear algebra method on `ExtendableSparseMatrix` starts with a `flush!` method which adds the LNK entries and the existing CSC entries into a new CSC struct and resets the LNK struct.
 
@@ -25,10 +63,10 @@ number types  not supported by Julia's standard implementation. Notably, this  i
 This package assumes that a  $m \times n$  matrix is sparse if *each* row and *each* column have less than $C$ entries with
 $C << n$ and $C << m$ . Adding a full matrix row will be a performance hit.
 
-
 ## Working with ForwardDiff
 
 In particular, it cooperates with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) when it comes to the assembly of a sparse jacobian. For a function 'f!(y,x)' returning it's result in a vector `y`, one can use e.g.
+
 ````
 x=...
 y=zeros(n)
@@ -38,46 +76,46 @@ jac=DiffResults.jacobian(dresult)
 h=jac\x
 ````
 
-However, without a priori information on sparsity, ForwardDiff calls element insertion for the full range of n^2 indices, 
+However, without a priori information on sparsity, ForwardDiff calls element insertion for the full range of n^2 indices,
 leading to a O(n^2) scaling behavior due to the nevertheless necessary search operations, see  this [discourse thread](https://discourse.julialang.org/t/non-sorted-sparsematrixcsc/37133).
 
 ## updateindex!
-In addition, the package provides a method `updateindex!(A,op,v,i,j)` for both `SparseMatrixCSC` and for `ExtendableSparse` which allows to update a matrix element with one index search instead of two. It allows to replace e.g. `A[i,j]+=v` by `updateindex!(A,+,v,i,j)`. The former operation is lowered to 
+
+In addition, the package provides a method `updateindex!(A,op,v,i,j)` for both `SparseMatrixCSC` and for `ExtendableSparse` which allows to update a matrix element with one index search instead of two. It allows to replace e.g. `A[i,j]+=v` by `updateindex!(A,+,v,i,j)`. The former operation is lowered to
+
 ````
 %1 = Base.getindex(A, 1, 2)
 %2 = %1 + 3
 Base.setindex!(A, %2, 1, 2)
 ````
+
 triggering two index searches, one for `getindex!` and another one for `setindex!`.
 
 See [Julia issue #15630](https://github.com/JuliaLang/julia/issues/15630) for a discussion on this.
 
-
-
 ## Factorizations and Preconditioners
+
 The package provides a common API for factorizations and preconditioners supporting
 series of solutions of similar problem as they occur during nonlinear and transient solves.
 For details, see the [corresponding documentation](https://j-fu.github.io/ExtendableSparse.jl/stable/iter/).
 
+With the advent of LinearSolve.jl, this functionality probably will be reduced to some core cases.
 
 ### Interfaces to other packages
 
-Compatibility with [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) is under development.
-
-The package provides interfaces to other sparse matrix solvers and preconditioners. Dependencies on these
+The package directly provides interfaces to other sparse matrix solvers and preconditioners. Dependencies on these
 packages are handeled via [Requires.jl](https://github.com/JuliaPackaging/Requires.jl).
 Currently, support includes:
 
-- [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl) (both ["project Pardiso"](https://pardiso-project.org)
-  and [MKL Pardiso](https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-interface.html))
-- [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl)
-- [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl) (Ruge-Stüben AMG)
+  - [Pardiso.jl](https://github.com/JuliaSparse/Pardiso.jl) (both ["project Pardiso"](https://pardiso-project.org)
+    and [MKL Pardiso](https://software.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-interface.html))
+  - [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl)
+  - [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl) (Ruge-Stüben AMG)
 
 For a similar approach, see [Preconditioners.jl](https://github.com/mohamed82008/Preconditioners.jl)
 
 ## Alternatives
 
 You may also evaluate alternatives to this package:
 
-- [DynamicSparseArrays.jl](https://github.com/atoptima/DynamicSparseArrays.jl)
-  
+  - [DynamicSparseArrays.jl](https://github.com/atoptima/DynamicSparseArrays.jl)

docs/Project.toml

@@ -5,8 +5,10 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
+ILUZero = "88f59080-6952-5380-9ea5-54057fb9a43f"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
 

docs/make.jl

@@ -1,24 +1,24 @@
-push!(LOAD_PATH,"../src/")
-using Documenter, ExtendableSparse,Pardiso,AlgebraicMultigrid,IncompleteLU,Sparspak
+push!(LOAD_PATH, "../src/")
+using Documenter, ExtendableSparse, Pardiso, AlgebraicMultigrid, IncompleteLU, Sparspak
 
 function mkdocs()
-    makedocs(sitename="ExtendableSparse.jl",
+    makedocs(; sitename = "ExtendableSparse.jl",
              modules = [ExtendableSparse],
              doctest = true,
              clean = false,
              authors = "J. Fuhrmann",
-             repo="https://github.com/j-fu/ExtendableSparse.jl",
-             pages=[
-                 "Home"=>"index.md",
+             repo = "https://github.com/j-fu/ExtendableSparse.jl",
+             pages = [
+                 "Home" => "index.md",
                  "example.md",
+                 "linearsolve.md",
                  "extsparse.md",
                  "iter.md",
-#                 "linearsolve.md",
                  "internal.md",
                  "changes.md",
              ])
 end
 
 mkdocs()
 
-deploydocs(repo = "github.com/j-fu/ExtendableSparse.jl.git")
+deploydocs(; repo = "github.com/j-fu/ExtendableSparse.jl.git")

docs/src/changes.md

@@ -1,67 +1,101 @@
 # Changes
+
+## v0.9.6, Jan 22, 2023
+
+  - support for LinearSolve.jl
+
+## v0.9.1 ... v0.9.5, Jan 5, 2023
+
+  - update interface with Sparspak
+  - Addition, multiplication etc. for ExtendableSparse
+  - ILUZero.jl as dependency for ILUZeroPreconditioner
+  - fix ILU0Preconditioner backsubstitution for unsymmetric matrices
+  - Add LinearAlgebra.lu etc. for ExtendableSparse
+
 ## v0.9, Sept 28, 2022
-- `\` support for general number types and non-GPL system images based on Sparspak.jl 0.3.0
+
+  - `\` support for general number types and non-GPL system images based on Sparspak.jl 0.3.0
+
 ## v0.8, Sept 1, 2022
-- Remove LinearSolve compatibility in favor of (future) interfacing via AbstractSparseMatrixCSC
-- Add Sparspak LU factorization
-- Add handling of GPL free sysimage build
+
+  - Remove LinearSolve compatibility in favor of (future) interfacing via AbstractSparseMatrixCSC
+  - Add Sparspak LU factorization
+  - Add handling of GPL free sysimage build
+
 ## v0.7, August 19, 2022
-- Require Julia 1.6
-- first steps to compatibility with LinearSolve.jl
+
+  - Require Julia 1.6
+  - first steps to compatibility with LinearSolve.jl
+
 ## v0.6, April 20, 2021
-- use type parameters to describe factorizations
+
+  - use type parameters to describe factorizations
+
 ## v0.5, April 10, 2021
-- Introduce lu/lu! , factorize/factorize!, unifying LU factorizations and preconditioners
-- Interface packages: Pardiso, AlgebraicMultigrid, IncompleteLU via Requires.jl
+
+  - Introduce lu/lu! , factorize/factorize!, unifying LU factorizations and preconditioners
+  - Interface packages: Pardiso, AlgebraicMultigrid, IncompleteLU via Requires.jl
+
 ## v0.4, March 2021
-- Fix handling of Symmetrix matrices
-- `rawupdateindex` does not check for entering zeros
-- Compare with `COO` method
-- Benchmarks in documentation
+
+  - Fix handling of Symmetrix matrices
+  - `rawupdateindex` does not check for entering zeros
+  - Compare with `COO` method
+  - Benchmarks in documentation
+
 ## v0.3.7, March 20, 2021
-- Added parallel jacobi preconditioner (thanks, @jkr)
-- Fixes ldiv
-- Added simple iterative solver
-- Documentation update
-- Tests for precondioners, fdrand
+
+  - Added parallel jacobi preconditioner (thanks, @jkr)
+  - Fixes ldiv
+  - Added simple iterative solver
+  - Documentation update
+  - Tests for precondioners, fdrand
 
 ## v0.3.0, April 10, 2020
-- Don't create new entry if the value to be assigned is zero, making things consistent with SparseMatrixCSC and ForwardDiff 
-  as suggested by @MaximilianJHuber
+
+  - Don't create new entry if the value to be assigned is zero, making things consistent with SparseMatrixCSC and ForwardDiff
+    as suggested by @MaximilianJHuber
 
 ## v0.2.5, Jan 26, 2020
-- fixed allocations in  Base.+
-- added updateindex! method 
-- provide fdrand and fdrand! matrix constructors
-- automatic benchmarks in examples
+
+  - fixed allocations in  Base.+
+  - added updateindex! method
+  - provide fdrand and fdrand! matrix constructors
+  - automatic benchmarks in examples
 
 ## v0.2.4, Jan 19, 2020
-- Allow preconditioner creation directly from CSC Matrix
-- Rename AbstractPreconditioner to AbstractExtendablePreconditioner
+
+  - Allow preconditioner creation directly from CSC Matrix
+  - Rename AbstractPreconditioner to AbstractExtendablePreconditioner
 
 ## v0.2.3, Jan 15, 2020
-- Started to introduce preconditioners (undocumented)
+
+  - Started to introduce preconditioners (undocumented)
 
 ## v0.2.3, Jan 8, 2020
-- added norm, cond, opnorm methods
-- resize! instead of push! when adding entries should trigger less allocation operations
+
+  - added norm, cond, opnorm methods
+  - resize! instead of push! when adding entries should trigger less allocation operations
 
 ## v0.2.2. Dec 23, 2019
-- What used to be `_splice`  is now `+` and allows now real addition (resulting in a CSC matrix)
-- Added constructors of LNK matrix from CSC matrix and vice versa
-- reorganized tests
+
+  - What used to be `_splice`  is now `+` and allows now real addition (resulting in a CSC matrix)
+  - Added constructors of LNK matrix from CSC matrix and vice versa
+  - reorganized tests
 
 ## v0.2.1 Dec 22, 2019
-- Tried to track down the source from which I learned the linked list based struct in order
-  to document this. Ended up with SPARSEKIT of Y.Saad, however I believe this 
-  already was in SPARSEPAK by Chu,George,Liu.
-- Internal rename of SparseMatrixExtension to SparseMatrixLNK. 
+
+  - Tried to track down the source from which I learned the linked list based struct in order
+    to document this. Ended up with SPARSEKIT of Y.Saad, however I believe this
+    already was in SPARSEPAK by Chu,George,Liu.
+  - Internal rename of SparseMatrixExtension to SparseMatrixLNK.
 
 ## v0.2 Dec 21, 2019
-- more interface methods delegating to csc, in particular mul! and ldiv!
-- lazy creation of extendable part: don't create idle memory
-- nicer constructors
-  
+
+  - more interface methods delegating to csc, in particular mul! and ldiv!
+  - lazy creation of extendable part: don't create idle memory
+  - nicer constructors
+
 ## V0.1, July 2019
-- Initial release
 
+  - Initial release