[petsc-users] [11 SEGV: Segmentation Violation] petsc4py.TAO

[Paganini, Alberto D.M. (Dr.)]

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 <bsmith at petsc.dev>
Date: Monday, 27 May 2024 at 22:15
To: Paganini, Alberto D.M. (Dr.) <a.paganini at leicester.ac.uk>
Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
Subject: Re: [petsc-users] [11 SEGV: Segmentation Violation] petsc4py.TAO
You don't often get email from bsmith at petsc.dev. 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.) <a.paganini at 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20240527/b50dc130/attachment-0001.html>