DOLFIN

PETScKrylovMatrix multiplication in Python

Asked by Christian Clason

I am trying to solve a nonlinear pde-constrained optimization problem with a reduced Newton method, where instead of constructing the full Hessian, I compute the action of the Hessian on a function for use in a Krylov solver such as CG.

I assume that this can be done by constructing a PETScKrylovMatrix and defining a suitable mult function. (If there is a more high-level way of doing this, I'd be grateful for a hint.) However, I'm a bit stumped by the interface between PETSc and Dolfin here. The core step in the matrix multiplication would be something like

def mult(self, *args) :
    du = Function(V,args[0])
    L = du*v*dx
    problem = VariationalProblem(a, L, bc)
    y = problem.solve()
    args[1] = y.vector()

(assume that the function space V, test function v, boundary conditions bc and bilinear form a have all been defined previously). This however doesn't work, since args[0] and args[1] are not of proper type. How do I need to massage args[0] and args[1] to make this work? Any hints would be appreciated, especially a template for the use of PETScKrylovMatrix in Python.

(If I can get this working, I'd be happy to turn this into a documented demo. I'm sure there are more people interested in using FEniCS as a framework for pde-constrained optimization.)

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#1

You can use the slice operator [:] or the cpp.down_cast functions, depending
on what
args are.

Eg. you could do:

du = Function(V)
du.vector()[:] = args[0][:]

Hope this helps, or else let me know.

Kent

On 21 March 2011 18:12, Christian Clason <
<email address hidden>> wrote:

> New question #149919 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/149919
>
> I am trying to solve a nonlinear pde-constrained optimization problem with
> a reduced Newton method, where instead of constructing the full Hessian, I
> compute the action of the Hessian on a function for use in a Krylov solver
> such as CG.
>
> I assume that this can be done by constructing a PETScKrylovMatrix and
> defining a suitable mult function. (If there is a more high-level way of
> doing this, I'd be grateful for a hint.) However, I'm a bit stumped by the
> interface between PETSc and Dolfin here. The core step in the matrix
> multiplication would be something like
>
> def mult(self, *args) :
> du = Function(V,args[0])
> L = du*v*dx
> problem = VariationalProblem(a, L, bc)
> y = problem.solve()
> args[1] = y.vector()
>
> (assume that the function space V, test function v, boundary conditions bc
> and bilinear form a have all been defined previously). This however doesn't
> work, since args[0] and args[1] are not of proper type. How do I need to
> massage args[0] and args[1] to make this work? Any hints would be
> appreciated, especially a template for the use of PETScKrylovMatrix in
> Python.
>
> (If I can get this working, I'd be happy to turn this into a documented
> demo. I'm sure there are more people interested in using FEniCS as a
> framework for pde-constrained optimization.)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Christian Clason (christian-clason) said : 		#2

Thanks for your quick reply. I'm afraid this doesn't work:

   du.vector()[:] = args[0][:]
TypeError: 'SwigPyObject' object is not subscriptable

Unfortunately, I know next to nothing about swig.

Thanks,

Christian

Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#3

Could you do
print type(args[0])
print type(args|1])

(or show the complete code)

Kent

On 21 March 2011 18:49, Christian Clason <
<email address hidden>> wrote:

> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
>
> Status: Answered => Open
>
> Christian Clason is still having a problem:
> Thanks for your quick reply. I'm afraid this doesn't work:
>
> du.vector()[:] = args[0][:]
> TypeError: 'SwigPyObject' object is not subscriptable
>
> Unfortunately, I know next to nothing about swig.
>
> Thanks,
>
> Christian
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Christian Clason (christian-clason) said : 		#4

Sure:

<type 'SwigPyObject'>
<type 'SwigPyObject'>

I'll try to construct a minimal example.

Christian

Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#5

Are the args PETSc vectors created by PETSc ?

Kent

On 21 March 2011 18:53, Christian Clason <
<email address hidden>> wrote:

> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
>
> Christian Clason posted a new comment:
> Sure:
>
> <type 'SwigPyObject'>
> <type 'SwigPyObject'>
>
> I'll try to construct a minimal example.
>
> Christian
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Christian Clason (christian-clason) said : 		#6

Yes, via PETScKrylovSolver. Below is a minimal example for what I'm trying to do.

Christian

from dolfin import *

def boundary(x,on_boundary):
    return on_boundary

mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, 'CG', 1)

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx

class NewtonMatrix(PETScKrylovMatrix) :
    def __init__(self) :
        PETScKrylovMatrix.__init__(self, V.dim(), V.dim())

    def mult(self, *args) :
        du = Function(V)
        du.vector()[:] = args[0][:]
        f = du*v*dx
        problem = VariationalProblem(a, L, bc)
        dy = problem.solve()
        args[1][:] = dy.vector()[:]

H = NewtonMatrix()

E = project(Constant('1'),V)
b = PETScVector(E.vector().size())
b.set_local(E.vector().array())

x = PETScVector(E.vector().size())
x.zero()
NewtonSolver = PETScKrylovSolver('cg', 'none');
NewtonSolver.solve(H, x, b);

Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#7

