a test for upper triangular matrix with non trivial block matrices

diff --git a/tests/lac/block_linear_operator_07.cc b/tests/lac/block_linear_operator_07.cc
new file mode 100644 (file)
index 0000000..d902ad8
--- /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 <deal.II/lac/block_linear_operator.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/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#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<double> &u1, BlockVector<double> &u2, BlockVector<double> &v1, BlockVector<double> &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<double> 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<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<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 (file)
index 0000000..9d3ee42
--- /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