Alberto,
You need to construct a matrix and pass it in as the second argument to
solver.setJacobianEquality(myEqualityJacobian, JE)
Barry
> On May 27, 2024, at 11:58 AM, Paganini, Alberto D.M. (Dr.)
> <a.pagan...@leicester.ac.uk> wrote:
>
> This Message Is From an External Sender
> This message came from outside your organization.
> Dear PETSc experts,
>
> I’m trying to set up a simple constrained optimization problem with
> petsc4py.TAO and I’m hitting a
>
> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
> probably memory access out of range
>
> Below is a minimal failing example (which also highlights how little I know
> about PETSc :) ). I hit this error because I was trying to understand how to
> assign values to the Jacobian of the equality constraint in the function
> myEqualityJacobian. Any advice on how to avoid the segmentation violation
> (and more generally, any advice about how to code this up better) are warmly
> welcomed.
>
> import petsc4py.PETSc as PETSc
> # initial guess
> x = PETSc.Vec().createSeq(2)
> x.set(0)
> x.assemble()
> print("initial guess: ", x.array_r)
> # create TAO solver and set solution
> solver = PETSc.TAO().create()
> solver.setType("almm")
> solver.setFromOptions()
> # add objective and gradient routines
> solver.setSolution(x)
> def myObjGrad(tao, x, g):
> J = x[0]**2 + x[1]**2
> dJ = 2*x
> dJ.copy(g)
> return J
> solver.setObjectiveGradient(myObjGrad, None)
> # add equality constraint and its derivative
> def myEqualityConstraint(tao, x, ce):
> ce.set((x[0]-2)**2 + x[1]**2)
> ce = PETSc.Vec().createSeq(1)
> ce.assemble()
> solver.setEqualityConstraints(myEqualityConstraint, ce)
> def myEqualityJacobian(tao, x, JE, JEpre):
> v = PETSc.Vec().createSeq(2)
> v.setValue(0, 4)
> v.assemble()
> vals = 2*x - v
> # How should I enter vals into JE?
> col = JE.getColumnVector(0)
> vals.copy(col)
> #col = JE.getDenseColumnVec(0)
> #JE.restoreDenseColumnVec(0)
> solver.setJacobianEquality(myEqualityJacobian)
> # solve the problem
> solver.solve()
> print("found optimum at: ", solver.getSolution().array_r)
> solver.destroy()
>
> Many thanks in advance for your help.
>
> Best,
> Alberto