Nice!
There may be some SWIG director methods missing for the matrices. I will
have a look
to check how to do it with or without these director methods.

Kent

On 21 March 2011 19:17, Christian Clason <
<email address hidden>> wrote:

> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
>
> Status: Answered => Open
>
> Christian Clason is still having a problem:
> Yes, via PETScKrylovSolver. Below is a minimal example for what I'm
> trying to do.
>
> Christian
>
> from dolfin import *
>
> def boundary(x,on_boundary):
> return on_boundary
>
> mesh = UnitSquare(32, 32)
> V = FunctionSpace(mesh, 'CG', 1)
>
> u0 = Constant(0.0)
> bc = DirichletBC(V, u0, boundary)
> u = TrialFunction(V)
> v = TestFunction(V)
> a = inner(grad(u), grad(v))*dx
>
> class NewtonMatrix(PETScKrylovMatrix) :
> def __init__(self) :
> PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
>
> def mult(self, *args) :
> du = Function(V)
> du.vector()[:] = args[0][:]
> f = du*v*dx
> problem = VariationalProblem(a, L, bc)
> dy = problem.solve()
> args[1][:] = dy.vector()[:]
>
> H = NewtonMatrix()
>
> E = project(Constant('1'),V)
> b = PETScVector(E.vector().size())
> b.set_local(E.vector().array())
>
> x = PETScVector(E.vector().size())
> x.zero()
> NewtonSolver = PETScKrylovSolver('cg', 'none');
> NewtonSolver.solve(H, x, b);
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Johan Hake (johan-hake) said : 		#8

On Monday March 21 2011 11:24:28 Kent-Andre Mardal wrote:
> Nice!
> There may be some SWIG director methods missing for the matrices. I will
> have a look
> to check how to do it with or without these director methods.

I think the director methods works fine. These are the guys that reports the
error.

But I think we need to change order of the #includes in dolfin_la.h so
PETScVector come before PETScKrylovSolver.

Can you check that Kent?

We need some sort of testing for this...

Johan

> Kent
>
> On 21 March 2011 19:17, Christian Clason <
>
> <email address hidden>> wrote:
> > Question #149919 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/149919
> >
> > Status: Answered => Open
> >
> > Christian Clason is still having a problem:
> > Yes, via PETScKrylovSolver. Below is a minimal example for what I'm
> > trying to do.
> >
> > Christian
> >
> > from dolfin import *
> >
> > def boundary(x,on_boundary):
> > return on_boundary
> >
> > mesh = UnitSquare(32, 32)
> > V = FunctionSpace(mesh, 'CG', 1)
> >
> > u0 = Constant(0.0)
> > bc = DirichletBC(V, u0, boundary)
> > u = TrialFunction(V)
> > v = TestFunction(V)
> > a = inner(grad(u), grad(v))*dx
> >
> > class NewtonMatrix(PETScKrylovMatrix) :
> > def __init__(self) :
> > PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
> >
> > def mult(self, *args) :
> > du = Function(V)
> > du.vector()[:] = args[0][:]
> >
> > f = du*v*dx
> >
> > problem = VariationalProblem(a, L, bc)
> >
> > dy = problem.solve()
> >
> > args[1][:] = dy.vector()[:]
> >
> > H = NewtonMatrix()
> >
> > E = project(Constant('1'),V)
> > b = PETScVector(E.vector().size())
> > b.set_local(E.vector().array())
> >
> > x = PETScVector(E.vector().size())
> > x.zero()
> > NewtonSolver = PETScKrylovSolver('cg', 'none');
> > NewtonSolver.solve(H, x, b);
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp

Revision history for this message
Garth Wells (garth-wells) said : 		#9

PETScKrylovMatrix is desperately untested. Some simple unit tests are needed.

Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#10

On 21 March 2011 19:33, Johan Hake <email address hidden>wrote:

> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
>
> Johan Hake proposed the following answer:
> On Monday March 21 2011 11:24:28 Kent-Andre Mardal wrote:
> > Nice!
> > There may be some SWIG director methods missing for the matrices. I will
> > have a look
> > to check how to do it with or without these director methods.
>
> I think the director methods works fine. These are the guys that reports
> the
> error.
>
> But I think we need to change order of the #includes in dolfin_la.h so
> PETScVector come before PETScKrylovSolver.
>
> Can you check that Kent?
>
> We need some sort of testing for this...
>
> Johan
>
>
Simply changing the order of includes did not fix it.

Kent

