documentation

indentation and an assert

layout

a small error

documentation

small bugs + rebase

small change in the comment

notation

a small comment

renaming the function

a -> block_matrix

diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h
index 6cadd62cd240b057bee327f3f8a62b155ba4129e..73b0598c63ba641e63a6fe61aa10b4ea6c9ee7eb 100644 (file)
--- a/include/deal.II/lac/linear_operator.h
+++ b/include/deal.II/lac/linear_operator.h
@@ -2127,234 +2127,286 @@ operator*(const PackagedOperation<Range> &comp,
 }
 
 
-
+/**
+ * @relates LinearOperator
+ *
+ * A function that takes a block matrix @p a and returns its associated
+ * lower triangular matrix operator (diagonal is not included).
+ *
+ * @code
+ *  a00 | a01 | a02                 |     |
+ *  ---------------             ---------------
+ *  a10 | a11 | a12     ->      a10 |     |
+ *  ---------------             ---------------
+ *  a20 | a21 | a22             a20 | a21 |
+ * @endcode
+ *
+ * @ingroup LAOperators
+ */
 template <typename Range = BlockVector<double>,
           typename Domain = Range,
           typename BlockMatrix>
 LinearOperator<Range, Domain>
-lower_triangular_operator(BlockMatrix &a)
+lower_triangular_operator(const BlockMatrix &block_matrix)
 {
-  Assert( a.n_block_rows() == a.n_block_cols(),
-          ExcDimensionMismatch(a.n_block_rows(),a.n_block_cols()) );
+  Assert( block_matrix.n_block_rows() == block_matrix.n_block_cols(),
+          ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) );
 
   LinearOperator<Range, Domain> return_op;
 
-  return_op.reinit_range_vector = [&a](Range &v, bool fast)
+  return_op.reinit_range_vector = [&block_matrix](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());
+    // equal to the number of row blocks of the matrix block_matrix.
+    v.reinit(block_matrix.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);
+    for (unsigned int i = 0; i < block_matrix.n_block_rows(); ++i)
+      linear_operator(block_matrix.block(i,0)).reinit_range_vector(v.block(i), fast);
 
     v.collect_sizes();
   };
 
-  return_op.reinit_domain_vector = [&a](Domain &v, bool fast)
+  return_op.reinit_domain_vector = [&block_matrix](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());
+    // equal to the number of coloumn blocks of the matrix block_matrix.
+    v.reinit(block_matrix.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);
+    for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
+      linear_operator(block_matrix.block(0,i)).reinit_domain_vector(v.block(i), fast);
 
     v.collect_sizes();
   };
 
-  return_op.vmult = [&a](Range &v, const Domain &u)
+  return_op.vmult = [&block_matrix](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()));
+    Assert( v.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows()));
+    Assert( u.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols()));
 
     v.block(0) *= 0;
-    for (unsigned int i = 1; i < a.n_block_rows(); ++i)
+    for (unsigned int i = 1; i < block_matrix.n_block_rows(); ++i)
       {
-        a.block(i,0).vmult(v.block(i), u.block(0));
+        block_matrix.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));
+          block_matrix.block(i,j).vmult_add(v.block(i), u.block(j));
       }
   };
 
-  return_op.vmult_add = [&a](Range &v, const Domain &u)
+  return_op.vmult_add = [&block_matrix](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()));
+    Assert( v.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows()));
+    Assert( u.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols()));
 
-    for (unsigned int i = 1; i < a.n_block_rows(); ++i)
+    for (unsigned int i = 1; i < block_matrix.n_block_rows(); ++i)
       {
-        a.block(i,0).vmult_add(v.block(i), u.block(0));
+        block_matrix.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));
+          block_matrix.block(i,j).vmult_add(v.block(i), u.block(j));
       }
   };
 
-
-  return_op.Tvmult = [&a](Domain &v, const Range &u)
+  return_op.Tvmult = [&block_matrix](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()));
-
+    Assert( v.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols()));
+    Assert( u.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows()));
 
-    for (unsigned int i = 0; i < a.n_block_cols(); ++i)
+    for (unsigned int i = 0; i < block_matrix.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));
+        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));
       }
   };
 
-  return_op.Tvmult_add = [&a](Domain &v, const Range &u)
+  return_op.Tvmult_add = [&block_matrix](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()));
+    Assert( v.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols()));
+    Assert( u.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.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));
+    for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
+      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));
   };
 
   return return_op;
 }
 
 
