From b2a40cee78cefc1d50dc9f709a0b7c9de23b0d38 Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Wed, 3 Jun 2015 10:32:36 +0200 Subject: [PATCH] a test for upper triangular matrix with non trivial block matrices --- tests/lac/block_linear_operator_07.cc | 145 ++++++++++++++++++ ...ck_linear_operator_07.with_cxx11=on.output | 37 +++++ 2 files changed, 182 insertions(+) create mode 100644 tests/lac/block_linear_operator_07.cc create mode 100644 tests/lac/block_linear_operator_07.with_cxx11=on.output diff --git a/tests/lac/block_linear_operator_07.cc b/tests/lac/block_linear_operator_07.cc new file mode 100644 index 0000000000..d902ad86d5 --- /dev/null +++ b/tests/lac/block_linear_operator_07.cc @@ -0,0 +1,145 @@ +// --------------------------------------------------------------------- +// +// 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 upper_triangular_block_op and lower_triangular_operator +// simultaneously: + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include + +#define PRINTME(name, var) \ + deallog << "Block vector: " name << ":" << std::endl; \ + for (unsigned int j = 0; j < 2; ++j) \ + deallog << "[block " << j << " ] " << var.block(j); + +void reinit_vec(BlockVector &u1, BlockVector &u2, BlockVector &v1, BlockVector &v2) +{ + for(unsigned int i = 0; i<2; ++i) + { + for(unsigned int j = 0; j<2; ++j) + { + u1.block(i)[j] = i+j; + u2.block(i)[j] = i+j; + v1.block(i)[j] = 1; + v2.block(i)[j] = 1; + } + } +} + +using namespace dealii; + +int main() +{ + initlog(); + deallog << std::setprecision(2); + + // BlockSparseMatrix: + { + BlockSparsityPattern dsp(2,2); + // set sizes + for (unsigned int i=0; i<2; ++i) + for (unsigned int j=0; j<2; ++j) + dsp.block(i,j).reinit ( 2, 2, 5); + dsp.collect_sizes (); + + for (unsigned int row=0; row<4; ++row) + for (unsigned int i=0; i<4; ++i) + dsp.add (row, i); + dsp.compress (); + + BlockSparseMatrix a (dsp); + + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < 2; ++j) + for (unsigned int h = 0; h < 2; ++h) + for (unsigned int k = 0; k < 2; ++k) + a.block(i,j).set(h, k, j + i + h + k); + + auto op_a = linear_operator>(a); + + auto lower_triangular_block_op = lower_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); + auto upper_triangular_block_op = upper_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); + + + BlockVector u1; + op_a.reinit_domain_vector(u1, false); + BlockVector v1; + op_a.reinit_range_vector(v1, false); + BlockVector u2; + op_a.reinit_domain_vector(u2, false); + BlockVector v2; + op_a.reinit_range_vector(v2, false); + + + // check 1 + // for(unsigned int i = 0; i<2; ++i) + // { + // for(unsigned int j = 0; j<2; ++j) + // { + // u1.block(i)[j] = i+1; + // u2.block(i)[j] = i+1; + // v1.block(i)[j] = 1; + // v2.block(i)[j] = 1; + // } + // } + reinit_vec(u1,u2,v1,v2); + PRINTME("u1", u1); + PRINTME("u2", u2); + PRINTME("v1", v1); + PRINTME("v2", v2); + + + lower_triangular_block_op.vmult(v1, u1); + upper_triangular_block_op.Tvmult(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + + + // check 2 + reinit_vec(u1,u2,v1,v2); + + lower_triangular_block_op.Tvmult(v1, u1); + upper_triangular_block_op.vmult(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + + // check 3 + reinit_vec(u1,u2,v1,v2); + + lower_triangular_block_op.vmult_add(v1, u1); + upper_triangular_block_op.Tvmult_add(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + + + // check 4 + reinit_vec(u1,u2,v1,v2); + + lower_triangular_block_op.Tvmult_add(v1, u1); + upper_triangular_block_op.vmult_add(v2, u2); + + PRINTME("v1", v1); + PRINTME("v2", v2); + } +} diff --git a/tests/lac/block_linear_operator_07.with_cxx11=on.output b/tests/lac/block_linear_operator_07.with_cxx11=on.output new file mode 100644 index 0000000000..9d3ee42e60 --- /dev/null +++ b/tests/lac/block_linear_operator_07.with_cxx11=on.output @@ -0,0 +1,37 @@ + +DEAL::Block vector: u1: +DEAL::[block 0 ] 0.0 1.0 +DEAL::[block 1 ] 1.0 2.0 +DEAL::Block vector: u2: +DEAL::[block 0 ] 0.0 1.0 +DEAL::[block 1 ] 1.0 2.0 +DEAL::Block vector: v1: +DEAL::[block 0 ] 1.0 1.0 +DEAL::[block 1 ] 1.0 1.0 +DEAL::Block vector: v2: +DEAL::[block 0 ] 1.0 1.0 +DEAL::[block 1 ] 1.0 1.0 +DEAL::Block vector: v1: +DEAL::[block 0 ] 0.0 0.0 +DEAL::[block 1 ] 2.0 3.0 +DEAL::Block vector: v2: +DEAL::[block 0 ] 0.0 0.0 +DEAL::[block 1 ] 2.0 3.0 +DEAL::Block vector: v1: +DEAL::[block 0 ] 5.0 8.0 +DEAL::[block 1 ] 0.0 0.0 +DEAL::Block vector: v2: +DEAL::[block 0 ] 5.0 8.0 +DEAL::[block 1 ] 0.0 0.0 +DEAL::Block vector: v1: +DEAL::[block 0 ] 1.0 1.0 +DEAL::[block 1 ] 3.0 4.0 +DEAL::Block vector: v2: +DEAL::[block 0 ] 1.0 1.0 +DEAL::[block 1 ] 3.0 4.0 +DEAL::Block vector: v1: +DEAL::[block 0 ] 6.0 9.0 +DEAL::[block 1 ] 1.0 1.0 +DEAL::Block vector: v2: +DEAL::[block 0 ] 6.0 9.0 +DEAL::[block 1 ] 1.0 1.0 -- 2.39.5