> > Kent
> >
> > On 21 March 2011 19:17, Christian Clason <
> >
> > <email address hidden>> wrote:
> > > Question #149919 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/149919
> > >
> > > Status: Answered => Open
> > >
> > > Christian Clason is still having a problem:
> > > Yes, via PETScKrylovSolver. Below is a minimal example for what I'm
> > > trying to do.
> > >
> > > Christian
> > >
> > > from dolfin import *
> > >
> > > def boundary(x,on_boundary):
> > > return on_boundary
> > >
> > > mesh = UnitSquare(32, 32)
> > > V = FunctionSpace(mesh, 'CG', 1)
> > >
> > > u0 = Constant(0.0)
> > > bc = DirichletBC(V, u0, boundary)
> > > u = TrialFunction(V)
> > > v = TestFunction(V)
> > > a = inner(grad(u), grad(v))*dx
> > >
> > > class NewtonMatrix(PETScKrylovMatrix) :
> > > def __init__(self) :
> > > PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
> > >
> > > def mult(self, *args) :
> > > du = Function(V)
> > > du.vector()[:] = args[0][:]
> > >
> > > f = du*v*dx
> > >
> > > problem = VariationalProblem(a, L, bc)
> > >
> > > dy = problem.solve()
> > >
> > > args[1][:] = dy.vector()[:]
> > >
> > > H = NewtonMatrix()
> > >
> > > E = project(Constant('1'),V)
> > > b = PETScVector(E.vector().size())
> > > b.set_local(E.vector().array())
> > >
> > > x = PETScVector(E.vector().size())
> > > x.zero()
> > > NewtonSolver = PETScKrylovSolver('cg', 'none');
> > > NewtonSolver.solve(H, x, b);
> > >
> > > --
> > > You received this question notification because you are a member of
> > > DOLFIN Team, which is an answer contact for DOLFIN.
> > >
> > > _______________________________________________
> > > Mailing list: https://launchpad.net/~dolfin
> > > Post to : <email address hidden>
> > > Unsubscribe : https://launchpad.net/~dolfin
> > > More help : https://help.launchpad.net/ListHelp
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Johan Hake (johan-hake) said : 		#11

On Monday March 21 2011 13:03:16 Kent-Andre Mardal wrote:
> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
>
> Kent-Andre Mardal proposed the following answer:
> On 21 March 2011 19:33, Johan Hake
>
> <email address hidden>wrote:
> > Question #149919 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/149919
> >
> > Johan Hake proposed the following answer:
> >
> > On Monday March 21 2011 11:24:28 Kent-Andre Mardal wrote:
> > > Nice!
> > > There may be some SWIG director methods missing for the matrices. I
> > > will have a look
> > > to check how to do it with or without these director methods.
> >
> > I think the director methods works fine. These are the guys that reports
> > the
> > error.
> >
> > But I think we need to change order of the #includes in dolfin_la.h so
> > PETScVector come before PETScKrylovSolver.
> >
> > Can you check that Kent?
> >
> > We need some sort of testing for this...
> >
> > Johan
>
> Simply changing the order of includes did not fix it.

Sorry fo the iterations. It should of course be that PETScVector needs to come
before PETScKrylovMatrix...

Not sure how we in general should address that SWIG is sensitive on
declaration order, as it can become quite tedious to resolve. But, as in C++
we can forward declare some types within foo_pre.i directly. In this way we
can just put all problematic types into something like:

// Forward declaration of types
namespace dolfin {
  PETScVector;
}

We already do this for some types. We could "standardize" this so we have such
a section in all foo_pre.i files.

Johan

> Kent
>
> > > Kent
> > >
> > > On 21 March 2011 19:17, Christian Clason <
> > >
> > > <email address hidden>> wrote:
> > > > Question #149919 on DOLFIN changed:
> > > > https://answers.launchpad.net/dolfin/+question/149919
> > > >
> > > > Status: Answered => Open
> > > >
> > > > Christian Clason is still having a problem:
> > > > Yes, via PETScKrylovSolver. Below is a minimal example for what I'm
> > > > trying to do.
> > > >
> > > > Christian
> > > >
> > > > from dolfin import *
> > > >
> > > > def boundary(x,on_boundary):
> > > > return on_boundary
> > > >
> > > > mesh = UnitSquare(32, 32)
> > > > V = FunctionSpace(mesh, 'CG', 1)
> > > >
> > > > u0 = Constant(0.0)
> > > > bc = DirichletBC(V, u0, boundary)
> > > > u = TrialFunction(V)
> > > > v = TestFunction(V)
> > > > a = inner(grad(u), grad(v))*dx
> > > >
> > > > class NewtonMatrix(PETScKrylovMatrix) :
> > > > def __init__(self) :
> > > > PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
> > > >
> > > > def mult(self, *args) :
> > > > du = Function(V)
> > > > du.vector()[:] = args[0][:]
> > > >
> > > > f = du*v*dx
> > > >
> > > > problem = VariationalProblem(a, L, bc)
> > > >
> > > > dy = problem.solve()
> > > >
> > > > args[1][:] = dy.vector()[:]
> > > >
> > > > H = NewtonMatrix()
> > > >
> > > > E = project(Constant('1'),V)
> > > > b = PETScVector(E.vector().size())
> > > > b.set_local(E.vector().array())
> > > >
> > > > x = PETScVector(E.vector().size())
> > > > x.zero()
> > > > NewtonSolver = PETScKrylovSolver('cg', 'none');
> > > > NewtonSolver.solve(H, x, b);
> > > >
> > > > --
> > > > You received this question notification because you are a member of
> > > > DOLFIN Team, which is an answer contact for DOLFIN.
> > > >
> > > > _______________________________________________
> > > > Mailing list: https://launchpad.net/~dolfin
> > > > Post to : <email address hidden>
> > > > Unsubscribe : https://launchpad.net/~dolfin
> > > > More help : https://help.launchpad.net/ListHelp
> >
> > --
> > You received this question notification because you are a member of
> > DOLFIN Team, which is an answer contact for DOLFIN.
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp

Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#12

