Add support for constrained systems

[MarcelKoch]

This PR adds support for linear systems Ax=b with constraints of the form x[i] = g[i], i in Idxs. Since the PR already got somewhat large, I've only implemented reference kernels. Additionally, an example which shows the constraints handling in action is supplied.

The usage of the handler looks like this:

ConstraintsHandler handler(idxs, A, values, b, /*optional*/ x, /*optional*/ std::make_shared<ZeroRowStrategy>());

auto cons_A = handler.get_operator();  // cons_A is a shared_ptr
auto cons_b = handler.get_right_hand_side();  //cons_b is a bare ptr
auto cons_x = handler.get_get_initial_guess();  //cons_x is a bare_ptr

auto solver = factory->generate(lend(cons_A));
solver->apply(cons_b, cons_x);
handler.correct_solution(cons_x);

Alternatively, the handler can be constructed just from idxs, A only, which requires setting at least the right-hand-side and the constrained values using

handler.with_constrained_values(values)
  .with_right_hand_side(b)
  .with_initial_guess(x);  // the last one is optional

With this approach the constrained vectors are only created when queried from get_*. Also, this allows for repeatedly solving the system with the same operator and the same index set, but with changing right-hand-sides, constrained values or both. An use-case for this could be time-dependent simulations.

If the initial guess has not been specified, then a default one, filled with zero and constrained values, is created.

Some topics that might be discussed:

• Should the ConstraintsHandler use pointer or value semantics? By that I mean, if it should use a create method similarly to many of ginkgo's classes.
• If pointer semantics are used, should the class hold an executor to use the EnableCreateMethod? Currently, that would be pointless, since I'm always using the executors of the vectors or matrix.
• I'm throwing an NotSupported error, if the get_* methods are used while the handler is in an incomplete state. Should a more descriptive exception (e.g. InvalidState) be added?
• Ideally, I would like to change the type of the constraint indices to the index set from Add an IndexSet class. #676
• How should different matrix formats be supported? Currently, only CSR is supported, and some implementations are really specific to it, since they access the underlying arrays.

---

Old Description

This PR will add support for linear systems Ax=b with constraints of the form x[i] = g[i], i in Idxs. Currently, it only contains a first interface for a constraint handler, and acts more as a base for future discussions.

The usage of the handler looks lilke this:

ConstrainedHandler cons_handler(idxs, values);

cons_handler.setup(A, b, x);

auto cons_A = cons_handler.get_operator();
auto cons_b = cons_handler.get_right_hand_side();
auto cons_x = cons_handler.get_get_initial_guess();

auto solver = factory->generate(lend(cons_A));
solver->apply(lend(cons_b), lend(cons_x));

I've also added setters, which I called update_* to (hopefully) indicate that they do more than simply setting a member.

I'm not sure if this interface will become abstract, or if one switches between different behaviors with a tag.

I'm also thinking about adding something like ConstrainedSolver which combines the handler and a solver for easier usage.

[sonarqubecloud]

Kudos, SonarCloud Quality Gate passed!   

0 Bugs
0 Vulnerabilities
0 Security Hotspots
0 Code Smells

No Coverage information
0.0% Duplication

[fritzgoebel]

This looks very useful! The interface looks quite similar to what is there in opencarp for handling Dirichlet boundary conditions. They have a dbc_manager which is tied to a matrix and some stimuli that are either active or inactive. Depending on which stimuli are active, the boundary conditions are updated with an update function. It looks very close to what you're suggesting.

[MarcelKoch]

I've been thinking a bit about what the invariants of the ConstraintHandler should be. The current interface would suggest that only the constraint indices are invariants. But I think that is too general. I would suggest that both the indices and the system matrix are invariant.

For the right-hand-side, the initial guess, and the constrained values I can think of scenarios, where they change, while the system matrix and the indices remain constant.
For the matrix, I don't think there is a reasonable scenario where only the system matrix changes, while the rest stays the same. Additionally, changing the system matrix requires recomputing everything else, so one might as well just create a new handler.

[MarcelKoch]

I've updated the interface according to the previous comment.

Here is an example which shows how the new interface could be used.

ConstraintHandler ch(idxs, A);

auto cons_A = ch.get_operator();
auto solver = factory->generate(const_A);

for(int time_step=0; time_step < end; ++ time_step){
  // compute new right-hand-side, dirichlet values, etc.

  ch.with_constrained_values(dir_vals)
    .with_right_hand_side(rhs);

  auto cons_rhs = ch.get_right_hand_side();
  auto const_x0 = ch.get_initial_guess();

  solver->apply(lend(cons_rhs), lend(cons_x0));

  ch.correct_solution(lend(cons_x0));
}

[MarcelKoch]

I'm currently a bit unsure how to integrate this into ginkgo's class hierarchy. Using the EnablePolymorphicObject mixin seems reasonable, with the objection that this would inject an executor into the class. But I don't think that an executor is useful for this interface. The only use I can see is to check that it's using the same executor as the operator, indices, etc. Other than that, I don't know what to do with it. Nevertheless, I would be nice to have the same semantics (i.e. pointer semantic) as the rest of ginkgo.

[MarcelKoch]

I think the interface is now where I'm satisfied with it. @fritzgoebel, @greole do you think you can use that in your implementations?

[fritzgoebel]

I think the interface is now where I'm satisfied with it. @fritzgoebel, @greole do you think you can use that in your implementations?

Yes, I think I will be able to use this.
Also: I agree that we can assume the system matrix to be invariant. The constraint handler looks like it will be quite light-weight, so when a new system matrix has to be computed creating a new constraint handler won't matter too much.

[MarcelKoch]

@fritzgoebel, @greole would you mind looking into this, if it could be used in your cases?

[greole]

That example looks broken, maybe it is just the missing line break between matrix); and handler . But also the use of the factory parameters like .with_constrained_values() etc on the already existing instance of the ConstraintHandler are somewhat confusing to me.

[greole]

If this changes the state of the underlying system, is there any function to get whether the underlying system is reconstructed, ie bool is_reconstructed();

[MarcelKoch]

When this is called, all internal vectors are reset, so it forces the recomputation of the constrained system (except for the operator).

[greole]

Is there a specific reason to return a raw pointer vs std::shared_ptr<const Dense> here?

[MarcelKoch]

Mostly to be more inline with how you would pass these vector along to other functions, e.g. apply. Also I don't think there is any need to use shared_ptr here, as the lifetime of these objects is always handled by this class.

[ginkgo-bot]

Error: The following files need to be formatted:

core/constraints/utility_kernels.hpp
include/ginkgo/ginkgo.hpp
reference/test/constraints/constraints_handler.cpp
reference/test/constraints/zero_rows.cpp

You can find a formatting patch under Artifacts here or run format! if you have write access to Ginkgo

[ginkgo-bot]

Note: This PR changes the Ginkgo ABI:

Functions changes summary: 106 Removed, 2 Changed, 524 Added functions
Variables changes summary: 0 Removed, 0 Changed, 0 Added variable

For details check the full ABI diff under Artifacts here

[MarcelKoch]

I'm closing this PR, since I will revisit this with a reduced approach. Instead of providing a single interface that does everything, only the underlying building blocks, necessary to enforce the constraints will be added.