+/**
+ * @relates LinearOperator
+ *
+ * A function that takes a block matrix @p a and returns its associated
+ * upper triangular matrix operator (diagonal is not included).
+ *
+ * @code
+ *  a00 | a01 | a02                 | a01 | a02
+ *  ---------------             ---------------
+ *  a10 | a11 | a12     ->          |     | a12
+ *  ---------------             ---------------
+ *  a20 | a21 | a22                 |     |
+ * @endcode
+ *
+ * @ingroup LAOperators
+ */
 template <typename Range = BlockVector<double>,
           typename Domain = Range,
           typename BlockMatrix>
 LinearOperator<Range, Domain>
-upper_triangular_operator(BlockMatrix &a)
+upper_triangular_operator(const BlockMatrix &block_matrix)
 {
-  Assert( a.n_block_rows() == a.n_block_cols(),
-          ExcDimensionMismatch(a.n_block_rows(),a.n_block_cols()) );
+  Assert( block_matrix.n_block_rows() == block_matrix.n_block_cols(),
+          ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) );
 
   LinearOperator<Range, Domain> return_op;
 
-  return_op.reinit_range_vector = [&a](Range &v, bool fast)
+  return_op.reinit_range_vector = [&block_matrix](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());
+    // equal to the number of row blocks of the matrix block_matrix.
+    v.reinit(block_matrix.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);
+    for (unsigned int i = 0; i < block_matrix.n_block_rows(); ++i)
+      linear_operator(block_matrix.block(i,0)).reinit_range_vector(v.block(i), fast);
 
     v.collect_sizes();
   };
 
-  return_op.reinit_domain_vector = [&a](Domain &v, bool fast)
+  return_op.reinit_domain_vector = [&block_matrix](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());
+    // equal to the number of coloumn blocks of the matrix block_matrix.
+    v.reinit(block_matrix.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);
+    for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
+      linear_operator(block_matrix.block(0,i)).reinit_domain_vector(v.block(i), fast);
 
     v.collect_sizes();
   };
 
-  return_op.vmult = [&a](Range &v, const Domain &u)
+  return_op.vmult = [&block_matrix](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()));
+    Assert( v.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows()));
+    Assert( u.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols()));
 
-    for (unsigned int i = 0; i < a.n_block_rows(); ++i)
+    for (unsigned int i = 0; i < block_matrix.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));
+        for (unsigned int j = i + 1; j < block_matrix.n_block_cols(); ++j)
+          block_matrix.block(i,j).Tvmult_add(v.block(i), u.block(j));
       }
 
   };
 
-  return_op.vmult_add = [&a](Range &v, const Domain &u)
+  return_op.vmult_add = [&block_matrix](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()));
+    Assert( v.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows()));
+    Assert( u.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.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));
+    for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
+      for (unsigned int j = i + 1; j < block_matrix.n_block_rows(); ++j)
+        block_matrix.block(i,j).Tvmult_add(v.block(i), u.block(j));
   };
 
+  return_op.Tvmult = [&block_matrix](Domain &v, const Range &u)
+  {
+    Assert( v.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols()));
+    Assert( u.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows()));
 
+    for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
+      {
+        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));
+      }
+  };
 
-    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)
+  return_op.Tvmult_add = [&block_matrix](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()));
+    Assert( v.n_blocks() == block_matrix.n_block_cols(),
+            ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols()));
+    Assert( u.n_blocks() == block_matrix.n_block_rows(),
+            ExcDimensionMismatch(u.n_blocks(), block_matrix.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));
+    for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i)
+      for (unsigned int j = 0; j < i; ++j)
+        block_matrix.block(j,i).vmult_add(v.block(i), u.block(j));
   };
 
   return return_op;
 }
 
 
+/**
+ * @relates LinearOperator
+ *
+ * Let M be a block matrix of the form Id + T made of nxn blocks
+ * and where T is a lower / upper triangular (without diagonal).
+ * Then, its inverse is of the form:
+ * @code
+ *  Id + sum_{i=1}^{n-1} (-1)^i T^i
+ * @endcode
+ * This formula can be used to invert all triangular matrices
+ * (diagonal included of course).
+ *
+ * This function takes a block matrix @p block_matrix (possibly full) and
+ * a linear block operator  @p inverse_diagonal made of the inverse
+ * of the diagonal blocks inverses.
+ * The output is the inverse of the matrix in the case of a triangular matrix and
+ * as inverse_diagonal its diagonal blocks inverses.
+ * Otherwise, the result is a preconditioner.
+ *
+ * The parameter @p lower is a bool that allows to specify if we want
+ * to use lower triangular part of @p block_matrix (true, this is the default value) or
+ * to use upper triangular part of @p block_matrix (false).
+ * @ingroup LAOperators
+ */
 template <typename Range = BlockVector<double>,
           typename Domain = Range,
           typename BlockMatrix>
 LinearOperator<Range, Domain>
