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.

Sign up for GitHub

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

Sym storage not working #16

Merged
merged 4 commits into from
Feb 6, 2020
Merged

Sym storage not working #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ecc9f5d
non working femdriver for symmetric storage
fverdugo Feb 5, 2020
282900c
Fix test/femdriver
victorsndvg Feb 6, 2020
017cd64
Add feature suggested in https://github.com/gridap/GridapPardiso.jl/i…
victorsndvg Feb 6, 2020
d88a2a0
Implement comments to the previous commit
victorsndvg Feb 6, 2020
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 README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ t_Ω = AffineFETerm(

op = AffineFEOperator(SparseMatrixCSR{1,Float64,Int},V,U,t_Ω)

ls = PardisoSolver()
ls = PardisoSolver(op)
solver = LinearFESolver(ls)

uh = solve(solver,op)
Expand Down
2 changes: 2 additions & 0 deletions src/GridapPardiso.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@ module GridapPardiso
using Libdl
using SparseArrays
using Gridap.Algebra
using Gridap.FESpaces
using Gridap.Helpers

import Gridap.Algebra: LinearSolver
import Gridap.Algebra: SymbolicSetup, NumericalSetup
Expand Down
54 changes: 51 additions & 3 deletions src/LinearSolver.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,23 +83,39 @@ PardisoSolver() = PardisoSolver(MTYPE_REAL_NON_SYMMETRIC,
PardisoSolver constructor overloading with default values.
Returns a PardisoSolver given its matrix type.
"""
PardisoSolver(mtype) =
PardisoSolver(mtype::Int) =
PardisoSolver(mtype, new_iparm(), MSGLVL_QUIET, new_pardiso_handle())

PardisoSolver(::Type{<:AbstractSparseMatrix{Tv,Ti}}) where {Tv,Ti}=
PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC)

PardisoSolver(::Type{SymSparseMatrixCSR{Bi,Tv,Ti}}) where {Bi,Tv,Ti}=
PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE)

PardisoSolver(op::FEOperator) = PardisoSolver(get_matrix_type(op))

function get_matrix_type(op::FEOperator)
@abstractmethod
end

function get_matrix_type(op::AffineFEOperator)
typeof(get_matrix(op))
end

"""
function PardisoSolver(mtype, iparm)
PardisoSolver constructor overloading with default values.
Returns a PardisoSolver given its matrix type and Pardiso parameters.
"""
PardisoSolver(mtype, iparm) =
PardisoSolver(mtype::Int, iparm) =
PardisoSolver(mtype, iparm, MSGLVL_QUIET, new_pardiso_handle())

"""
function PardisoSolver(mtype, iparm, msglvl)
PardisoSolver constructor overloading with default values.
Returns a PardisoSolver given its matrix type, Pardiso parameters and verbosity.
"""
PardisoSolver(mtype, iparm, msglvl) =
PardisoSolver(mtype::Int, iparm, msglvl) =
PardisoSolver(mtype, iparm, msglvl, new_pardiso_handle())


Expand All @@ -111,6 +127,10 @@ Returns a PardisoSymbolicSetup from a given AbstractSparseMatrix.
function build_PardisoSymbolicSetup(phase::Integer,
mat::AbstractSparseMatrix{T,Ti},
solver::PardisoSolver) where {T,Ti}
if !(solver.mtype in (MTYPE_REAL_STRUCTURALLY_SYMMETRIC,
MTYPE_REAL_NON_SYMMETRIC))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
PardisoSymbolicSetup(phase,SparseMatrixCSC{T,Ti}(mat),solver)
end

Expand All @@ -122,6 +142,10 @@ Returns a PardisoNumericalSetup from a given AbstractSparseMatrix.
function build_PardisoNumericalSetup(phase::Integer,
mat::AbstractSparseMatrix{T,Ti},
solver::PardisoSolver) where {T,Ti}
if !(solver.mtype in (MTYPE_REAL_STRUCTURALLY_SYMMETRIC,
MTYPE_REAL_NON_SYMMETRIC))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
PardisoNumericalSetup(phase,SparseMatrixCSC{T,Ti}(mat),solver)
end

Expand All @@ -133,6 +157,10 @@ Returns a PardisoSymbolicSetup from a given SparseMatrixCSC.
function build_PardisoSymbolicSetup(phase::Integer,
mat::SparseMatrixCSC{T,Ti},
solver::PardisoSolver) where {T,Ti}
if !(solver.mtype in (MTYPE_REAL_STRUCTURALLY_SYMMETRIC,
MTYPE_REAL_NON_SYMMETRIC))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
PardisoSymbolicSetup(phase,mat,solver)
end

Expand All @@ -144,6 +172,10 @@ Returns a PardisoNumericalSetup from a given SparseMatrixCSC.
function build_PardisoNumericalSetup(phase::Integer,
mat::SparseMatrixCSC{T,Ti},
solver::PardisoSolver) where {T,Ti}
if !(solver.mtype in (MTYPE_REAL_STRUCTURALLY_SYMMETRIC,
MTYPE_REAL_NON_SYMMETRIC))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
PardisoNumericalSetup(phase,mat,solver)
end

Expand All @@ -155,6 +187,10 @@ Returns a PardisoSymbolicSetup from a given SparseMatrixCSR.
function build_PardisoSymbolicSetup(phase::Integer,
mat::SparseMatrixCSR{Bi},
solver::PardisoSolver) where {Bi}
if !(solver.mtype in (MTYPE_REAL_STRUCTURALLY_SYMMETRIC,
MTYPE_REAL_NON_SYMMETRIC))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING
PardisoSymbolicSetup(phase,mat,solver)
end
Expand All @@ -167,6 +203,10 @@ Returns a PardisoNumericalSetup from a given SparseMatrixCSR.
function build_PardisoNumericalSetup(phase::Integer,
mat::SparseMatrixCSR{Bi},
solver::PardisoSolver) where {Bi}
if !(solver.mtype in (MTYPE_REAL_STRUCTURALLY_SYMMETRIC,
MTYPE_REAL_NON_SYMMETRIC))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING
PardisoNumericalSetup(phase,mat,solver)
end
Expand All @@ -179,6 +219,10 @@ Returns a PardisoSymbolicSetup from a given SymSparseMatrixCSR.
function build_PardisoSymbolicSetup(phase::Integer,
mat::SymSparseMatrixCSR{Bi},
solver::PardisoSolver) where {Bi}
if !(solver.mtype in (MTYPE_REAL_SYMMETRIC_POSITIVE_DEFINITE,
MTYPE_REAL_SYMMETRIC_INDEFINITE))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING
PardisoSymbolicSetup(phase,mat.uppertrian,solver)
end
Expand All @@ -191,6 +235,10 @@ Returns a PardisoNumericalSetup from a given SparseMatrixCSR.
function build_PardisoNumericalSetup(phase::Integer,
mat::SymSparseMatrixCSR{Bi},
solver::PardisoSolver) where {Bi}
if !(solver.mtype in (MTYPE_REAL_SYMMETRIC_POSITIVE_DEFINITE,
MTYPE_REAL_SYMMETRIC_INDEFINITE))
@warn "Pardiso matrix type ($(solver.mtype)) does not match with $(typeof(mat))."
end
solver.iparm[IPARM_ONE_OR_ZERO_BASED_INDEXING] = Bi == 0 ? PARDISO_ZERO_BASED_INDEXING : PARDISO_ONE_BASED_INDEXING
PardisoNumericalSetup(phase,mat.uppertrian,solver)
end
Expand Down
32 changes: 28 additions & 4 deletions test/femdriver.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,13 @@ using GridapPardiso

tol = 1e-10

model = CartesianDiscreteModel((0,1,0,1,0,1), (10,10,10))
domain = (0,1,0,1,0,1)
partition = (10,10,10)

# Simple 2D data for debugging. TODO: remove when fixed.
domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

V = TestFESpace(reffe=:Lagrangian, order=1, valuetype=Float64,
conformity=:H1, model=model, dirichlet_tags="boundary")
Expand All @@ -21,16 +27,34 @@ t_Ω = AffineFETerm(
(v) -> inner(v, (x) -> x[1]*x[2] ),
trian, quad)

# With non-symmetric storage

op = AffineFEOperator(SparseMatrixCSR{1,Float64,Int},V,U,t_Ω)

ls = PardisoSolver()
ls = PardisoSolver(op)
solver = LinearFESolver(ls)

uh = solve(solver,op)

x = get_free_values(uh)
A = get_matrix(op)
b = get_vector(op)

r = A*x - b
@test maximum(abs.(r)) < tol

# With symmetric storage

op = AffineFEOperator(SymSparseMatrixCSR{1,Float64,Int},V,U,t_Ω)

ls = PardisoSolver(op)
solver = LinearFESolver(ls)

uh = solve(solver,op)

x = get_free_values(uh)
A = op.op.matrix
b = op.op.vector
A = get_matrix(op)
b = get_vector(op)

r = A*x - b
@test maximum(abs.(r)) < tol
Expand Down