From 7d1a007091783cfae4310a9bbd1a778acdbf2b4f Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Tue, 16 Jun 2015 00:05:50 +0200 Subject: [PATCH] test for gauss elimination argument --- tests/lac/block_linear_operator_09.cc | 215 ++++++++++++++++++ ...ck_linear_operator_09.with_cxx11=on.output | 79 +++++++ 2 files changed, 294 insertions(+) create mode 100644 tests/lac/block_linear_operator_09.cc create mode 100644 tests/lac/block_linear_operator_09.with_cxx11=on.output diff --git a/tests/lac/block_linear_operator_09.cc b/tests/lac/block_linear_operator_09.cc new file mode 100644 index 0000000000..49c726f56f --- /dev/null +++ b/tests/lac/block_linear_operator_09.cc @@ -0,0 +1,215 @@ +// --------------------------------------------------------------------- +// +// 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 gauss_upper_block_inverse and gauss_lower_block_inverse: + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include + +#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(12); + + // 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 a (sparsity_pattern); + for (unsigned int i = 0; i < 3; ++i) + { + a.block(i,i).set(0, 0, i+i +1); + for (unsigned int j = 0; j < i; ++j) + a.block(i,j).set(0, 0, 10); + } + + BlockSparseMatrix d(sparsity_pattern); + for (unsigned int i = 0; i < 3; ++i) + d.block(i,i).set(0,0, 1.0 / (i+i +1) ); + + auto op_a = linear_operator< BlockVector >(a); + auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); + auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); + auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); + // auto diagonal_inv = block_diagonal_operator(diag); + auto inverse_op_a = gauss_lower_block_inverse<3, BlockVector, BlockVector, BlockSparseMatrix >(a, {{d00, d11, d22}}); + // + auto identity = inverse_op_a * op_a; + + BlockVector u; + BlockVector v; + + deallog << " -- Matrix -- " << std::endl; + op_a.reinit_domain_vector(u, false); + op_a.reinit_range_vector(v, false); + for(unsigned int j = 0; j<3; ++j) + { + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0;; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1; + + op_a.vmult(v, u); + + PRINTME("v", v); + } + + deallog << " -- Inverse -- " << std::endl; + inverse_op_a.reinit_domain_vector(u, false); + inverse_op_a.reinit_range_vector(v, true); + for(unsigned int j = 0; j<3; ++j) + { + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1;; + + inverse_op_a.vmult(v, u); + + PRINTME("v", v); + } + + deallog << " -- Identity -- " << std::endl; + identity.reinit_domain_vector(u, false); + identity.reinit_range_vector(v, false); + for(unsigned int j = 0; j<3; ++j) + { + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0;; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1;; + + identity.vmult(v, u); + + PRINTME("v", v); + } + } + + + { + 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 a (sparsity_pattern); + for (unsigned int i = 0; i < 3; ++i) + { + a.block(i,i).set(0, 0, i+i +1); + for (unsigned int j = i+1; j < 3; ++j) + a.block(i,j).set(0, 0, 10); + } + + BlockSparseMatrix d(sparsity_pattern); + for (unsigned int i = 0; i < 3; ++i) + d.block(i,i).set(0,0, 1.0 / (i+i +1) ); + + auto op_a = linear_operator< BlockVector >(a); + auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); + auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); + auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); + // auto diagonal_inv = block_diagonal_operator(diag); + auto inverse_op_a = gauss_upper_block_inverse<3, BlockVector, BlockVector, BlockSparseMatrix >(a, {{d00, d11, d22}}); + // + auto identity = inverse_op_a * op_a; + + BlockVector u; + BlockVector v; + + deallog << " -- Matrix -- " << std::endl; + op_a.reinit_domain_vector(u, false); + op_a.reinit_range_vector(v, false); + for(unsigned int j = 0; j<3; ++j) + { + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0;; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1; + + op_a.vmult(v, u); + + PRINTME("v", v); + } + + deallog << " -- Inverse -- " << std::endl; + inverse_op_a.reinit_domain_vector(u, false); + inverse_op_a.reinit_range_vector(v, true); + for(unsigned int j = 0; j<3; ++j) + { + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1;; + + inverse_op_a.vmult(v, u); + + PRINTME("v", v); + } + + deallog << " -- Identity -- " << std::endl; + identity.reinit_domain_vector(u, false); + identity.reinit_range_vector(v, false); + for(unsigned int j = 0; j<3; ++j) + { + for(unsigned int i = 0; i<3; ++i) + { + u.block(i)[0] = 0;; + v.block(i)[0] = 0; + } + u.block(j)[0] = 1;; + + identity.vmult(v, u); + + PRINTME("v", v); + } + } +} diff --git a/tests/lac/block_linear_operator_09.with_cxx11=on.output b/tests/lac/block_linear_operator_09.with_cxx11=on.output new file mode 100644 index 0000000000..1e80e98cd0 --- /dev/null +++ b/tests/lac/block_linear_operator_09.with_cxx11=on.output @@ -0,0 +1,79 @@ + +DEAL:: -- Matrix -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] 10.0000000000 +DEAL::[block 2 ] 10.0000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 3.00000000000 +DEAL::[block 2 ] 10.0000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 5.00000000000 +DEAL:: -- Inverse -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] -3.33333333333 +DEAL::[block 2 ] 4.66666666667 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.333333333333 +DEAL::[block 2 ] -0.666666666667 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 0.200000000000 +DEAL:: -- Identity -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 1.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 1.00000000000 +DEAL:: -- Matrix -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 10.0000000000 +DEAL::[block 1 ] 3.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 10.0000000000 +DEAL::[block 1 ] 10.0000000000 +DEAL::[block 2 ] 5.00000000000 +DEAL:: -- Inverse -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] -3.33333333333 +DEAL::[block 1 ] 0.333333333333 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 4.66666666667 +DEAL::[block 1 ] -0.666666666667 +DEAL::[block 2 ] 0.200000000000 +DEAL:: -- Identity -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 1.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 1.00000000000 -- 2.39.5