-inverse_operator( BlockMatrix &a,
-                  const LinearOperator<Range, Domain> &inverse_diagonal,
-                  bool lower = true)
+block_triangular_inverse( const BlockMatrix &block_matrix,
+                          const LinearOperator<Range, Domain> &inverse_diagonal,
+                          bool lower = true)
 {
+  Assert( block_matrix.n_block_rows() == block_matrix.n_block_cols(),
+          ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) );
 
   LinearOperator<Domain, Range> op_a;
 
-  if(lower)
-  {
-    op_a = lower_triangular_operator(a);
-  }
+  if (lower)
+    {
+      op_a = lower_triangular_operator(block_matrix);
+    }
   else
-  {
-    op_a = upper_triangular_operator(a);
-  }
+    {
+      op_a = upper_triangular_operator(block_matrix);
+    }
 
   auto id     = identity_operator(op_a.reinit_range_vector);
   auto result = identity_operator(op_a.reinit_range_vector);
 
-
-  for(unsigned int i = 0; i < a.n_block_cols() - 1; ++i)
+  // Notice that the following formula is recursive.
+  //  we are evaluating:
+  //  Id - T + T^2 - T^3 ... (- T)^block_matrix.n_block_cols()
+  for (unsigned int i = 0; i < block_matrix.n_block_cols() - 1; ++i)
     result = id - inverse_diagonal * op_a * result;
 
   return result * inverse_diagonal;

diff --git a/tests/lac/linear_operator_15.cc b/tests/lac/linear_operator_11.cc
similarity index 88%
rename from tests/lac/linear_operator_15.cc
rename to tests/lac/linear_operator_11.cc
index 18491944f40d92f2eb822def1ccc71a38732172d..02c674bd7b349190a27b3556c2da0fca0037241a 100644 (file)
--- a/tests/lac/linear_operator_15.cc
+++ b/tests/lac/linear_operator_11.cc
@@ -13,7 +13,7 @@
 //
 // ---------------------------------------------------------------------
 
-// Test non symmetric variants:
+// Test lower_triangular_operator:
 
 #include "../tests.h"
 
@@ -24,7 +24,8 @@
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/vector.h>
 
-#define PRINTME(name, var) \
+#define PRINTME(title, name, var) \
+  deallog << title << std::endl; \
   deallog << "Block vector: " name << ":" << std::endl; \
   for (unsigned int i = 0; i < var.n_blocks(); ++i) \
     deallog << "[block " << i << " ]  " << var.block(i);
@@ -58,7 +59,6 @@ int main()
         a.block(i,j).set(0, 0, 2 + i - j);
 
     auto op_a = linear_operator<BlockVector<double>>(a);
-
     auto triangular_block_op = lower_triangular_operator< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a);
 
     BlockVector<double> u;
@@ -75,31 +75,27 @@ int main()
     }
 
     triangular_block_op.vmult(v, u);
-    deallog << " -- vmult --" << std::endl;
-    PRINTME("v", v);
+    PRINTME(" -- vmult --", "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);
+    PRINTME(" -- vmult_add --", "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);
+    PRINTME(" -- Tvmult --", "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);
+    PRINTME(" -- Tvmult_add --", "u", u);
   }
 }

diff --git a/tests/lac/linear_operator_12.cc b/tests/lac/linear_operator_12.cc
index f76d503b95ef0bdd238da94365fe6c97b70b96db..1b6601d624c6f38766fbb1933a3c5123564d86c5 100644 (file)
--- a/tests/lac/linear_operator_12.cc
+++ b/tests/lac/linear_operator_12.cc
@@ -13,7 +13,7 @@
 //
 // ---------------------------------------------------------------------
 
-// Test non symmetric variants:
+// Test upper_triangular_operator:
 
 #include "../tests.h"
 
@@ -24,12 +24,12 @@
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/vector.h>
 
-#define PRINTME(name, var) \
+#define PRINTME(title, name, var) \
+  deallog << title << std::endl; \
   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()
@@ -75,31 +75,27 @@ int main()
     }
 
     triangular_block_op.vmult(v, u);
