Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

NullSpaces #88

Draft
wants to merge 2 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/CI.yml
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ jobs:
cd $SOURCES_DIR
./configure --prefix=$PETSC_INSTALL --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 \
--download-mumps --download-scalapack --download-parmetis --download-metis \
--download-ptscotch --with-debugging --with-x=0 --with-shared=1 \
--download-fblaslapack --download-ptscotch --with-debugging --with-x=0 --with-shared=1 \
--with-mpi=1 --with-64-bit-indices
make
make install
Expand Down
2 changes: 2 additions & 0 deletions src/LinearSolvers/LinearSolvers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ export SchurComplementSolver
export SchwarzLinearSolver

export CallbackSolver
export NullspaceSolver

# Wrappers for IterativeSolvers.jl
export IS_ConjugateGradientSolver
Expand Down Expand Up @@ -66,5 +67,6 @@ include("MatrixSolvers.jl")
include("SchwarzLinearSolvers.jl")

include("CallbackSolver.jl")
include("NullspaceSolvers.jl")

end
122 changes: 122 additions & 0 deletions src/LinearSolvers/NullspaceSolvers.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,122 @@

"""
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
display(mat)
display(A)
display(K)
else
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)
x .+= w1
return x
end
139 changes: 139 additions & 0 deletions src/SolverInterfaces/NullSpaces.jl
Original file line number Diff line number Diff line change
@@ -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
5 changes: 5 additions & 0 deletions src/SolverInterfaces/SolverInterfaces.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
module SolverInterfaces

using LinearAlgebra

using Gridap
using Gridap.Helpers
using Gridap.Algebra
Expand All @@ -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
66 changes: 66 additions & 0 deletions test/_dev/nullspaces.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@

using Gridap
using Gridap.Algebra

using GridapSolvers
using GridapSolvers.SolverInterfaces, GridapSolvers.LinearSolvers

using GridapSolvers.SolverInterfaces: NullSpace, is_orthonormal, is_orthogonal, gram_schmidt!, modified_gram_schmidt!
using GridapSolvers.SolverInterfaces: project, project!, make_orthogonal!, reconstruct, reconstruct!

A = [1.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
N = NullSpace(A)

is_orthonormal(N)
is_orthogonal(N,[1.0,0.0,0.0])
is_orthogonal(N,[0.0,1.0,0.0])

v = [1.0,2.0,3.0]
p, α = project(N,v)
w, β = make_orthogonal!(N,copy(v))
u = reconstruct(N,w,α)
u ≈ v

V1 = [[2.0,1.0,1.0],[1.0,2.0,1.0],[1.0,1.0,1.0]]
gram_schmidt!(V1)
is_orthonormal(NullSpace(V1))

V2 = [[2.0,1.0,1.0],[1.0,2.0,1.0],[1.0,1.0,1.0]]
modified_gram_schmidt!(V2)
is_orthonormal(NullSpace(V2))
V1 ≈ V2


############################################################################################

model = CartesianDiscreteModel((0,1,0,1),(4,4))
Ω = Triangulation(model)
Γ = BoundaryTriangulation(model)

dΩ = Measure(Ω,4)
dΓ = Measure(Γ,4)
n = get_normal_vector(Γ)

V = FESpace(model, ReferenceFE(lagrangian,Float64,1))

u_exact(x) = x[1]*x[2]
f(x) = Δ(u_exact)(x)
a(u,v) = ∫(∇(u)⋅∇(v))*dΩ
l(v) = ∫(f*v)*dΩ + ∫(v*(∇(u_exact)⋅n))*dΓ

op = AffineFEOperator(a,l,V,V)
A, b = get_matrix(op), get_vector(op)

N = NullSpace(ones(size(A,2)))
K = SolverInterfaces.matrix_representation(N)
norm(A*K)

s = GMRESSolver(10;verbose=true)
s = LUSolver()
solver = NullspaceSolver(s,N;constrain_matrix=true)
ns = numerical_setup(symbolic_setup(solver,A),A)

x = allocate_in_domain(A)
fill!(x,0.0)
solve!(x,ns,b)
norm(A*x-b)
Loading