NLSolve + GridapPardiso + SparseMatricesCSR error

[amartinhuertas]

Hi @fverdugo,

when running the following WE, I get the error below the example code. The error disappears if one uses SparseMatrixCSC instead of SparseMatrixCSR. I think the error is related to the NLsolver cache, which holds an object of type ::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}. Shouldn't Array{Float64,2} be SparseMatrixCSR. Do u know where the source of the error can be? (if not I will debug carefully).

@amartinhuertas

module PLaplacianTests

using Test
using Gridap
using GridapPardiso
import Gridap: ∇
using LinearAlgebra
using SparseMatricesCSR

domain = (0,1,0,1)
cells = (4,4)
model = CartesianDiscreteModel(domain,cells)

u(x) = x[1] + x[2]
∇u(x) = VectorValue(1,1)
f(x) = 0
∇(::typeof(u)) = ∇u

const p = 3
flux(∇u) = norm(∇u)^(p-2) * ∇u
dflux(∇du,∇u) = (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2)*∇du

order = 3
T = Float64

V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order),dirichlet_tags="boundary")
U = TrialFESpace(V,u)

Ω = Triangulation(model)

degree = 2*order
dΩ = Measure(Ω,degree)

res(u,v) = ∫( ∇(v)⋅(flux∘∇(u)) )*dΩ
jac(u,du,v) = ∫( ∇(v)⋅(dflux∘(∇(du),∇(u))) )*dΩ

assem = SparseMatrixAssembler(SparseMatrixCSR{1,Float64,Int},Vector{Float64},U,V)

op = FEOperator(res,jac,U,V,assem)

linear_solver=PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC,
                            GridapPardiso.new_iparm(),
                            GridapPardiso.MSGLVL_VERBOSE,
                            GridapPardiso.new_pardiso_handle())

nls = NLSolver(linear_solver;show_trace=false, method=:newton)
solver = FESolver(nls)

x = rand(T,num_free_dofs(U))
uh0 = FEFunction(U,x)
uh, = solve!(uh0,solver,op)
end # module

This function belongs to an interface definition and cannot be used.
Stacktrace:
 [1] newton_(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::LineSearches.Static, ::Gridap.Algebra.var"#linsolve!#4"{GridapPardiso.PardisoNumericalSetup{Float64,Int64}}, ::NLsolve.NewtonCache{Array{Float64,1}}) at /home/amartin/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:104
 [2] #newton#7 at /home/amartin/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:146 [inlined]
 [3] nlsolve(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::Gridap.Algebra.var"#linsolve!#4"{GridapPardiso.PardisoNumericalSetup{Float64,Int64}}, factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64) at /home/amartin/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:23
 [4] _nlsolve_with_updated_cache!(::Array{Float64,1}, ::Gridap.Algebra.NLSolver, ::Gridap.FESpaces.AlgebraicOpFromFEOp, ::Gridap.Algebra.NLSolversCache) at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/NLSolvers.jl:158
 [5] solve! at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/NLSolvers.jl:139 [inlined]
 [6] solve!(::Array{Float64,1}, ::Gridap.Algebra.NLSolver, ::Gridap.FESpaces.AlgebraicOpFromFEOp) at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/NonlinearSolvers.jl:22
 [7] solve!(::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{Gridap.CellData.ReferenceDomain}}, ::Gridap.FESpaces.NonlinearFESolver, ::Gridap.FESpaces.FEOperatorFromWeakForm, ::Nothing) at /home/amartin/.julia/packages/Gridap/J87be/src/FESpaces/FESolvers.jl:139
 [8] solve!(::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{Gridap.CellData.ReferenceDomain}}, ::Gridap.FESpaces.NonlinearFESolver, ::Gridap.FESpaces.FEOperatorFromWeakForm) at /home/amartin/.julia/packages/Gridap/J87be/src/FESpaces/FESolvers.jl:14
 [9] top-level scope at /home/amartin/git-repos/GridapPardiso.jl/mwe.jl:54
in expression starting at /home/amartin/git-repos/GridapPardiso.jl/mwe.jl:54