On 21 March 2011 21:10, Johan Hake <email address hidden> wrote:

> On Monday March 21 2011 13:03:16 Kent-Andre Mardal wrote:
> > Question #149919 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/149919
> >
> > Kent-Andre Mardal proposed the following answer:
> > On 21 March 2011 19:33, Johan Hake
> >
> > <email address hidden>wrote:
> > > Question #149919 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/149919
> > >
> > > Johan Hake proposed the following answer:
> > >
> > > On Monday March 21 2011 11:24:28 Kent-Andre Mardal wrote:
> > > > Nice!
> > > > There may be some SWIG director methods missing for the matrices. I
> > > > will have a look
> > > > to check how to do it with or without these director methods.
> > >
> > > I think the director methods works fine. These are the guys that
> reports
> > > the
> > > error.
> > >
> > > But I think we need to change order of the #includes in dolfin_la.h so
> > > PETScVector come before PETScKrylovSolver.
> > >
> > > Can you check that Kent?
> > >
> > > We need some sort of testing for this...
> > >
> > > Johan
> >
> > Simply changing the order of includes did not fix it.
>
> Sorry fo the iterations. It should of course be that PETScVector needs to
> come
> before PETScKrylovMatrix...
>
>
I did (in kernel_modules) but it did not help.

> Not sure how we in general should address that SWIG is sensitive on
> declaration order, as it can become quite tedious to resolve. But, as in
> C++
> we can forward declare some types within foo_pre.i directly. In this way we
> can just put all problematic types into something like:
>
> // Forward declaration of types
> namespace dolfin {
> PETScVector;
> }
>
> We already do this for some types. We could "standardize" this so we have
> such
> a section in all foo_pre.i files.
>
> Johan

>
> > Kent
> >
> > > > Kent
> > > >
> > > > On 21 March 2011 19:17, Christian Clason <
> > > >
> > > > <email address hidden>> wrote:
> > > > > Question #149919 on DOLFIN changed:
> > > > > https://answers.launchpad.net/dolfin/+question/149919
> > > > >
> > > > > Status: Answered => Open
> > > > >
> > > > > Christian Clason is still having a problem:
> > > > > Yes, via PETScKrylovSolver. Below is a minimal example for what I'm
> > > > > trying to do.
> > > > >
> > > > > Christian
> > > > >
> > > > > from dolfin import *
> > > > >
> > > > > def boundary(x,on_boundary):
> > > > > return on_boundary
> > > > >
> > > > > mesh = UnitSquare(32, 32)
> > > > > V = FunctionSpace(mesh, 'CG', 1)
> > > > >
> > > > > u0 = Constant(0.0)
> > > > > bc = DirichletBC(V, u0, boundary)
> > > > > u = TrialFunction(V)
> > > > > v = TestFunction(V)
> > > > > a = inner(grad(u), grad(v))*dx
> > > > >
> > > > > class NewtonMatrix(PETScKrylovMatrix) :
> > > > > def __init__(self) :
> > > > > PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
> > > > >
> > > > > def mult(self, *args) :
> > > > > du = Function(V)
> > > > > du.vector()[:] = args[0][:]
> > > > >
> > > > > f = du*v*dx
> > > > >
> > > > > problem = VariationalProblem(a, L, bc)
> > > > >
> > > > > dy = problem.solve()
> > > > >
> > > > > args[1][:] = dy.vector()[:]
> > > > >
> > > > > H = NewtonMatrix()
> > > > >
> > > > > E = project(Constant('1'),V)
> > > > > b = PETScVector(E.vector().size())
> > > > > b.set_local(E.vector().array())
> > > > >
> > > > > x = PETScVector(E.vector().size())
> > > > > x.zero()
> > > > > NewtonSolver = PETScKrylovSolver('cg', 'none');
> > > > > NewtonSolver.solve(H, x, b);
> > > > >
> > > > > --
> > > > > You received this question notification because you are a member of
> > > > > DOLFIN Team, which is an answer contact for DOLFIN.
> > > > >
> > > > > _______________________________________________
> > > > > Mailing list: https://launchpad.net/~dolfin
> > > > > Post to : <email address hidden>
> > > > > Unsubscribe : https://launchpad.net/~dolfin
> > > > > More help : https://help.launchpad.net/ListHelp
> > >
> > > --
> > > You received this question notification because you are a member of
> > > DOLFIN Team, which is an answer contact for DOLFIN.
> > >
> > > _______________________________________________
> > > Mailing list: https://launchpad.net/~dolfin
> > > Post to : <email address hidden>
> > > Unsubscribe : https://launchpad.net/~dolfin
> > > More help : https://help.launchpad.net/ListHelp
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Johan Hake (johan-hake) said : 		#13

