Many thanks Barry, that really helped! In my ignorance I was confused by
petsc4py’s signature
setJacobianEquality(self, jacobian_equality, Mat J=None, Mat P=None, args=None,
kargs=None)
which I misread as setting J was optional. Here is the complete working code
for the public archive in case someone else has the same question.
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 - 1)
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
JE.setValues(0, [0, 1], vals.array_r)
JE.assemble()
JE = PETSc.Mat().createDense([1,2])
JE.assemble()
solver.setJacobianEquality(myEqualityJacobian, JE)
# solve the problem
solver.solve()
print("found optimum at: ", solver.getSolution().array_r)
solver.destroy()
Best,
Alberto
From: Barry Smith <[email protected]>
Date: Monday, 27 May 2024 at 22:15
To: Paganini, Alberto D.M. (Dr.) <[email protected]>
Cc: [email protected] <[email protected]>
Subject: Re: [petsc-users] [11 SEGV: Segmentation Violation] petsc4py.TAO
You don't often get email from [email protected]. Learn why this is
important<https://urldefense.us/v3/__https://aka.ms/LearnAboutSenderIdentification__;!!G_uCfscf7eWS!cDFJhPD3KVsSrnNZRx_WfSLoEpdk6LP9YcsDXPd-sM98dT4l-Q6CIB9hUzRyBpdcJsKTS_PA6MC0nUKpAtMVp3oeDi-EyTnt$
>
***CAUTION:*** This email was sent from an EXTERNAL source. Think before
clicking links or opening attachments.
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.)
<[email protected]> 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
