upper and lower triangular operator for array of array of linear operator

diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h
index d17e290f33ba2307ea0aa3d9854daa71e1f1c165..9118f12ce47456da4ba76ef506279116346db1b8 100644 (file)
--- a/include/deal.II/lac/block_linear_operator.h
+++ b/include/deal.II/lac/block_linear_operator.h
@@ -399,7 +399,7 @@ lower_triangular_operator(const BlockMatrix &block_matrix)
     Assert(u.n_blocks() == block_matrix.n_block_cols(),
            ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols()));
 
-    v.block(0) *= 0;
+    v.block(0) = 0;
     for (unsigned int i = 1; i < block_matrix.n_block_rows(); ++i)
       {
         block_matrix.block(i,0).vmult(v.block(i), u.block(0));
@@ -432,7 +432,7 @@ lower_triangular_operator(const BlockMatrix &block_matrix)
 
     for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
       {
-        v.block(i) *= 0;
+        v.block(i) = 0;
         for (unsigned int j = i + 1; j < block_matrix.n_block_rows(); ++j)
           block_matrix.block(j,i).Tvmult_add(v.block(i), u.block(j));
       }
@@ -454,6 +454,90 @@ lower_triangular_operator(const BlockMatrix &block_matrix)
 }
 
 
+/**
+ * @relates LinearOperator
+ *
+ * This function is a specification of the above function that
+ * allows to work with std::array of std::array of LinearOperator
+ *
+ * @ingroup LAOperators
+ */
+
+template <unsigned int n,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+LinearOperator<Range, Domain>
+lower_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &);
+
+
+template <unsigned int n, typename Range, typename Domain>
+LinearOperator<Range, Domain>
+lower_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &block_matrix)
+{
+  LinearOperator<Range, Domain> return_op;
+
+  return_op.reinit_range_vector = [block_matrix](Range &v, bool fast)
+  {
+    v.reinit(n);
+
+    for (unsigned int i = 0; i < n; ++i)
+      block_matrix[i][0].reinit_range_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.reinit_domain_vector = [block_matrix](Domain &v, bool fast)
+  {
+    v.reinit(n);
+
+    for (unsigned int i = 0; i < n; ++i)
+      block_matrix[0][i].reinit_domain_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.vmult = [block_matrix](Range &v, const Domain &u)
+  {
+    v.block(0) = 0;
+    for (unsigned int i = 1; i < n; ++i)
+      {
+        block_matrix[i][0].vmult(v.block(i), u.block(0));
+        for (unsigned int j = 1; j < i; ++j)
+          block_matrix[i][j].vmult_add(v.block(i), u.block(j));
+      }
+  };
+
+  return_op.vmult_add = [block_matrix](Range &v, const Domain &u)
+  {
+    for (unsigned int i = 1; i < n; ++i)
+      {
+        block_matrix[i][0].vmult_add(v.block(i), u.block(0));
+        for (unsigned int j = 1; j < i; ++j)
+          block_matrix[i][j].vmult_add(v.block(i), u.block(j));
+      }
+  };
+
+  return_op.Tvmult = [block_matrix](Domain &v, const Range &u)
+  {
+    for (unsigned int i = 0; i < n; ++i)
+      {
+        v.block(i) = 0;
+        for (unsigned int j = i+1; j < n; ++j)
+          block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
+      }
+  };
+
+  return_op.Tvmult_add = [block_matrix](Domain &v, const Range &u)
+  {
+    for (unsigned int i = 0; i < n; ++i)
+      for (unsigned int j = i+1; j < n; ++j)
+        block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
+  };
+
+  return return_op;
+}
+
+
 /**
  * @relates LinearOperator
  *
@@ -561,7 +645,7 @@ upper_triangular_operator(const BlockMatrix &block_matrix)
 
     for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
       {
-        v.block(i) *= 0;
+        v.block(i) = 0;
         for (unsigned int j = 0; j < i; ++j)
           block_matrix.block(j,i).vmult_add(v.block(i), u.block(j));
       }
@@ -583,6 +667,87 @@ upper_triangular_operator(const BlockMatrix &block_matrix)
 }
 
 
+/**
+ * @relates LinearOperator
+ *
+ * This function is a specification of the above function that
+ * allows to work with std::array of std::array of LinearOperator
+ *
+ * @ingroup LAOperators
+ */
+
+template <unsigned int n,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+LinearOperator<Range, Domain>
+upper_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &);
+
+
+template <unsigned int n, typename Range, typename Domain>
+LinearOperator<Range, Domain>
+upper_triangular_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, n> &block_matrix)
+{
+  LinearOperator<Range, Domain> return_op;
+
+  return_op.reinit_range_vector = [block_matrix](Range &v, bool fast)
+  {
+    v.reinit(n);
+
+    for (unsigned int i = 0; i < n; ++i)
+      block_matrix[i][0].reinit_range_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.reinit_domain_vector = [block_matrix](Domain &v, bool fast)
+  {
+    v.reinit(n);
+
+    for (unsigned int i = 0; i < n; ++i)
+      block_matrix[0][i].reinit_domain_vector(v.block(i), fast);
+
+    v.collect_sizes();
+  };
+
+  return_op.vmult = [block_matrix](Range &v, const Domain &u)
+  {
+    for (unsigned int i = 0; i < n - 1; ++i)
+      {
+        block_matrix[i][n - 1].vmult(v.block(i), u.block(n - 1));
+        for (unsigned int j = n - 2; j > i; --j)
+          block_matrix[i][j].vmult_add(v.block(i), u.block(j));
+      }
+    v.block(n - 1) = 0;
+  };
+
+  return_op.vmult_add = [block_matrix](Range &v, const Domain &u)
+  {
+    for (unsigned int i = 0; i < n; ++i)
+      for (unsigned int j = i + 1; j < n; ++j)
+        block_matrix[i][j].vmult_add(v.block(i), u.block(j));
+  };
+
+  return_op.Tvmult = [block_matrix](Domain &v, const Range &u)
+  {
+    for (unsigned int i = 0; i < n; ++i)
+      {
+        v.block(i) = 0;
+        for (unsigned int j = 0; j < i; ++j)
+          block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
+      }
+  };
+
+  return_op.Tvmult_add = [block_matrix](Domain &v, const Range &u)
+  {
+    for (unsigned int i = 0; i < n; ++i)
+      for (unsigned int j = 0; j < i; ++j)
+        block_matrix[j][i].Tvmult_add(v.block(i), u.block(j));
+  };
+
+  return return_op;
+}
+
+
 /**
  * @relates LinearOperator
  *

diff --git a/tests/lac/block_linear_operator_07.cc b/tests/lac/block_linear_operator_07.cc
index d902ad86d5b31e1602b71190978e083fc7c42a83..d0c778628f18ef7f248611d02c29e4cca8860f59 100644 (file)
--- a/tests/lac/block_linear_operator_07.cc
+++ b/tests/lac/block_linear_operator_07.cc
@@ -90,16 +90,6 @@ int main()
 
 
     // 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);