On Monday March 21 2011 14:08:35 Kent-Andre Mardal wrote:
> On 21 March 2011 21:10, Johan Hake <email address hidden> wrote:
> > On Monday March 21 2011 13:03:16 Kent-Andre Mardal wrote:
> > > Question #149919 on DOLFIN changed:
> > > https://answers.launchpad.net/dolfin/+question/149919
> > >
> > > Kent-Andre Mardal proposed the following answer:
> > > On 21 March 2011 19:33, Johan Hake
> > >
> > > <email address hidden>wrote:
> > > > Question #149919 on DOLFIN changed:
> > > > https://answers.launchpad.net/dolfin/+question/149919
> > > >
> > > > Johan Hake proposed the following answer:
> > > >
> > > > On Monday March 21 2011 11:24:28 Kent-Andre Mardal wrote:
> > > > > Nice!
> > > > > There may be some SWIG director methods missing for the matrices. I
> > > > > will have a look
> > > > > to check how to do it with or without these director methods.
> > > >
> > > > I think the director methods works fine. These are the guys that
> >
> > reports
> >
> > > > the
> > > > error.
> > > >
> > > > But I think we need to change order of the #includes in dolfin_la.h
> > > > so PETScVector come before PETScKrylovSolver.
> > > >
> > > > Can you check that Kent?
> > > >
> > > > We need some sort of testing for this...
> > > >
> > > > Johan
> > >
> > > Simply changing the order of includes did not fix it.
> >
> > Sorry fo the iterations. It should of course be that PETScVector needs to
> > come
> > before PETScKrylovMatrix...
>
> I did (in kernel_modules) but it did not help.

Then I am not sure what will work. Need to spend time debugging which I
unfortunately does not have for now. File a bug and attach the test script.

Johan

> > Not sure how we in general should address that SWIG is sensitive on
> > declaration order, as it can become quite tedious to resolve. But, as in
> > C++
> > we can forward declare some types within foo_pre.i directly. In this way
> > we can just put all problematic types into something like:
> >
> > // Forward declaration of types
> > namespace dolfin {
> >
> > PETScVector;
> >
> > }
> >
> > We already do this for some types. We could "standardize" this so we have
> > such
> > a section in all foo_pre.i files.
> >
> > Johan
> >
> > > Kent
> > >
> > > > > Kent
> > > > >
> > > > > On 21 March 2011 19:17, Christian Clason <
> > > > >
> > > > > <email address hidden>> wrote:
> > > > > > Question #149919 on DOLFIN changed:
> > > > > > https://answers.launchpad.net/dolfin/+question/149919
> > > > > >
> > > > > > Status: Answered => Open
> > > > > >
> > > > > > Christian Clason is still having a problem:
> > > > > > Yes, via PETScKrylovSolver. Below is a minimal example for what
> > > > > > I'm trying to do.
> > > > > >
> > > > > > Christian
> > > > > >
> > > > > > from dolfin import *
> > > > > >
> > > > > > def boundary(x,on_boundary):
> > > > > > return on_boundary
> > > > > >
> > > > > > mesh = UnitSquare(32, 32)
> > > > > > V = FunctionSpace(mesh, 'CG', 1)
> > > > > >
> > > > > > u0 = Constant(0.0)
> > > > > > bc = DirichletBC(V, u0, boundary)
> > > > > > u = TrialFunction(V)
> > > > > > v = TestFunction(V)
> > > > > > a = inner(grad(u), grad(v))*dx
> > > > > >
> > > > > > class NewtonMatrix(PETScKrylovMatrix) :
> > > > > > def __init__(self) :
> > > > > > PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
> > > > > >
> > > > > > def mult(self, *args) :
> > > > > > du = Function(V)
> > > > > > du.vector()[:] = args[0][:]
> > > > > >
> > > > > > f = du*v*dx
> > > > > >
> > > > > > problem = VariationalProblem(a, L, bc)
> > > > > >
> > > > > > dy = problem.solve()
> > > > > >
> > > > > > args[1][:] = dy.vector()[:]
> > > > > >
> > > > > > H = NewtonMatrix()
> > > > > >
> > > > > > E = project(Constant('1'),V)
> > > > > > b = PETScVector(E.vector().size())
> > > > > > b.set_local(E.vector().array())
> > > > > >
> > > > > > x = PETScVector(E.vector().size())
> > > > > > x.zero()
> > > > > > NewtonSolver = PETScKrylovSolver('cg', 'none');
> > > > > > NewtonSolver.solve(H, x, b);
> > > > > >
> > > > > > --
> > > > > > You received this question notification because you are a member
> > > > > > of DOLFIN Team, which is an answer contact for DOLFIN.
> > > > > >
> > > > > > _______________________________________________
> > > > > > Mailing list: https://launchpad.net/~dolfin
> > > > > > Post to : <email address hidden>
> > > > > > Unsubscribe : https://launchpad.net/~dolfin
> > > > > > More help : https://help.launchpad.net/ListHelp
> > > >
> > > > --
> > > > You received this question notification because you are a member of
> > > > DOLFIN Team, which is an answer contact for DOLFIN.
> > > >
> > > > _______________________________________________
> > > > Mailing list: https://launchpad.net/~dolfin
> > > > Post to : <email address hidden>
> > > > Unsubscribe : https://launchpad.net/~dolfin
> > > > More help : https://help.launchpad.net/ListHelp
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp

Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#14

