LinearOperator of BlockVector issue and potential fix

[Tao Jin]

Dear all,

Background: I have a BlockVector<double> u that needs to be updated via a series of rank-2 update to get BlockVector<double> v. For notation convenience, let's call this series of rank-2 update as P. I can define P as either a LinearOperator or BlockLinearOperator, since I need to solve a linear system later that involves P using an iterative solver (for example, conjugate-gradient method). Defining P as a BlockLinearOperator does not work in this case, since I only know how to perform the rank-2 updates on u but not the block structure of the operator P. 

Solution: My solution is to define P as a LinearOperator that operates on the source BlockVector u to get the destination BlockVector v. I know this is unconventional, because LinearOperator usually operates on Vector but not BlockVector. But theoretically, a LinearOperator should still be able to operator on a BlockVector to produce another BlockVector.

Issue: I can define the needed LinearOperator P1 and P2 that operate on a BlockVector without a problem. I can call P1.vmult(v, u) and P2.vmult(v, u) correctly. However, when I define a composite operator such as (for demonstration purpose)

auto total_op = P1 * P2; 

and then call

total_op.vmult(v, u);

I get the following error message:

"""

--------------------------------------------------------
An error occurred in line <1487> of file </home/taojin/dealii-9.4.0/bin/include/deal.II/lac/block_vector_base.h> in function
    dealii::BlockVectorBase<VectorType>::BlockType& dealii::BlockVectorBase<VectorType>::block(unsigned int) [with VectorType = dealii::Vector<double>; dealii::BlockVectorBase<VectorType>::BlockType = dealii::Vector<double>]
The violated condition was:
    ::dealii::deal_II_exceptions::internals::compare_less_than(i, n_blocks())
Additional information:
    Index 0 is not in the half-open range [0,0). In the current case, this
    half-open range is in fact empty, suggesting that you are accessing an
    element of an empty collection such as a vector that has not been set
    to the correct size.

"""

Fix: This error message is due to the definition of operator* (Line 582 in linear_operator.h)

template <typename Range,
          typename Intermediate,
          typename Domain,
          typename Payload>
LinearOperator<Range, Domain, Payload>
operator*(const LinearOperator<Range, Intermediate, Payload> & first_op,
          const LinearOperator<Intermediate, Domain, Payload> &second_op)
{

......

      return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
        GrowingVectorMemory<Intermediate> vector_memory;

        typename VectorMemory<Intermediate>::Pointer i(vector_memory);
        second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
        second_op.vmult(*i, u);
        first_op.vmult(v, *i);

......

}

It looks like that Intermediate only has the correct memory size but not the block information. My fix is to add one line of code 

        (*i).reinit(u.n_blocks());

in my own definition operator*(). This fixes the error and works correctly.

Questions:

1. Does it even make sense to define a LinearOperator operating on a BlockVector? If yes, should we expand the dealii capability to consider this option?

2. For my rank-2 update operation P, I cannot define it as a BlockLinearOperator. Even though P operates on a BlockVector, but  P itself does not possess blocks. Any thoughts?

3. Another solution is to convert BlockVector u to a Vector and define  P as a LinearOperator operating on Vector. After the rank-2 update, convert the Vector back to a BlockVector. However, it seems rather awkward. Does the community have better ideas?

I have attached my code and a test case here.

Best,

Tao

[Tao Jin]

Dear all,

The actual fix to this issue is much simpler than what I originally thought. I just need to make sure the destination BlockVector v has the proper block sizes when I define the vmult(),

vmult_add(),  

Tvmult(), 

Tvmult_add(), 

reinit_range_vector,

and reinit_domain_vector for the LinearOperator.

