Merge pull request #88 from gridap/develop

NullSpaces

NEWS.md

@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added Vanka-like smoothers in serial. Since PR[#68](https://github.com/gridap/GridapSolvers.jl/pull/68).
 - Added `StaggeredFEOperators` and `StaggeredFESolvers`. Since PR[#84](https://github.com/gridap/GridapSolvers.jl/pull/84).
 - Added `RichardsonLinearSolver`. Since PR[#87](https://github.com/gridap/GridapSolvers.jl/pull/87).
+- Added `NullspaceSolver` for serial. Since PR[#88](https://github.com/gridap/GridapSolvers.jl/pull/88).
 
 ## Previous versions
 

docs/src/LinearSolvers.md

@@ -54,6 +54,12 @@ Given a linear system ``Ax = b``, a **preconditioner** is an operator that takes
   GMGLinearSolver
 ```
 
+## Nullspaces
+
+```@docs
+  NullspaceSolver
+```
+
 ## Wrappers
 
 ### PETSc

src/LinearSolvers/LinearSolvers.jl

@@ -29,6 +29,7 @@ export SchwarzLinearSolver
 export RichardsonLinearSolver
 
 export CallbackSolver
+export NullspaceSolver
 
 # Wrappers for IterativeSolvers.jl
 export IS_ConjugateGradientSolver
@@ -68,5 +69,6 @@ include("MatrixSolvers.jl")
 include("SchwarzLinearSolvers.jl")
 
 include("CallbackSolver.jl")
+include("NullspaceSolvers.jl")
 
 end

src/LinearSolvers/NullspaceSolvers.jl

@@ -0,0 +1,120 @@
+
+"""
+    struct NullspaceSolver{A,B}
+      solver :: A <: LinearSolver
+      kernel :: B <: Union{AbstractMatrix, Vector{<:AbstractVector}}
+      constrain_matrix :: Bool = true
+    end
+
+Solver that computes the solution of a linear system `A⋅x = b` with a kernel constraint,
+i.e the returned solution is orthogonal to the provided kernel.
+
+We assume the kernel is provided as a matrix `K` of dimensions `(k,n)` with `n` the dimension
+of the original system and `k` the number of kernel vectors.
+
+Two modes of operation are supported:
+
+- If `constrain_matrix` is `true`, the solver will explicitly constrain the system matrix. 
+  I.e we consider the augmented system `Â⋅x̂ = b̂` where
+  -  `Ak = [A, K'; K, 0]`
+  -  `x̂ = [x; λ]`
+  -  `b̂ = [b; 0]`
+  This is often the only option for direct solvers, which require the system matrix to be
+  invertible. This is only supported in serial (its performance bottleneck in parallel).
+
+- If `constrain_matrix` is `false`, the solver preserve the original system and simply 
+  project the initial guess and the solution onto the orthogonal complement of the kernel. 
+  This option is more suitable for iterative solvers, which usually do not require the 
+  system matrix to be invertible (e.g. GMRES, BiCGStab, etc).
+"""
+struct NullspaceSolver{A,B}
+  solver    :: A
+  nullspace :: B
+  constrain_matrix :: Bool
+
+  function NullspaceSolver(
+    solver::LinearSolver,
+    nullspace::NullSpace; 
+    constrain_matrix::Bool = true
+  )
+    A, B = typeof(solver), typeof(nullspace)
+    new{A,B}(solver, nullspace, constrain_matrix)
+  end
+end
+
+struct NullspaceSolverSS{A} <: Algebra.SymbolicSetup
+  solver :: A
+end
+
+function Gridap.Algebra.symbolic_setup(solver::NullspaceSolver,A)
+  return NullspaceSolverSS(solver)
+end
+
+struct NullspaceSolverNS{S}
+  solver
+  ns
+  caches
+end
+
+function Algebra.numerical_setup(ss::NullspaceSolverSS, A::AbstractMatrix)
+  solver = ss.solver
+  N = solver.nullspace
+  nK, nV = size(N)
+  @assert size(A,1) == size(A,2) == nV
+  if solver.constrain_matrix
+    K = SolverInterfaces.matrix_representation(N)
+    mat = [A K; K' zeros(nK,nK)] # TODO: Explore reusing storage for A
+  else
+    SolverInterfaces.make_orthonormal!(N)
+    mat = A
+  end
+  S = ifelse(solver.constrain_matrix, :constrained, :projected)
+  ns = numerical_setup(symbolic_setup(solver.solver, mat), mat)
+  caches = (allocate_in_domain(mat), allocate_in_domain(mat))
+  return NullspaceSolverNS{S}(solver, ns, caches)
+end
+
+function Algebra.numerical_setup!(ns::NullspaceSolverNS, A::AbstractMatrix)
+  solver = ns.solver
+  N = solver.nullspace
+  nK, nV = size(N)
+  @assert size(A,1) == size(A,2) == nV
+  if solver.constrain_matrix
+    K = SolverInterfaces.matrix_representation(N)
+    mat = [A K; K' zeros(nK,nK)] # TODO: Explore reusing storage for A
+  else
+    mat = A
+  end
+  numerical_setup!(ns.ns, mat)
+  return ns
+end
+
+function Algebra.solve!(x::AbstractVector, ns::NullspaceSolverNS{:constrained}, b::AbstractVector)
+  solver = ns.solver
+  A_ns, caches = ns.ns, ns.caches
+  N = solver.nullspace
+  w1, w2 = caches
+  nK, nV = size(N)
+  @assert length(x) == nV
+  @assert length(b) == nV
+  w1[1:nV] .= x
+  w1[nV+1:nV+nK] .= 0.0
+  w2[1:nV] .= b
+  w2[nV+1:nV+nK] .= 0.0
+  solve!(w1, A_ns, w2)
+  x .= w1[1:nV]
+  return x
+end
+
+function Algebra.solve!(x::AbstractVector, ns::NullspaceSolverNS{:projected}, b::AbstractVector)
+  solver = ns.solver
+  A_ns, caches = ns.ns, ns.caches
+  N = solver.nullspace
+  w1, w2 = caches
+
+  w1, α = SolverInterfaces.project!(w1, N, x)
+  x .-= w1
+  solve!(x, A_ns, b)
+
+  return x
+end

src/SolverInterfaces/NullSpaces.jl

@@ -0,0 +1,139 @@
+struct NullSpace{T,VT}
+  V :: VT
+  function NullSpace(
+    V::AbstractVector{<:AbstractVector{T}}
+  ) where T
+    n = length(first(V))
+    @assert all(length(v) == n for v in V)
+    VT = typeof(V)
+    new{T,VT}(V)
+  end
+end
+
+Base.size(N::NullSpace) = (length(N.V),length(first(N.V)))
+Base.size(N::NullSpace, i::Int) = size(N)[i]
+Base.merge(a::NullSpace{T}, b::NullSpace{T}) where T = NullSpace(vcat(a.V, b.V))
+
+function matrix_representation(N::NullSpace)
+  return stack(N.V)
+end
+
+NullSpace(v::AbstractVector{<:Number}) = NullSpace([v])
+
+function NullSpace(A::Matrix)
+  V = eachcol(nullspace(A))
+  return NullSpace(V)
+end
+
+function NullSpace(f::Function,X::FESpace,Λ::FESpace)
+  A = assemble_matrix(f, X, Λ)
+  return NullSpace(A)
+end
+
+function is_orthonormal(N::NullSpace; tol = 1.e-12)
+  for w in N.V
+    !(abs(norm(w) - 1.0) < tol) && return false
+  end
+  return is_orthogonal(N, tol = tol)
+end
+
+function is_orthogonal(N::NullSpace; tol = 1.e-12)
+  for (k,w) in enumerate(N.V)
+    for v in N.V[k+1:end]
+      !(abs(dot(w, v)) < tol) && return false
+    end
+  end
+  return true
+end
+
+function is_orthogonal(N::NullSpace, v::AbstractVector; tol = 1.e-12)
+  @assert length(v) == size(N,2)
+  for w in N.V
+    !(abs(dot(w, v)) < tol) && return false
+  end
+  return true
+end
+
+function is_orthogonal(N::NullSpace, A::AbstractMatrix; tol = 1.e-12)
+  @assert length(v) == size(N,2)
+  v = allocate_in_range(A)
+  for w in N.V
+    mul!(v, A, w)
+    !(abs(norm(v)) < tol) && return false
+  end
+  return true
+end
+
+function make_orthonormal!(N::NullSpace; method = :gram_schmidt)
+  if method == :gram_schmidt
+    gram_schmidt!(N.V)
+  elseif method == :modified_gram_schmidt
+    modified_gram_schmidt!(N.V)
+  else
+    error("Unknown method: $method")
+  end
+  return N
+end
+
+function gram_schmidt!(V::AbstractVector{<:AbstractVector{T}}) where T
+  n = length(V)
+  for j in 1:n
+    for i in 1:j-1
+      α = dot(V[j], V[i])
+      V[j] .-= α .* V[i]
+    end
+    V[j] ./= norm(V[j])
+  end
+  return V
+end
+
+function modified_gram_schmidt!(V::AbstractVector{<:AbstractVector{T}}) where T
+  n = length(V)
+  for j in 1:n
+    V[j] ./= norm(V[j])
+    for i in j+1:n
+      α = dot(V[j], V[i])
+      V[i] .-= α .* V[j]
+    end
+  end
+  return V
+end
+
+function project(N::NullSpace,v)
+  p = similar(v)
+  project!(p,N,v)
+end
+
+function project!(p,N::NullSpace{T},v) where T
+  @assert length(v) == size(N,2)
+  α = zeros(T, size(N,1))
+  fill!(p,0.0)
+  for (k,w) in enumerate(N.V)
+    α[k] = dot(v, w)
+    p .+= α[k] .* w
+  end
+  return p, α
+end
+
+function make_orthogonal!(N::NullSpace{T},v) where T
+  @assert length(v) == size(N,2)
+  α = zeros(T, size(N,1))
+  for (k,w) in enumerate(N.V)
+    α[k] = dot(v, w)
+    v .-= α[k] .* w
+  end
+  return v, α
+end
+
+function reconstruct(N::NullSpace,v,α)
+  w = copy(v)
+  reconstruct!(N,w,α)
+  return w
+end
+
+function reconstruct!(N::NullSpace,v,α)
+  for (k,w) in enumerate(N.V)
+    v .+= α[k] .* w
+  end
+  return v
+end

src/SolverInterfaces/SolverInterfaces.jl

@@ -1,5 +1,7 @@
 module SolverInterfaces
 
+using LinearAlgebra
+
 using Gridap
 using Gridap.Helpers
 using Gridap.Algebra
@@ -15,11 +17,14 @@ include("GridapExtras.jl")
 include("SolverTolerances.jl")
 include("ConvergenceLogs.jl")
 include("SolverInfos.jl")
+include("NullSpaces.jl")
 
 export SolverVerboseLevel, SolverConvergenceFlag
 export SolverTolerances, get_solver_tolerances, set_solver_tolerances!
 export ConvergenceLog, init!, update!, finalize!, reset!, print_message, set_depth!
 
 export SolverInfo
 
+export NullSpace
+
 end