Christian,

For now, you should write your Newton solver in Python. It shouldn´t be many
extra lines.
This would not require that cross language inheritance works properly for
matrices
and vectors, which there is no easy fix for at the moment.
Please tell me if this is not a feasible approach for you.

Kent

On 21 March 2011 18:12, Christian Clason <
<email address hidden>> wrote:

> New question #149919 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/149919
>
> I am trying to solve a nonlinear pde-constrained optimization problem with
> a reduced Newton method, where instead of constructing the full Hessian, I
> compute the action of the Hessian on a function for use in a Krylov solver
> such as CG.
>
> I assume that this can be done by constructing a PETScKrylovMatrix and
> defining a suitable mult function. (If there is a more high-level way of
> doing this, I'd be grateful for a hint.) However, I'm a bit stumped by the
> interface between PETSc and Dolfin here. The core step in the matrix
> multiplication would be something like
>
> def mult(self, *args) :
> du = Function(V,args[0])
> L = du*v*dx
> problem = VariationalProblem(a, L, bc)
> y = problem.solve()
> args[1] = y.vector()
>
> (assume that the function space V, test function v, boundary conditions bc
> and bilinear form a have all been defined previously). This however doesn't
> work, since args[0] and args[1] are not of proper type. How do I need to
> massage args[0] and args[1] to make this work? Any hints would be
> appreciated, especially a template for the use of PETScKrylovMatrix in
> Python.
>
> (If I can get this working, I'd be happy to turn this into a documented
> demo. I'm sure there are more people interested in using FEniCS as a
> framework for pde-constrained optimization.)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Christian Clason (christian-clason) said : 		#15

Kent,

since the problem I'm really interested in requires GMRES for the Newton solver, it's probably easier to switch to C++ for now. Since a Python version would be really useful as well, I'll file a bug.

Thanks for your efforts!

Christian

Revision history for this message
Christian Clason (christian-clason) said : 		#16

Just a small update: The problem seems to be indeed with PETScKrylovSolver, since PETScKrylovMatrix by itself works fine even with standard vectors:

from dolfin import *

def boundary(x,on_boundary):
    return on_boundary

mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, 'CG', 1)

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx

class NewtonMatrix(PETScKrylovMatrix) :
    def __init__(self) :
        PETScKrylovMatrix.__init__(self, V.dim(), V.dim())

    def mult(self, *args) :
        du = Function(V,args[0])
        f = du*v*dx
        problem = VariationalProblem(a, f, bc)
        dy = problem.solve()
        args[1][:] = dy.vector()[:]

H = NewtonMatrix()

y = Function(V)
b = project(Constant(1.0),V)

H.mult(b.vector(),y.vector())
print(y.vector().norm('l2'))

However,

NewtonSolver = PETScKrylovSolver('cg', 'none')
NewtonSolver.solve(H, b.vector(), y.vector())

throws the error

TypeError: in method 'PETScKrylovSolver_solve', argument 3 of type 'dolfin::PETScVector &'

Christian

Revision history for this message
Johan Hake (johan-hake) said : 		#17

On Tuesday March 22 2011 03:31:18 Christian Clason wrote:
> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
>
> Christian Clason posted a new comment:
> Just a small update: The problem seems to be indeed with
> PETScKrylovSolver, since PETScKrylovMatrix by itself works fine even
> with standard vectors:
>
> from dolfin import *
>
> def boundary(x,on_boundary):
> return on_boundary
>
> mesh = UnitSquare(32, 32)
> V = FunctionSpace(mesh, 'CG', 1)
>
> u0 = Constant(0.0)
> bc = DirichletBC(V, u0, boundary)
> u = TrialFunction(V)
> v = TestFunction(V)
> a = inner(grad(u), grad(v))*dx
>
> class NewtonMatrix(PETScKrylovMatrix) :
> def __init__(self) :
> PETScKrylovMatrix.__init__(self, V.dim(), V.dim())
>
> def mult(self, *args) :
> du = Function(V,args[0])
> f = du*v*dx
> problem = VariationalProblem(a, f, bc)
> dy = problem.solve()
> args[1][:] = dy.vector()[:]
>
> H = NewtonMatrix()
>
> y = Function(V)
> b = project(Constant(1.0),V)
>
> H.mult(b.vector(),y.vector())

This is expected to work, as everything is handled in Python. But when you use
the KrylovSolver we take a detour through C++. Your python mult function is
actually called from within C++, through a tricky callback machinery that SWIG
provides.

