add LDLFactorizations tutorial (#32)

tutorials/introduction-to-ldlfactorizations/Project.toml

@@ -0,0 +1,24 @@
+[deps]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixMarket = "4d4711f2-db25-561a-b6b3-d35e7d4047d3"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QPSReader = "10f199a5-22af-520b-b891-7ce84a7b1bd0"
+QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
+RipQP = "1e40b3f8-35eb-4cd8-8edd-3e515bb9de08"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SparseMatricesCOO = "fa32481b-f100-4b48-8dc8-c62f61b13870"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+
+[compat]
+JSON = "0.21"
+LDLFactorizations = "0.8"
+MatrixMarket = "0.3"
+Plots = "1"
+QPSReader = "0.2"
+QuadraticModels = "0.7"
+RipQP = "0.3"
+SparseMatricesCOO = "0.1"
+TimerOutputs = "0.5"

tutorials/introduction-to-ldlfactorizations/index.jmd

@@ -0,0 +1,77 @@
+---
+title: "LDLFactorizations tutorial"
+tags: ["linear", "factorization", "ldlt"]
+author: "Geoffroy Leconte"
+---
+
+LDLFactorizations.jl is a translation of Tim Davis's Concise LDLᵀ Factorization, part of [SuiteSparse](http://faculty.cse.tamu.edu/davis/suitesparse.html) with several improvements.
+
+This package is appropriate for matrices A that possess a factorization of the
+form LDLᵀ without pivoting, where L is unit lower triangular and D is *diagonal* (indefinite in general), including definite and quasi-definite matrices.
+
+## A basic example
+
+```julia
+using LinearAlgebra
+n = 10
+A0 = rand(n, n)
+A = A0 * A0' + I # A is symmetric positive definite
+b = rand(n)
+```
+
+We solve the system $A x = b$ using LDLFactorizations.jl:
+
+```julia 
+using LDLFactorizations
+Au = Symmetric(triu(A), :U) # get upper triangle and apply Symmetric wrapper
+LDL = ldl(Au)
+x = LDL \ b
+```
+
+## A more performance-focused example
+
+We build a problem with sparse arrays.
+
+```julia
+using SparseArrays
+n = 100
+# create create a SQD matrix A:
+A0 = sprand(Float64, n, n, 0.1)
+A1 = A0 * A0' + I
+A = [A1   A0;
+     A0' -A1]
+b = rand(2 * n)
+```
+
+Now if we want to use the factorization to solve multiple systems that have 
+the same sparsity pattern as A, we only have to use `ldl_analyze` once.
+
+```julia 
+Au = Symmetric(triu(A), :U) # get upper triangle and apply Symmetric wrapper
+x = similar(b)
+
+LDL = ldl_analyze(Au) # symbolic analysis
+ldl_factorize!(Au, LDL) # factorization
+ldiv!(x, LDL, b) # solve in-place (we could use ldiv!(LDL, b) if we want to overwrite b)
+
+Au.data.nzval .+= 1.0 # modify Au without changing the sparsity pattern
+ldl_factorize!(Au, LDL) 
+ldiv!(x, LDL, b)
+```
+
+## Dynamic Regularization
+
+When the matrix to factorize is (nearly) singular and the factorization encounters (nearly) zero pivots, 
+if we know the signs of the pivots and if they are clustered by signs (for example, the 
+`n_d` first pivots are positive and the other pivots are negative before permuting), we can use:
+
+```julia
+ϵ = sqrt(eps())
+Au = Symmetric(triu(A), :U)
+LDL = ldl_analyze(Au)
+LDL.tol = ϵ
+LDL.n_d = 10
+LDL.r1 = 2 * ϵ # if any of the n_d first pivots |D[i]| < ϵ, then D[i] = sign(LDL.r1) * max(abs(D[i] + LDL.r1), abs(LDL.r1))
+LDL.r2 = -ϵ # if any of the n - n_d last pivots |D[i]| < ϵ, then D[i] = sign(LDL.r2) * max(abs(D[i] + LDL.r2), abs(LDL.r2))
+ldl_factorize!(Au, LDL)
+```