caused by: This function belongs to an interface definition and cannot be used.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] macro expansion at /home/amartin/.julia/packages/Gridap/J87be/src/Helpers/Macros.jl:9 [inlined]
 [3] numerical_setup!(::GridapPardiso.PardisoNumericalSetup{Float64,Int64}, ::Array{Float64,2}) at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/LinearSolvers.jl:121
 [4] (::Gridap.Algebra.var"#linsolve!#4"{GridapPardiso.PardisoNumericalSetup{Float64,Int64}})(::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/NLSolvers.jl:155
 [5] newton_(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::LineSearches.Static, ::Gridap.Algebra.var"#linsolve!#4"{GridapPardiso.PardisoNumericalSetup{Float64,Int64}}, ::NLsolve.NewtonCache{Array{Float64,1}}) at /home/amartin/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:94
 [6] #newton#7 at /home/amartin/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:146 [inlined]
 [7] nlsolve(::NLSolversBase.OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::Gridap.Algebra.var"#linsolve!#4"{GridapPardiso.PardisoNumericalSetup{Float64,Int64}}, factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64) at /home/amartin/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:23
 [8] _nlsolve_with_updated_cache!(::Array{Float64,1}, ::Gridap.Algebra.NLSolver, ::Gridap.FESpaces.AlgebraicOpFromFEOp, ::Gridap.Algebra.NLSolversCache) at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/NLSolvers.jl:158
 [9] solve! at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/NLSolvers.jl:139 [inlined]
 [10] solve!(::Array{Float64,1}, ::Gridap.Algebra.NLSolver, ::Gridap.FESpaces.AlgebraicOpFromFEOp) at /home/amartin/.julia/packages/Gridap/J87be/src/Algebra/NonlinearSolvers.jl:22
 [11] solve!(::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{Gridap.CellData.ReferenceDomain}}, ::Gridap.FESpaces.NonlinearFESolver, ::Gridap.FESpaces.FEOperatorFromWeakForm, ::Nothing) at /home/amartin/.julia/packages/Gridap/J87be/src/FESpaces/FESolvers.jl:139
 [12] solve!(::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{Gridap.CellData.ReferenceDomain}}, ::Gridap.FESpaces.NonlinearFESolver, ::Gridap.FESpaces.FEOperatorFromWeakForm) at /home/amartin/.julia/packages/Gridap/J87be/src/FESpaces/FESolvers.jl:14
 [13] top-level scope at /home/amartin/git-repos/GridapPardiso.jl/mwe.jl:54

package versions

(GridapPardiso) pkg> status
Project GridapPardiso v0.5.0
Status `~/git-repos/GridapPardiso.jl/Project.toml`
  [56d4f2e9] Gridap v0.16.5
  [a0a7dd2c] SparseMatricesCSR v0.6.1
  [8f399da3] Libdl
  [2f01184e] SparseArrays

[fverdugo]

I have never used Pardiso within a non-linear solver. Perhaps similar is not properly implemented for SparseMatrixCSR.

[amartinhuertas]

I have never used Pardiso within a non-linear solver.

Ok, I have now learnt that it is not actually a Pardiso issue. It is a SparseMatricesCSR issue, triggered from NLSolve.jl. In particular, in this line:

https://github.com/JuliaNLSolvers/NLSolversBase.jl/blob/master/src/objective_types/oncedifferentiable.jl#L231

copy(DF), with DF a SparseMatrixCSR returns a full matrix. For SparseMatrixCSR, copy is implemented here:

which(copy,(SparseMatrixCSR{1,Float64,Int64},))
copy(a::AbstractArray) in Base at abstractarray.jl:913

I guess that, as you say, we should overload similar for SparseMatrixCSR. This is dangerous even not using pardiso, as you end up factorizing a full matrix with backslash.

[amartinhuertas]

I just realized that, for SparseMatrixCSC, there is an specialized copy method

julia> which(copy,(SparseMatrixCSC{Float64,Int64},))
copy(S::SparseArrays.AbstractSparseMatrixCSC) in SparseArrays at /home/amartin/software_installers/julia-1.5.4/share/julia/stdlib/v1.5/SparseArrays/src/sparsematrix.jl:306

How do you think we should proceed to solve this issue? I guess that ideally one should somehow replicate what it is done in SparseMAtrixCSC for SparseMatrixCSR ?

[amartinhuertas]

ok, I have double checked that by adding this method

function Base.copy(a::SparseMatrixCSR{Bi}) where Bi
  SparseMatrixCSR{Bi}(a.m,a.n,copy(a.rowptr),copy(a.colval),copy(a.nzval))
end

The example above fully works. I will add this method to SparseMatrixCsr in a PR.

[amartinhuertas]

I will add this method to SparseMatrixCsr in a PR.

Done in PR gridap/SparseMatricesCSR.jl#12. Discussion can continue there. Clossing this issue.