> print(y.vector().norm('l2'))
>
> However,
>
> NewtonSolver = PETScKrylovSolver('cg', 'none')
> NewtonSolver.solve(H, b.vector(), y.vector())
>
> throws the error
>
> TypeError: in method 'PETScKrylovSolver_solve', argument 3 of type
> 'dolfin::PETScVector &'

That is because the solve method that takes a PETScKrylovMatrix as first
argument needs a PETScVector and not a GenericVector as you provides above.

I really suspect that this is related to shared_ptr. I do not think the
shared_ptr functionality provides dereference typemaps for shared_ptr vectors.
It might be worth adding a shared_ptr version of all calls in
PETScKrylov{Matrix, Solve}.

Christian could you link to this thread when you file the bug report?

Johan

> Christian

Revision history for this message
Christian Clason (christian-clason) said : 		#18

Johan,

thank you for the explanation. That means I should call the PETScKrylovSolver with

y_petsc = PETScVector(V.dim()); y_petsc[:] = y.vector()[:]
x_petsc = PETScVector(V.dim())

NewtonSolver = PETScKrylovSolver('cg', 'none')
NewtonSolver.solve(H, x_petsc, y_petsc)

right? (I'm trying to get a correct minimal example for the bug report.)

Christian

Revision history for this message
Johan Hake (johan-hake) said : 		#19

On Tuesday March 22 2011 08:40:42 Christian Clason wrote:
> Question #149919 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/149919
>
> Christian Clason posted a new comment:
> Johan,
>
> thank you for the explanation. That means I should call the
> PETScKrylovSolver with
>
> y_petsc = PETScVector(V.dim()); y_petsc[:] = y.vector()[:]
> x_petsc = PETScVector(V.dim())

Sure or just:

  y_petsc = down_cast(y.vector())

Johan

> NewtonSolver = PETScKrylovSolver('cg', 'none')
> NewtonSolver.solve(H, x_petsc, y_petsc)
>
> right? (I'm trying to get a correct minimal example for the bug report.)
>
> Christian

Revision history for this message
Best Johan Hake (johan-hake) said : 		#20

This should be fixed now.

Johan

On Monday March 21 2011 10:12:52 Christian Clason wrote:
> New question #149919 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/149919
>
> I am trying to solve a nonlinear pde-constrained optimization problem with
> a reduced Newton method, where instead of constructing the full Hessian, I
> compute the action of the Hessian on a function for use in a Krylov solver
> such as CG.
>
> I assume that this can be done by constructing a PETScKrylovMatrix and
> defining a suitable mult function. (If there is a more high-level way of
> doing this, I'd be grateful for a hint.) However, I'm a bit stumped by the
> interface between PETSc and Dolfin here. The core step in the matrix
> multiplication would be something like
>
> def mult(self, *args) :
> du = Function(V,args[0])
> L = du*v*dx
> problem = VariationalProblem(a, L, bc)
> y = problem.solve()
> args[1] = y.vector()
>
> (assume that the function space V, test function v, boundary conditions bc
> and bilinear form a have all been defined previously). This however
> doesn't work, since args[0] and args[1] are not of proper type. How do I
> need to massage args[0] and args[1] to make this work? Any hints would be
> appreciated, especially a template for the use of PETScKrylovMatrix in
> Python.
>
> (If I can get this working, I'd be happy to turn this into a documented
> demo. I'm sure there are more people interested in using FEniCS as a
> framework for pde-constrained optimization.)

Revision history for this message
Kent-Andre Mardal (kent-and) said : 		#21

The SWIG wizard!

Kent

On 23 March 2011 22:01, Johan Hake <email address hidden> wrote:

> This should be fixed now.
>
> Johan
>
> On Monday March 21 2011 10:12:52 Christian Clason wrote:
> > New question #149919 on DOLFIN:
> > https://answers.launchpad.net/dolfin/+question/149919
> >
> > I am trying to solve a nonlinear pde-constrained optimization problem
> with
> > a reduced Newton method, where instead of constructing the full Hessian,
> I
> > compute the action of the Hessian on a function for use in a Krylov
> solver
> > such as CG.
> >
> > I assume that this can be done by constructing a PETScKrylovMatrix and
> > defining a suitable mult function. (If there is a more high-level way of
> > doing this, I'd be grateful for a hint.) However, I'm a bit stumped by
> the
> > interface between PETSc and Dolfin here. The core step in the matrix
> > multiplication would be something like
> >
> > def mult(self, *args) :
> > du = Function(V,args[0])
> > L = du*v*dx
> > problem = VariationalProblem(a, L, bc)
> > y = problem.solve()
> > args[1] = y.vector()
> >
> > (assume that the function space V, test function v, boundary conditions
> bc
> > and bilinear form a have all been defined previously). This however
> > doesn't work, since args[0] and args[1] are not of proper type. How do I
> > need to massage args[0] and args[1] to make this work? Any hints would be
> > appreciated, especially a template for the use of PETScKrylovMatrix in
> > Python.
> >
> > (If I can get this working, I'd be happy to turn this into a documented
> > demo. I'm sure there are more people interested in using FEniCS as a
> > framework for pde-constrained optimization.)
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
>

Revision history for this message
Johan Hake (johan-hake) said : 		#22

On Wednesday March 23 2011 14:20:16 Kent-Andre Mardal wrote:
> The SWIG wizard!

Well... copy/paste mostly :)

Johan

> Kent
>
> On 23 March 2011 22:01, Johan Hake <email address hidden> wrote:
> > This should be fixed now.
> >
> > Johan
> >
> > On Monday March 21 2011 10:12:52 Christian Clason wrote:
> > > New question #149919 on DOLFIN:
> > > https://answers.launchpad.net/dolfin/+question/149919
> > >
> > > I am trying to solve a nonlinear pde-constrained optimization problem
> >
> > with
> >
> > > a reduced Newton method, where instead of constructing the full
> > > Hessian,
> >
> > I
> >
> > > compute the action of the Hessian on a function for use in a Krylov
> >
> > solver
> >
> > > such as CG.
> > >
> > > I assume that this can be done by constructing a PETScKrylovMatrix and
> > > defining a suitable mult function. (If there is a more high-level way
> > > of doing this, I'd be grateful for a hint.) However, I'm a bit stumped
> > > by
> >
> > the
> >
> > > interface between PETSc and Dolfin here. The core step in the matrix
> > > multiplication would be something like
> > >
> > > def mult(self, *args) :
> > > du = Function(V,args[0])
> > > L = du*v*dx
> > > problem = VariationalProblem(a, L, bc)
> > > y = problem.solve()
> > > args[1] = y.vector()
> > >
> > > (assume that the function space V, test function v, boundary conditions
> >
> > bc
> >
> > > and bilinear form a have all been defined previously). This however
> > > doesn't work, since args[0] and args[1] are not of proper type. How do
> > > I need to massage args[0] and args[1] to make this work? Any hints
> > > would be appreciated, especially a template for the use of
> > > PETScKrylovMatrix in Python.
> > >
> > > (If I can get this working, I'd be happy to turn this into a documented
> > > demo. I'm sure there are more people interested in using FEniCS as a
> > > framework for pde-constrained optimization.)
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp

Revision history for this message
Anders Logg (logg) said : 		#23

On Wed, Mar 23, 2011 at 02:24:45PM -0700, Johan Hake wrote:
> On Wednesday March 23 2011 14:20:16 Kent-Andre Mardal wrote:
> > The SWIG wizard!
>
> Well... copy/paste mostly :)

But you copy/paste from your own magic SWIG code which still qualifies
and magic. :-)

--
Anders

> Johan
>
> > Kent
> >
> > On 23 March 2011 22:01, Johan Hake <email address hidden> wrote:
> > > This should be fixed now.
> > >
> > > Johan
> > >
> > > On Monday March 21 2011 10:12:52 Christian Clason wrote:
> > > > New question #149919 on DOLFIN:
> > > > https://answers.launchpad.net/dolfin/+question/149919
> > > >
> > > > I am trying to solve a nonlinear pde-constrained optimization problem
> > >
> > > with
> > >
> > > > a reduced Newton method, where instead of constructing the full
> > > > Hessian,
> > >
> > > I
> > >
> > > > compute the action of the Hessian on a function for use in a Krylov
> > >
> > > solver
> > >
> > > > such as CG.
> > > >
> > > > I assume that this can be done by constructing a PETScKrylovMatrix and
> > > > defining a suitable mult function. (If there is a more high-level way
> > > > of doing this, I'd be grateful for a hint.) However, I'm a bit stumped
> > > > by
> > >
> > > the
> > >
> > > > interface between PETSc and Dolfin here. The core step in the matrix
> > > > multiplication would be something like
> > > >
> > > > def mult(self, *args) :
> > > > du = Function(V,args[0])
> > > > L = du*v*dx
> > > > problem = VariationalProblem(a, L, bc)
> > > > y = problem.solve()
> > > > args[1] = y.vector()
> > > >
> > > > (assume that the function space V, test function v, boundary conditions
> > >
> > > bc
> > >
> > > > and bilinear form a have all been defined previously). This however
> > > > doesn't work, since args[0] and args[1] are not of proper type. How do
> > > > I need to massage args[0] and args[1] to make this work? Any hints
> > > > would be appreciated, especially a template for the use of
> > > > PETScKrylovMatrix in Python.
> > > >
> > > > (If I can get this working, I'd be happy to turn this into a documented
> > > > demo. I'm sure there are more people interested in using FEniCS as a
> > > > framework for pde-constrained optimization.)
> > >
> > > _______________________________________________
> > > Mailing list: https://launchpad.net/~dolfin
> > > Post to : <email address hidden>
> > > Unsubscribe : https://launchpad.net/~dolfin
> > > More help : https://help.launchpad.net/ListHelp
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

Revision history for this message
Christian Clason (christian-clason) said : 		#24

It certainly looks like magic to me :)

Thank you very much for fixing this so quickly!

Christian

Revision history for this message
Christian Clason (christian-clason) said : 		#25

Thanks Johan Hake, that solved my question.