Constrained BlockLinearOperator

[Bruno Turcksin]

Hi all,

I want to use BlockLinearOperator for one of my code but I have two questions:
  - there is a constrained_linear_operator function for LinearOperator but none for BlockLinearOperator. Would there be any difficulty to implement it?
  - for the constrained_linear_operator, the documentation says: "Currently, this function may not work correctly for distributed data structures." What's the problem with distributed data structures?

Best,

Bruno

[luca.heltai]

Ciao Bruno,

BlockLinearOperator is a LinearOperator, so the contrained_linear_operator will work just fine.

The issue with distributed data structure is that when we implemented it, we were not able to make it scale very well, and we stopped its development for distributed data structures. If I recall correctly, there have been some changes in the code, but I’m not sure these changes were done with distributed data structures in mind.

You can try it out, by wrapping the matrix in step-40 with a linear operator and see how far you get…

L.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Martin Kronbichler]

Hi Bruno,

> - for the constrained_linear_operator, the documentation says:
> "Currently, this function may not work correctly for distributed data
> structures." What's the problem with distributed data structures?

This comment was requested by me when reviewing the patch I think. As
far as I know, the code uses distributed_constraints_linear_operator in
the vmult_add and Tvmult_add versions which implements the matrix-vector
product by explicitly going through the 'get_constraint_entries' of the
ConstraintMatrix object. However, the loop that goes through the
constraints usually involves remote vector entries from other
processors. Unless you have a vector that has imported the appropriate
ghosts, it will lead to invalid access. (Invalid read for vmult_add and
invalid write for Tvmult_add.)

We should probably fix this at some point but I don't see a trivial
solution (we could use the ConstraintMatrix::distribute/condense
operations that implement these operations in a distributed sense but
they create+copy some temporary vectors which is quite expensive
efficient). Looking at the code we have now, I wonder how this has
worked in parallel at all, though. Luca, did you have different code
when evaluating these things or were you just lucky to not have any
constraints that cross processor boundaries?

Best,
Martin

[Bruno Turcksin]

Martin,

This comment was requested by me when reviewing the patch I think. As far as I know, the code uses distributed_constraints_linear_operator in the vmult_add and Tvmult_add versions which implements the matrix-vector product by explicitly going through the 'get_constraint_entries' of the ConstraintMatrix object. However, the loop that goes through the constraints usually involves remote vector entries from other processors. Unless you have a vector that has imported the appropriate ghosts, it will lead to invalid access. (Invalid read for vmult_add and invalid write for Tvmult_add.)

Could I solve this problem by using a distributed matrix to store the constraints instead of a ConstraintMatrix? The matrix would be very sparse but it should scale.

Best,

Bruno

[Martin Kronbichler]

Hi Bruno,

> Could I solve this problem by using a distributed matrix to store the
constraints instead of a ConstraintMatrix? The matrix would be very
sparse but it should scale.

I think that would be perfectly viable. However, I have not really clear
for me how to construct/keep that matrix. It's going to require some
work. You would need to attach a distributed matrix that fits your
vector (i.e., identify the right matrix template for a vector template
if it ought to be generic), fill it with ones on the diagonal for the
unconstrained case or the constraint_entries in the constrained case,
and finally pass it to the correct linear operator object for accessing
the respective vmult() method. And ideally keep the matrix between
different calls.

Best,
Martin