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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/4f0dfa54-cec3-4725-ae17-a2ef155e94a3n%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2019 - 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
* This tutorial program was contributed by Martin Kronbichler
*/
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/vector_memory.templates.h>
#include <deal.II/lac/linear_operator_tools.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/block_linear_operator.h>
using namespace dealii;
// define my own operator*() to make the following line work:
// const auto total_op = op_1 * op_2;
// total_op.vmult(v, u);
// where v and u are BlockVector, and op_1 and op_2 are LinearOperator that operates on BlockVector
LinearOperator<BlockVector<double>, BlockVector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload>
operator*(const LinearOperator<BlockVector<double>, BlockVector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload> & first_op,
const LinearOperator<BlockVector<double>, BlockVector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload> &second_op)
{
std::cout << "calling usr-defined operator*() function" << std::endl;
if (first_op.is_null_operator || second_op.is_null_operator)
{
LinearOperator<BlockVector<double>, BlockVector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload> return_op;
return_op.reinit_domain_vector = second_op.reinit_domain_vector;
return_op.reinit_range_vector = first_op.reinit_range_vector;
return null_operator(return_op);
}
else
{
LinearOperator<BlockVector<double>, BlockVector<double>, dealii::internal::LinearOperatorImplementation::EmptyPayload> return_op{
static_cast<const dealii::internal::LinearOperatorImplementation::EmptyPayload &>(first_op) *
static_cast<const dealii::internal::LinearOperatorImplementation::EmptyPayload &>(second_op)};
return_op.reinit_domain_vector = second_op.reinit_domain_vector;
return_op.reinit_range_vector = first_op.reinit_range_vector;
// ensure to have valid computation objects by catching first_op and
// second_op by value
return_op.vmult = [first_op, second_op](BlockVector<double> &v, const BlockVector<double> &u) {
GrowingVectorMemory<BlockVector<double>> vector_memory;
typename VectorMemory<BlockVector<double>>::Pointer i(vector_memory);
(*i).reinit(u.n_blocks());
second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
second_op.vmult(*i, u);
first_op.vmult(v, *i);
};
return_op.vmult_add = [first_op, second_op](BlockVector<double> &v, const BlockVector<double> &u) {
GrowingVectorMemory<BlockVector<double>> vector_memory;
typename VectorMemory<BlockVector<double>>::Pointer i(vector_memory);
(*i).reinit(u.n_blocks());
second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
second_op.vmult(*i, u);
first_op.vmult_add(v, *i);
};
return_op.Tvmult = [first_op, second_op](BlockVector<double> &v, const BlockVector<double> &u) {
GrowingVectorMemory<BlockVector<double>> vector_memory;
typename VectorMemory<BlockVector<double>>::Pointer i(vector_memory);
(*i).reinit(u.n_blocks());
first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
first_op.Tvmult(*i, u);
second_op.Tvmult(v, *i);
};
return_op.Tvmult_add = [first_op, second_op](BlockVector<double> &v, const BlockVector<double> &u) {
GrowingVectorMemory<BlockVector<double>> vector_memory;
typename VectorMemory<BlockVector<double>>::Pointer i(vector_memory);
(*i).reinit(u.n_blocks());
first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
first_op.Tvmult(*i, u);
second_op.Tvmult_add(v, *i);
};
return return_op;
}
}
int
main()
{
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) {
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) {
BlockVector<double> temp_v(u);
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) {
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) {
BlockVector<double> temp_v(u);
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.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.block(0)).reinit(usr_matrix_1.n(), omit_zeroing_entries);
(v.block(1)).reinit(usr_matrix_2.n(), omit_zeroing_entries);
};
LinearOperator<BlockVector<double>, BlockVector<double>, Payload> usr_linearOperator_2;
FullMatrix<double> usr_matrix_3(3, 3);
usr_matrix_3(0,0) = 11;
usr_matrix_3(0,1) = 21;
usr_matrix_3(0,2) = 31;
usr_matrix_3(1,0) = 41;
usr_matrix_3(1,1) = 51;
usr_matrix_3(1,2) = 61;
usr_matrix_3(2,0) = 71;
usr_matrix_3(2,1) = 81;
usr_matrix_3(2,2) = 91;
FullMatrix<double> usr_matrix_4(2, 2);
usr_matrix_4(0,0) = -11;
usr_matrix_4(0,1) = -21;
usr_matrix_4(1,0) = -31;
usr_matrix_4(1,1) = -41;
usr_linearOperator_2.vmult = [& usr_matrix_3, & usr_matrix_4](BlockVector<double> &v, const BlockVector<double> &u) {
usr_matrix_3.vmult(v.block(0), u.block(0));
usr_matrix_4.vmult(v.block(1), u.block(1));
};
usr_linearOperator_2.vmult_add = [&usr_linearOperator_2](BlockVector<double> &v, const BlockVector<double> &u) {
BlockVector<double> temp_v(u);
usr_linearOperator_2.vmult(temp_v, u);
v += temp_v;
};
usr_linearOperator_2.Tvmult = [& usr_matrix_3, & usr_matrix_4](BlockVector<double> &v, const BlockVector<double> &u) {
usr_matrix_3.Tvmult(v.block(0), u.block(0));
usr_matrix_4.Tvmult(v.block(1), u.block(1));
};
usr_linearOperator_2.Tvmult_add = [usr_linearOperator_2](BlockVector<double> &v, const BlockVector<double> &u) {
BlockVector<double> temp_v(u);
usr_linearOperator_2.Tvmult(temp_v, u);
v += temp_v;
};
usr_linearOperator_2.reinit_range_vector = [& usr_matrix_3, & usr_matrix_4](BlockVector<double> &v, bool omit_zeroing_entries) {
(v.block(0)).reinit(usr_matrix_3.m(), omit_zeroing_entries);
(v.block(1)).reinit(usr_matrix_4.m(), omit_zeroing_entries);
};
usr_linearOperator_2.reinit_domain_vector = [& usr_matrix_3, & usr_matrix_4](BlockVector<double> &v, bool omit_zeroing_entries) {
(v.block(0)).reinit(usr_matrix_3.n(), omit_zeroing_entries);
(v.block(1)).reinit(usr_matrix_4.n(), omit_zeroing_entries);
};
std::vector<types::global_dof_index> block_vector_sizes(2);
block_vector_sizes[0] = 3;
block_vector_sizes[1] = 2;
BlockVector<double> u(block_vector_sizes);
u.block(0)[0] = 1;
u.block(0)[1] = 2;
u.block(0)[2] = 3;
u.block(1)[0] = -1;
u.block(1)[1] = -2;
BlockVector<double> v(block_vector_sizes);
auto op_total = usr_linearOperator_1 * usr_linearOperator_2;
op_total.vmult(v, u);
std::cout << "usr_linearOperator_1 * usr_linearOperator_2 =>" << std::endl;
std::cout << "v.block(0)[0] = " << v.block(0)[0] << std::endl;
std::cout << "v.block(0)[1] = " << v.block(0)[1] << std::endl;
std::cout << "v.block(0)[2] = " << v.block(0)[2] << std::endl;
std::cout << "v.block(1)[0] = " << v.block(1)[0] << std::endl;
std::cout << "v.block(1)[1] = " << v.block(1)[1] << std::endl;
op_total = usr_linearOperator_2 * usr_linearOperator_1;
op_total.vmult(v, u);
std::cout << "usr_linearOperator_2 * usr_linearOperator_1 =>" << std::endl;
std::cout << "v.block(0)[0] = " << v.block(0)[0] << std::endl;
std::cout << "v.block(0)[1] = " << v.block(0)[1] << std::endl;
std::cout << "v.block(0)[2] = " << v.block(0)[2] << std::endl;
std::cout << "v.block(1)[0] = " << v.block(1)[0] << std::endl;
std::cout << "v.block(1)[1] = " << v.block(1)[1] << std::endl;
return 0;
}