We don't need to touch the definition of operator*() at all. I attach the following test case as an example, in which I use two FullMatrix objects to construct a LinearOperator that operates on BlockVector.

  using namespace dealii;

  // LinearOperator
  using Payload = dealii::internal::LinearOperatorImplementation::EmptyPayload;
  LinearOperator<BlockVector<double>, BlockVector<double>, Payload> usr_linearOperator_1;

  // For demonstration purpose, define LinearOperator using two FullMatrix
  FullMatrix<double> usr_matrix_1(3, 3);
  usr_matrix_1(0,0) = 1;
  usr_matrix_1(0,1) = 2;
  usr_matrix_1(0,2) = 3;

  usr_matrix_1(1,0) = 4;
  usr_matrix_1(1,1) = 5;
  usr_matrix_1(1,2) = 6;

  usr_matrix_1(2,0) = 7;
  usr_matrix_1(2,1) = 8;
  usr_matrix_1(2,2) = 9;

  FullMatrix<double> usr_matrix_2(2, 2);
  usr_matrix_2(0,0) = -1;
  usr_matrix_2(0,1) = -2;

  usr_matrix_2(1,0) = -3;
  usr_matrix_2(1,1) = -4;

  usr_linearOperator_1.vmult = [& usr_matrix_1, & usr_matrix_2](BlockVector<double> &v, const BlockVector<double> &u) {
    std::vector<types::global_dof_index> sizes_per_block(u.n_blocks());
    for (unsigned int i = 0; i < u.n_blocks(); ++i)
      sizes_per_block[i] = u.block(i).size();
    v.reinit(sizes_per_block);

    usr_matrix_1.vmult(v.block(0), u.block(0));
    usr_matrix_2.vmult(v.block(1), u.block(1));
  };

  usr_linearOperator_1.vmult_add = [&usr_linearOperator_1](BlockVector<double> &v, const BlockVector<double> &u) {
    std::vector<types::global_dof_index> sizes_per_block(u.n_blocks());
    for (unsigned int i = 0; i < u.n_blocks(); ++i)
      sizes_per_block[i] = u.block(i).size();

    BlockVector<double> temp_v(sizes_per_block);
    usr_linearOperator_1.vmult(temp_v, u);
    v += temp_v;
  };

  usr_linearOperator_1.Tvmult = [& usr_matrix_1, & usr_matrix_2](BlockVector<double> &v, const BlockVector<double> &u) {
    std::vector<types::global_dof_index> sizes_per_block(u.n_blocks());
    for (unsigned int i = 0; i < u.n_blocks(); ++i)
      sizes_per_block[i] = u.block(i).size();
    v.reinit(sizes_per_block);

    usr_matrix_1.Tvmult(v.block(0), u.block(0));
    usr_matrix_2.Tvmult(v.block(1), u.block(1));
  };

  usr_linearOperator_1.Tvmult_add = [usr_linearOperator_1](BlockVector<double> &v, const BlockVector<double> &u) {
    std::vector<types::global_dof_index> sizes_per_block(u.n_blocks());
    for (unsigned int i = 0; i < u.n_blocks(); ++i)
      sizes_per_block[i] = u.block(i).size();

    BlockVector<double> temp_v(sizes_per_block);
    usr_linearOperator_1.Tvmult(temp_v, u);
    v += temp_v;
  };

  usr_linearOperator_1.reinit_range_vector = [& usr_matrix_1, & usr_matrix_2](BlockVector<double> &v, bool omit_zeroing_entries) {
    v.reinit(2);
    (v.block(0)).reinit(usr_matrix_1.m(), omit_zeroing_entries);
    (v.block(1)).reinit(usr_matrix_2.m(), omit_zeroing_entries);
  };

  usr_linearOperator_1.reinit_domain_vector = [& usr_matrix_1, & usr_matrix_2](BlockVector<double> &v, bool omit_zeroing_entries) {
    v.reinit(2);
    (v.block(0)).reinit(usr_matrix_1.n(), omit_zeroing_entries);
    (v.block(1)).reinit(usr_matrix_2.n(), omit_zeroing_entries);
  };

Best,

Tao