]> https://gitweb.dealii.org/ - dealii.git/commitdiff
upper triangular block matrix
authorESeNonFossiIo <esenonfossiio@gmail.com>
Wed, 20 May 2015 13:36:44 +0000 (15:36 +0200)
committerESeNonFossiIo <esenonfossiio@gmail.com>
Sat, 23 May 2015 10:16:24 +0000 (12:16 +0200)
upper triangular block matrix tests

comments

crossed tests

small typos

include/deal.II/lac/linear_operator.h
tests/lac/linear_operator_12.cc [new file with mode: 0644]
tests/lac/linear_operator_12.with_cxx11=on.output [new file with mode: 0644]
tests/lac/linear_operator_13.cc [new file with mode: 0644]
tests/lac/linear_operator_13.with_cxx11=on.output [new file with mode: 0644]
tests/lac/linear_operator_15.cc

index 0f39f75138e7b2e39c351ee53c3c3dd2432422ea..c46a0a44e53003cbcd10af8e44fec0281ca01ca7 100644 (file)
@@ -2229,6 +2229,106 @@ lower_triangular_operator(BlockMatrix &a)
 }
 
 
+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
diff --git a/tests/lac/linear_operator_12.cc b/tests/lac/linear_operator_12.cc
new file mode 100644 (file)
index 0000000..f76d503
--- /dev/null
@@ -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 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);
+  }
+}
diff --git a/tests/lac/linear_operator_12.with_cxx11=on.output b/tests/lac/linear_operator_12.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..822b437
--- /dev/null
@@ -0,0 +1,21 @@
+
+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.    
diff --git a/tests/lac/linear_operator_13.cc b/tests/lac/linear_operator_13.cc
new file mode 100644 (file)
index 0000000..199bd70
--- /dev/null
@@ -0,0 +1,135 @@
+// ---------------------------------------------------------------------
+//
+// 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);
+  }
+}
diff --git a/tests/lac/linear_operator_13.with_cxx11=on.output b/tests/lac/linear_operator_13.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..476923c
--- /dev/null
@@ -0,0 +1,33 @@
+
+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    
index 233b7cbdd7aa27e722fe2283cbbbfe260bdab4cd..18491944f40d92f2eb822def1ccc71a38732172d 100644 (file)
@@ -49,9 +49,9 @@ int main()
     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)
@@ -59,7 +59,7 @@ int main()
 
     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);
@@ -74,7 +74,7 @@ int main()
       v.block(i)[0] = 1;
     }
 
-    triangulat_block_op.vmult(v, u);
+    triangular_block_op.vmult(v, u);
     deallog << " -- vmult --" << std::endl;
     PRINTME("v", v);
 
@@ -82,7 +82,7 @@ int main()
     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);
 
@@ -90,7 +90,7 @@ int main()
     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);
 
@@ -98,7 +98,7 @@ int main()
     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);
   }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.