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
On Saturday, April 27, 2024 at 10:33:11 AM UTC-4 Tao Jin wrote:
> 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
>
>
>
--
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/cf5ea8ac-78d4-4d54-90b8-5fd3d972a8b6n%40googlegroups.com.
