}
+template <typename Range = BlockVector<double>,
+ typename Domain = Range,
+ typename BlockMatrix>
+LinearOperator<Range, Domain>
+upper_triangular_operator(BlockMatrix &a)
+{
+ Assert( a.n_block_rows() == a.n_block_cols(),
+ ExcDimensionMismatch(a.n_block_rows(),a.n_block_cols()) );
+
+ LinearOperator<Range, Domain> return_op;
+
+ return_op.reinit_range_vector = [&a](Range &v, bool fast)
+ {
+ // Reinitialize the block vector to have the number of blocks
+ // equal to the number of row blocks of the matrix a.
+ v.reinit(a.n_block_rows());
+
+ // And reinitialize every individual block with reinit_range_vectors:
+ for (unsigned int i = 0; i < a.n_block_rows(); ++i)
+ linear_operator(a.block(i,0)).reinit_range_vector(v.block(i), fast);
+
+ v.collect_sizes();
+ };
+
+ return_op.reinit_domain_vector = [&a](Domain &v, bool fast)
+ {
+ // Reinitialize the block vector to have the number of blocks
+ // equal to the number of coloumn blocks of the matrix a.
+ v.reinit(a.n_block_cols());
+
+ // And reinitialize every individual block with reinit_domain_vectors:
+ for (unsigned int i = 0; i < a.n_block_cols(); ++i)
+ linear_operator(a.block(0,i)).reinit_domain_vector(v.block(i), fast);
+
+ v.collect_sizes();
+ };
+
+ return_op.vmult = [&a](Range &v, const Domain &u)
+ {
+ Assert( v.n_blocks() == a.n_block_rows(),
+ ExcDimensionMismatch(v.n_blocks(), a.n_block_rows()));
+ Assert( u.n_blocks() == a.n_block_cols(),
+ ExcDimensionMismatch(u.n_blocks(), a.n_block_cols()));
+
+ for (unsigned int i = 0; i < a.n_block_rows(); ++i)
+ {
+ v.block(i) *= 0;
+ for (unsigned int j = i + 1; j < a.n_block_cols(); ++j)
+ a.block(i,j).Tvmult_add(v.block(i), u.block(j));
+ }
+
+ };
+
+ return_op.vmult_add = [&a](Range &v, const Domain &u)
+ {
+ Assert( v.n_blocks() == a.n_block_rows(),
+ ExcDimensionMismatch(v.n_blocks(), a.n_block_rows()));
+ Assert( u.n_blocks() == a.n_block_cols(),
+ ExcDimensionMismatch(u.n_blocks(), a.n_block_cols()));
+
+ for (unsigned int i = 0; i < a.n_block_cols(); ++i)
+ for (unsigned int j = i + 1; j < a.n_block_rows(); ++j)
+ a.block(i,j).Tvmult_add(v.block(i), u.block(j));
+ };
+
+
+
+ return_op.Tvmult = [&a](Domain &v, const Range &u)
+ {
+ Assert( v.n_blocks() == a.n_block_cols(),
+ ExcDimensionMismatch(v.n_blocks(), a.n_block_cols()));
+ Assert( u.n_blocks() == a.n_block_rows(),
+ ExcDimensionMismatch(u.n_blocks(), a.n_block_rows()));
+
+
+
+ for (unsigned int i = 0; i < a.n_block_cols(); ++i)
+ {
+ v.block(i) *= 0;
+ for (unsigned int j = 0; j < i; ++j)
+ a.block(j,i).vmult_add(v.block(i), u.block(j));
+ }
+ };
+
+
+ return_op.Tvmult_add = [&a](Domain &v, const Range &u)
+ {
+ Assert( v.n_blocks() == a.n_block_cols(),
+ ExcDimensionMismatch(v.n_blocks(), a.n_block_cols()));
+ Assert( u.n_blocks() == a.n_block_rows(),
+ ExcDimensionMismatch(u.n_blocks(), a.n_block_rows()));
+
+ for (unsigned int i = 0; i < a.n_block_cols(); ++i)
+ for (unsigned int j = 0; j < i; ++j)
+ a.block(j,i).vmult_add(v.block(i), u.block(j));
+ };
+
+ return return_op;
+}
+
//@}
DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test non symmetric variants:
+
+#include "../tests.h"
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#define PRINTME(name, var) \
+ deallog << "Block vector: " name << ":" << std::endl; \
+ for (unsigned int i = 0; i < var.n_blocks(); ++i) \
+ deallog << "[block " << i << " ] " << var.block(i);
+
+
+using namespace dealii;
+
+int main()
+{
+ initlog();
+ deallog << std::setprecision(2);
+
+ // BlockSparseMatrix:
+ {
+ BlockDynamicSparsityPattern dsp(3, 3);
+ for (unsigned int i = 0; i < 3; ++i)
+ for (unsigned int j = 0; j < 3; ++j)
+ dsp.block(i, j).reinit (1, 1);
+ dsp.collect_sizes ();
+
+ BlockSparsityPattern sparsity_pattern;
+ sparsity_pattern.copy_from(dsp);
+ sparsity_pattern.compress();
+
+ // | 2 3 4 |
+ // | 1 2 3 |
+ // | 0 1 2 |
+ BlockSparseMatrix<double> a (sparsity_pattern);
+ for (unsigned int i = 0; i < 3; ++i)
+ for (unsigned int j = 0; j < 3; ++j)
+ a.block(i,j).set(0, 0, 2 + j - i);
+
+ auto op_a = linear_operator<BlockVector<double>>(a);
+
+ auto triangular_block_op = upper_triangular_operator< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a);
+
+ BlockVector<double> u;
+ op_a.reinit_domain_vector(u, false);
+ BlockVector<double> v;
+ op_a.reinit_range_vector(v, false);
+
+
+ // vmult:
+ for(unsigned int i = 0; i<3; ++i)
+ {
+ u.block(i)[0] = i+1;
+ v.block(i)[0] = 1;
+ }
+
+ triangular_block_op.vmult(v, u);
+ deallog << " -- vmult --" << std::endl;
+ PRINTME("v", v);
+
+ // vmult_add:
+ for(unsigned int i = 0; i<3; ++i)
+ v.block(i)[0] = 1;
+
+ triangular_block_op.vmult_add(v, u);
+ deallog << " -- vmult_add --" << std::endl;
+ PRINTME("v", v);
+
+ // Tvmult
+ for(unsigned int i = 0; i<3; ++i)
+ v.block(i)[0] = i+1;
+
+ triangular_block_op.Tvmult(u, v);
+ deallog << " -- Tvmult --" << std::endl;
+ PRINTME("u", u);
+
+ // Tvmult_add
+ for(unsigned int i = 0; i<3; ++i)
+ u.block(i)[0] = 1;
+
+ triangular_block_op.Tvmult_add(u, v);
+ deallog << " -- Tvmult_add --" << std::endl;
+ PRINTME("u", u);
+ }
+}
--- /dev/null
+
+DEAL:: -- vmult --
+DEAL::Block vector: v:
+DEAL::[block 0 ] 18.
+DEAL::[block 1 ] 9.0
+DEAL::[block 2 ] 0.0
+DEAL:: -- vmult_add --
+DEAL::Block vector: v:
+DEAL::[block 0 ] 19.
+DEAL::[block 1 ] 10.
+DEAL::[block 2 ] 1.0
+DEAL:: -- Tvmult --
+DEAL::Block vector: u:
+DEAL::[block 0 ] 0.0
+DEAL::[block 1 ] 3.0
+DEAL::[block 2 ] 10.
+DEAL:: -- Tvmult_add --
+DEAL::Block vector: u:
+DEAL::[block 0 ] 1.0
+DEAL::[block 1 ] 4.0
+DEAL::[block 2 ] 11.
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test non symmetric variants:
+
+#include "../tests.h"
+
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#define PRINTME(name, var) \
+ deallog << "Block vector: " name << ":" << std::endl; \
+ for (unsigned int i = 0; i < var.n_blocks(); ++i) \
+ deallog << "[block " << i << " ] " << var.block(i);
+
+
+using namespace dealii;
+
+int main()
+{
+ initlog();
+ deallog << std::setprecision(2);
+
+ // BlockSparseMatrix:
+ {
+ BlockDynamicSparsityPattern dsp(3, 3);
+ for (unsigned int i = 0; i < 3; ++i)
+ for (unsigned int j = 0; j < 3; ++j)
+ dsp.block(i, j).reinit (1, 1);
+ dsp.collect_sizes ();
+
+ BlockSparsityPattern sparsity_pattern;
+ sparsity_pattern.copy_from(dsp);
+ sparsity_pattern.compress();
+
+ BlockSparseMatrix<double> a (sparsity_pattern);
+ for (unsigned int i = 0; i < 3; ++i)
+ for (unsigned int j = 0; j < 3; ++j)
+ a.block(i,j).set(0, 0, 2 + j + i);
+
+ auto op_a = linear_operator<BlockVector<double>>(a);
+
+ auto lower_triangular_block_op = lower_triangular_operator< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a);
+ auto upper_triangular_block_op = upper_triangular_operator< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a);
+
+
+ BlockVector<double> u1;
+ op_a.reinit_domain_vector(u1, false);
+ BlockVector<double> v1;
+ op_a.reinit_range_vector(v1, false);
+ BlockVector<double> u2;
+ op_a.reinit_domain_vector(u2, false);
+ BlockVector<double> v2;
+ op_a.reinit_range_vector(v2, false);
+
+
+ // check 1
+ for(unsigned int i = 0; i<3; ++i)
+ {
+ u1.block(i)[0] = i+1;
+ u2.block(i)[0] = i+1;
+ v1.block(i)[0] = 1;
+ v2.block(i)[0] = 1;
+ }
+
+ lower_triangular_block_op.vmult(v1, u1);
+ upper_triangular_block_op.Tvmult(v2, u2);
+
+ PRINTME("v1", v1);
+ PRINTME("v2", v2);
+
+
+ // check 2
+ for(unsigned int i = 0; i<3; ++i)
+ {
+ u1.block(i)[0] = i+1;
+ u2.block(i)[0] = i+1;
+ v1.block(i)[0] = 1;
+ v2.block(i)[0] = 1;
+ }
+
+ lower_triangular_block_op.Tvmult(v1, u1);
+ upper_triangular_block_op.vmult(v2, u2);
+
+ PRINTME("v1", v1);
+ PRINTME("v2", v2);
+
+ // check 3
+ for(unsigned int i = 0; i<3; ++i)
+ {
+ u1.block(i)[0] = i+1;
+ u2.block(i)[0] = i+1;
+ v1.block(i)[0] = 1;
+ v2.block(i)[0] = 1;
+ }
+
+ lower_triangular_block_op.vmult_add(v1, u1);
+ upper_triangular_block_op.Tvmult_add(v2, u2);
+
+ PRINTME("v1", v1);
+ PRINTME("v2", v2);
+
+
+ // check 4
+ for(unsigned int i = 0; i<3; ++i)
+ {
+ u1.block(i)[0] = i+1;
+ u2.block(i)[0] = i+1;
+ v1.block(i)[0] = 1;
+ v2.block(i)[0] = 1;
+ }
+
+ lower_triangular_block_op.Tvmult_add(v1, u1);
+ upper_triangular_block_op.vmult_add(v2, u2);
+
+ PRINTME("v1", v1);
+ PRINTME("v2", v2);
+ }
+}
--- /dev/null
+
+DEAL::Block vector: v1:
+DEAL::[block 0 ] 0.0
+DEAL::[block 1 ] 3.0
+DEAL::[block 2 ] 14.
+DEAL::Block vector: v2:
+DEAL::[block 0 ] 0.0
+DEAL::[block 1 ] 3.0
+DEAL::[block 2 ] 14.
+DEAL::Block vector: v1:
+DEAL::[block 0 ] 18.
+DEAL::[block 1 ] 15.
+DEAL::[block 2 ] 0.0
+DEAL::Block vector: v2:
+DEAL::[block 0 ] 18.
+DEAL::[block 1 ] 15.
+DEAL::[block 2 ] 0.0
+DEAL::Block vector: v1:
+DEAL::[block 0 ] 1.0
+DEAL::[block 1 ] 4.0
+DEAL::[block 2 ] 15.
+DEAL::Block vector: v2:
+DEAL::[block 0 ] 1.0
+DEAL::[block 1 ] 4.0
+DEAL::[block 2 ] 15.
+DEAL::Block vector: v1:
+DEAL::[block 0 ] 19.
+DEAL::[block 1 ] 16.
+DEAL::[block 2 ] 1.0
+DEAL::Block vector: v2:
+DEAL::[block 0 ] 19.
+DEAL::[block 1 ] 16.
+DEAL::[block 2 ] 1.0
sparsity_pattern.copy_from(dsp);
sparsity_pattern.compress();
- // | 2 1 0 |
- // | 3 2 1 |
- // | 4 3 2 |
+ // | 2 1 0 | | 0 0 0 |
+ // | 3 2 1 | -> | 3 0 0 |
+ // | 4 3 2 | | 4 3 0 |
BlockSparseMatrix<double> a (sparsity_pattern);
for (unsigned int i = 0; i < 3; ++i)
for (unsigned int j = 0; j < 3; ++j)
auto op_a = linear_operator<BlockVector<double>>(a);
- auto triangulat_block_op = lower_triangular_operator< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a);
+ auto triangular_block_op = lower_triangular_operator< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a);
BlockVector<double> u;
op_a.reinit_domain_vector(u, false);
v.block(i)[0] = 1;
}
- triangulat_block_op.vmult(v, u);
+ triangular_block_op.vmult(v, u);
deallog << " -- vmult --" << std::endl;
PRINTME("v", v);
for(unsigned int i = 0; i<3; ++i)
v.block(i)[0] = 1;
- triangulat_block_op.vmult_add(v, u);
+ triangular_block_op.vmult_add(v, u);
deallog << " -- vmult_add --" << std::endl;
PRINTME("v", v);
for(unsigned int i = 0; i<3; ++i)
v.block(i)[0] = i+1;
- triangulat_block_op.Tvmult(u, v);
+ triangular_block_op.Tvmult(u, v);
deallog << " -- Tvmult --" << std::endl;
PRINTME("u", u);
for(unsigned int i = 0; i<3; ++i)
u.block(i)[0] = 1;
- triangulat_block_op.Tvmult_add(u, v);
+ triangular_block_op.Tvmult_add(u, v);
deallog << " -- Tvmult_add --" << std::endl;
PRINTME("u", u);
}
projects / dealii.git / commitdiff