Bound Constrained Linear Solver

On 12 March 2012 17:41, corrado maurini <
<email address hidden>> wrote:

> New question #190450 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/190450
>
> Hi,
>
> As I declared in a related question (
> https://answers.launchpad.net/dolfin/+question/170426), I need solvers
> able to deal with bound constraints for some application to fracture and
> damage mechanics.
>
> Finally, I found the time to work on the subject and I implemented an
> interface in dolfin for TAO, a library of large scale optimization based
> on petsc (http://www.mcs.anl.gov/research/projects/tao/).
>
> Taking as example the interface for the PETSc and SLEPc solvers available
> in dolfin/la, I implemented I C++ class to interface with TAO in order to
> solve linear bound constrained problem, i.e. problem of the type
>
> A x =b, with xl < = x <= xu.
>
> You may have a look to a preliminary version of the source here:
> -
> http://www.lmm.jussieu.fr/~corrado/tmp/tao_dolfin/TAOLinearBoundSolver.h
> -
> http://www.lmm.jussieu.fr/~corrado/tmp/tao_dolfin/TAOLinearBoundSolver.cpp
>
> I enjoyed very much the automatic generation of the python interface
> within dolfin.
> Now I have some example (poisson with constraints or linear elasticity
> with constraints) running correctly both in cpp and python.
>
>
Nice! Is this an interface to the complete TAO API or is it just the
functions you needed?
It could also be a FEniCS app, see
http://fenicsproject.org/applications/

Kent

> To have an idea of how the interface is working, see the script of the two
> examples that I tested (do not try to run, they need some hacking of the
> dolfin installation to have them working):
> http://www.lmm.jussieu.fr/~corrado/tmp/tao_dolfin/demo_poisson.py
> http://www.lmm.jussieu.fr/~corrado/tmp/tao_dolfin/demo_LinearElasticity.py
>
> I wonder if you are interested to include this interface in dolfin.
>
> Note that such solvers may be useful for several applications, for example
> contact mechanics, but they are generally missing in many finite element
> applications. Eventually I may provide my programs and help improving them,
> defining test cases and examples, etc ... For me it would be useful to have
> this directly in the dolfin distribution instead of manually setting it on
> my machines each time. I am very far from being an expert of C++
> programming, so there may be many things that should be improved.
>
> The compilation and linking of TAO (ver 2.0) with dolfin is quite easy.
> It works exactly as SLEPc, taking the configure generated by PETSc. I have
> written a FindTao.cmake, taking inspiration from FindSLEPc.cmake.
> You find the file at this link:
> http://www.lmm.jussieu.fr/~corrado/tmp/tao_dolfin/FindTAO.cmake
>
> Best,
>
> Corrado
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

This is an interface only for the solution of linear problems (i.e. quadratic functional) with bound constraints.
I think that the most useful TAO functions are included in this context. Additions of others is trivial.
For example, one should add options for the fine-tuning of the underlying PETSC-KSP used by TAO, which is set by default now.

TAO has also algorithms to solve non-linear problem with non-linear constraints (at least in the convex case).
I did not consider the interface to those solvers. It becomes more complex because it needs user-provided routines to assemble gradient and hessian (like for example a PETSC-SNES, for which the interface is not implemented in dolfin).

In the spirit of the current interface to solvers in FEniCS, non-linear problems may be solved by writing in FEniCS a Newton algorithm using external linear solvers (but

Now it is just a class (and few lines of code), so perhaps quite limited as an FEniCS app.
At the moment, I put this inside my dolfin installation (in the "la" folder) to take advantage of the automatic generation of the python interface from C++ code. I have no idea how to do it outside doflin.

The above comments are very important. PetSc 3.2 has solvers for variational inequalities as well as a host of improvements on the Non linear solvers.

I have been trying something along the lines mentioned in the comment above. However I am having problems interfacing with PetSc snes routines. In particular I am having trouble passing Jacobian routines. I have tried something like the following code snippets:

My solve function ...

...solve( ..., PETScVector& x, ....){
...
...
  PETScVector b1;
  PETScMatrix A1;
  nonlinear_function.J(A1,x); \\ routine that assembles forms and computes the Jacobian Matrix
  nonlinear_function.F(b1, x); \\ routine that assembles forms and computes the Residual Vector

  SNESSetFunction(snes,*b1.vec(),FormFunction,&nonlinear_function);
  SNESSetJacobian(snes,*A1.mat(),*A1.mat(),FormJacobian,&nonlinear_function);
  SNESSetFromOptions(snes);
  SNESSolve(snes,PETSC_NULL,*x.vec());
 ...
 ...
}

The form Jacobian function looks like below where I have used the void pointer to send the object that computes & assemples forms, computes residuals, Jacobians etc...

PetscErrorCode FormJacobian(SNES snes, Vec x, Mat *jac, Mat *B, MatStructure *flag, void *ptr)
{
 NonlinearProblem* nlp = (NonlinearProblem*)ptr;
 PETScMatrix A1;
 boost::shared_ptr<Vec> myxjptr(&x, NoDeleter());
 PETScVector sol(myxjptr);
 nlp->J(A1,sol);
 *jac = *A1.mat();
 myxjptr.reset();
 return 0;
}

This is not working. Although I can see that Jacobians are computed correctly and stored in the *jac. The snes solver does not accept it. I think it does not accept the earlier declaration in solve. Whenever the matrix needs to be accessed, it sends an error message that looks like the following.

[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Invalid argument!
[0]PETSC ERROR: Wrong type of object: Parameter # 2!
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Release Version 3.2.0, Patch 6, Wed Jan 11 09:28:45 CST 2012
...
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: KSPSetOperators() line 438 in src/ksp/ksp/interface/itcreate.c
[0]PETSC ERROR: SNESSetFromOptions() line 506 in src/snes/interface/snes.c
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Invalid argument!
[0]PETSC ERROR: Wrong type of object: Parameter # 2!
[0]PETSC ERROR: ------------------------------------------------------------------------
.....
[0]PETSC ERROR: KSPSetOperators() line 438 in src/ksp/ksp/interface/itcreate.c
[0]PETSC ERROR: SNESSolve_LS() line 189 in src/snes/impls/ls/ls.c
[0]PETSC ERROR: SNESSolve() line 2676 in src/snes/interface/snes.c

At the end I can destroy the *A1.mat() by calling MatDestroy&(*A1.mat()) but when the solver tries to destroy the Jacobian matrix, it cannot find it either:

[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Invalid argument!
[0]PETSC ERROR: Wrong type of object: Parameter # 1!
[0]PETSC ERROR: ------------------------------------------------------------------------
...
...
[0]PETSC ERROR: MatDestroy() line 1028 in src/mat/interface/matrix.c
[0]PETSC ERROR: SNESSetJacobian() line 1491 in src/snes/interface/snes.c

Why is the Jacobian matrix set in SNESSetJacobian and computed in FormJacobian not recognized by the petsc snes solver? Any insights or solutions would be greatly appreciated.

Have you tried declaring A1 as

boost::shared_ptr<const PETScMatrix> A1

?

This is the way it works in my code.

Corrado

Thanks for posting your code. I am trying to understand how it works. I think one should use the boost shared pointer as you suggest.
Meanwhile are you aware of FENICS being used for problems involving variational inequalities ?

Regards,
Rochan

I think I am the first one trying to do this.

My interface to TAO works in C++ and python for linear problem with bound constraints.

In a pure C++ code it should be straightforward to use also the non-linear contrained solvers of TAO.

If you need further information, feel free to contact me at <email address hidden>

For the Fenics developers, have you any plan to add variational inequality solvers to FEnics?

Corrado

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

I guess the answer is no. If you do not mind, please add a Blueprint and you can link this Question to it. Then it is easier to keep track of good ideas.

Ok. I will do.
Do you mind if I open a branch in launchpad with this feature added?
Currently I have this working in a personal repository.

On Wed, Apr 04, 2012 at 11:10:53AM -0000, corrado maurini wrote:
> Question #190450 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/190450
>
> Status: Expired => Open
>
> corrado maurini is still having a problem:
> Ok. I will do.
> Do you mind if I open a branch in launchpad with this feature added?
> Currently I have this working in a personal repository.

Yes, please do.

If/when you think your code is ready for merging, do a merge request.
I think this would be nice to include in DOLFIN. We already have
SLEPcEigenSolver so I don't think there is a problem if it is PETSc
specific.

I haven't looked closely at your code yet but a merge would need to
include: the interface, Python wrappers, a demo and a unit test.

--
Anders

Blueprint added: https://blueprints.launchpad.net/dolfin/+spec/tao-bound-constrained-solver

I will open a branch for this on launchpad

Corrado