-    deallog << " -- vmult --" << std::endl;
-    PRINTME("v", v);
+    PRINTME(" -- vmult --","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);
+    PRINTME(" -- vmult_add --", "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);
+    PRINTME(" -- Tvmult --", "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);
+    PRINTME(" -- Tvmult_add --", "u", u);
   }
 }

diff --git a/tests/lac/linear_operator_13.cc b/tests/lac/linear_operator_13.cc
index 199bd702d3231cd0f52cd6c871295bc2f280478e..61488e2f89d0e04865774d70dff265721824a6dc 100644 (file)
--- a/tests/lac/linear_operator_13.cc
+++ b/tests/lac/linear_operator_13.cc
@@ -13,7 +13,8 @@
 //
 // ---------------------------------------------------------------------
 
-// Test non symmetric variants:
+// Test upper_triangular_block_op and lower_triangular_operator
+// simultaneously:
 
 #include "../tests.h"
 

diff --git a/tests/lac/linear_operator_14.cc b/tests/lac/linear_operator_14.cc
index a722ad1d7a448812a42975d96ccb3e552e8cce18..7a49f4574bf7332ad672cb844c1207feb854f8f2 100644 (file)
--- a/tests/lac/linear_operator_14.cc
+++ b/tests/lac/linear_operator_14.cc
@@ -13,7 +13,7 @@
 //
 // ---------------------------------------------------------------------
 
-// Test non symmetric variants:
+// Test preconditioner_from_diagonal_inverse:
 
 #include "../tests.h"
 
@@ -60,13 +60,11 @@ int main()
     BlockSparseMatrix<double> d(sparsity_pattern);
     for (unsigned int i = 0; i < 3; ++i)
         d.set(i,i, 1.0 / (i+i +1) );
-    // std::array<decltype(a.block(0,0)), 3> diag {{a.block(0,0), a.block(1,1),a.block(2,2)}};
-
 
     auto op_a         = linear_operator< BlockVector<double> >(a);
     auto diagonal_inv = linear_operator< BlockVector<double> >(d);
     // auto diagonal_inv = block_diagonal_operator(diag);
-    auto inverse_op_a = inverse_operator< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a, diagonal_inv);
+    auto inverse_op_a = block_triangular_inverse< BlockVector<double>, BlockVector<double>, BlockSparseMatrix<double> >(a, diagonal_inv);
 
     auto identity = inverse_op_a * op_a;
 
@@ -78,7 +76,6 @@ int main()
     op_a.reinit_range_vector(v, false);
     for(unsigned int j = 0; j<3; ++j)
     {
-      // check 1
       for(unsigned int i = 0; i<3; ++i)
       {
         u.block(i)[0] = 0;;
@@ -86,9 +83,7 @@ int main()
       }
       u.block(j)[0] = 1;;
 
-
       op_a.vmult(v, u);
-      // identity.Tvmult(v2, u2);
 
       PRINTME("v", v);
     }
@@ -98,17 +93,14 @@ int main()
     inverse_op_a.reinit_range_vector(v, false);
     for(unsigned int j = 0; j<3; ++j)
     {
-      // check 1
       for(unsigned int i = 0; i<3; ++i)
       {
-        u.block(i)[0] = 0;;
+        u.block(i)[0] = 0;
         v.block(i)[0] = 0;
       }
       u.block(j)[0] = 1;;
 
-
       inverse_op_a.vmult(v, u);
-      // identity.Tvmult(v2, u2);
 
       PRINTME("v", v);
     }
@@ -118,7 +110,6 @@ int main()
     identity.reinit_range_vector(v, false);
     for(unsigned int j = 0; j<3; ++j)
     {
-      // check 1
       for(unsigned int i = 0; i<3; ++i)
       {
         u.block(i)[0] = 0;;
@@ -126,21 +117,9 @@ int main()
       }
       u.block(j)[0] = 1;;
 
-
       identity.vmult(v, u);
-      // identity.Tvmult(v2, u2);
 
       PRINTME("v", v);
     }
-
-    // op_a.vmult(v, u);
-    // // identity.Tvmult(v2, u2);
-    //
-    // PRINTME("v", v);
-    //
-    // inverse_op_a.vmult(v, u);
-    // // identity.Tvmult(v2, u2);
-    //
-    // PRINTME("v", v);
   }
 }

diff --git a/tests/lac/linear_operator_14.with_cxx11=on.output b/tests/lac/linear_operator_14.with_cxx11=on.output
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..f4633ca2e29012dce47dc46dbfabe92d779d0789 100644 (file)
--- a/tests/lac/linear_operator_14.with_cxx11=on.output
+++ b/tests/lac/linear_operator_14.with_cxx11=on.output
@@ -0,0 +1,40 @@
+
+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