lower triangular block matrix

lower triangular matrix tests

diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h
index d02f153fe81c52b91a72cd2dc6e7e742f98d2cd6..0f39f75138e7b2e39c351ee53c3c3dd2432422ea 100644 (file)
--- a/include/deal.II/lac/linear_operator.h
+++ b/include/deal.II/lac/linear_operator.h
@@ -2126,6 +2126,109 @@ operator*(const PackagedOperation<Range> &comp,
   return return_comp;
 }
 
+
+
+template <typename Range = BlockVector<double>,
+          typename Domain = Range,
+          typename BlockMatrix>
+LinearOperator<Range, Domain>
+lower_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()));
+
+    v.block(0) *= 0;
+    for (unsigned int i = 1; i < a.n_block_rows(); ++i)
+      {
+        a.block(i,0).vmult(v.block(i), u.block(0));
+        for (unsigned int j = 1; j < i; ++j)
+          a.block(i,j).vmult_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 = 1; i < a.n_block_rows(); ++i)
+      {
+        a.block(i,0).vmult_add(v.block(i), u.block(0));
+        for (unsigned int j = 1; j < i; ++j)
+          a.block(i,j).vmult_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 = i + 1; j < a.n_block_rows(); ++j)
+          a.block(j,i).Tvmult_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 = i + 1; j < a.n_block_rows(); ++j)
+        a.block(j,i).Tvmult_add(v.block(i), u.block(j));
+  };
+
+  return return_op;
+}
+
+
 //@}
 
 DEAL_II_NAMESPACE_CLOSE

diff --git a/tests/lac/linear_operator_11.with_cxx11=on.output b/tests/lac/linear_operator_11.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..1507567
--- /dev/null
+++ b/tests/lac/linear_operator_11.with_cxx11=on.output
@@ -0,0 +1,21 @@
+
+DEAL:: -- vmult --
+DEAL::Block vector: v:
+DEAL::[block 0 ]  0.0    
+DEAL::[block 1 ]  3.0    
+DEAL::[block 2 ]  10.    
+DEAL:: -- vmult_add --
+DEAL::Block vector: v:
+DEAL::[block 0 ]  1.0    
+DEAL::[block 1 ]  4.0    
+DEAL::[block 2 ]  11.    
+DEAL:: -- Tvmult --
+DEAL::Block vector: u:
+DEAL::[block 0 ]  18.    
+DEAL::[block 1 ]  9.0    
+DEAL::[block 2 ]  0.0    
+DEAL:: -- Tvmult_add --
+DEAL::Block vector: u:
+DEAL::[block 0 ]  19.    
+DEAL::[block 1 ]  10.    
+DEAL::[block 2 ]  1.0    

diff --git a/tests/lac/linear_operator_15.cc b/tests/lac/linear_operator_15.cc
new file mode 100644 (file)
index 0000000..233b7cb
--- /dev/null
+++ b/tests/lac/linear_operator_15.cc
@@ -0,0 +1,105 @@
+// ---------------------------------------------------------------------
+//
+// 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 1 0 |
+    //  | 3 2 1 |
+    //  | 4 3 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 + i - j);
+
+    auto op_a = linear_operator<BlockVector<double>>(a);
+
+    auto triangulat_block_op = lower_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;
+    }
+
+    triangulat_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;
+
+    triangulat_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;
+
+    triangulat_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;
+
+    triangulat_block_op.Tvmult_add(u, v);
+    deallog << " -- Tvmult_add --" << std::endl;
+    PRINTME("u", u);
+  }
+}