]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add LAPACKFullMatrix::reinit_preserve() and Vector::reinit_preserve()
authorDenis Davydov <davydden@gmail.com>
Fri, 15 Dec 2017 15:18:58 +0000 (16:18 +0100)
committerDenis Davydov <davydden@gmail.com>
Tue, 19 Dec 2017 12:27:13 +0000 (13:27 +0100)
to (partly) keep the old values.

doc/news/changes/minor/20171216DenisDavydov [new file with mode: 0644]
include/deal.II/lac/lapack_full_matrix.h
include/deal.II/lac/vector.h
include/deal.II/lac/vector.templates.h
source/lac/lapack_full_matrix.cc
tests/lapack/full_matrix_26.cc [new file with mode: 0644]
tests/lapack/full_matrix_26.output [new file with mode: 0644]
tests/vector/vector_58.cc [new file with mode: 0644]
tests/vector/vector_58.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20171216DenisDavydov b/doc/news/changes/minor/20171216DenisDavydov
new file mode 100644 (file)
index 0000000..abfcd93
--- /dev/null
@@ -0,0 +1,4 @@
+New: Add Vector::reinit_preserve() and LAPACKFullMatrix::reinit_preserve()
+to (partly) keep the previous values upon resizing.
+<br>
+(Denis Davydov, 2017/12/16)
index a1be34d185b241aeaeef99256f2e14047f846518..38ca6d253b1715a8ba5e2e0485088ce56b18bf31 100644 (file)
@@ -157,6 +157,33 @@ public:
    */
   void reinit (const size_type size);
 
+  /**
+   * Same as above but will preserve the values of matrix upon resizing.
+   * The original values
+   * of the matrix are kept on increasing the size
+   * \f[
+   * \mathbf A \rightarrow
+   * \left(
+   * \begin{array}{cc}
+   * \mathbf A & \mathbf 0 \\
+   * \mathbf 0 & \mathbf 0
+   * \end{array}
+   * \right)
+   * \f]
+   * Whereas if the new size is smaller, the matrix will contain the upper left block
+   * of the original one
+   * \f[
+   * \mathbf A_{11} \leftarrow
+   * \left(
+   * \begin{array}{cc}
+   * \mathbf A_{11} & \mathbf A_{12} \\
+   * \mathbf A_{21} & \mathbf A_{22}
+   * \end{array}
+   * \right)
+   * \f]
+   */
+  void reinit_preserve (const size_type size);
+
   /**
    * Regenerate the current matrix by one that has the same properties as if
    * it were created by the constructor of this class with the same argument
index acccbd1458a526d3e2e8c889929debeb4e897bca..eba1f004d8003e2166327981d63f5ed1ed97e935 100644 (file)
@@ -258,6 +258,32 @@ public:
   virtual void reinit (const size_type N,
                        const bool      omit_zeroing_entries=false);
 
+  /**
+   * Same as above, but will preserve the values of vector upon resizing.
+   * If we new size is bigger, we will have
+   * \f[
+   * \mathbf V \rightarrow
+   * \left(
+   * \begin{array}{c}
+   * \mathbf V   \\
+   * \mathbf 0
+   * \end{array}
+   * \right)
+   * \f]
+   * whereas if the desired size is smaller, then
+   * \f[
+   * \mathbf V_1 \leftarrow
+   * \left(
+   * \begin{array}{c}
+   * \mathbf V_1   \\
+   * \mathbf V_2
+   * \end{array}
+   * \right)
+   * \f]
+   */
+  void reinit_preserve (const size_type N);
+
+
   /**
    * Change the dimension to that of the vector @p V. The same applies as for
    * the other @p reinit function.
@@ -924,9 +950,10 @@ private:
 
   /**
    * Allocate and align @p values along 64-byte boundaries. The size of the
-   * allocated memory is determined by @p max_vec_size .
+   * allocated memory is determined by @p max_vec_size . Copy first
+   * @p copy_n_el from the old values.
    */
-  void allocate();
+  void allocate(const size_type copy_n_el = 0);
 };
 
 /*@}*/
index 19eaafcf05b010b16df9385549c03c4a91388fc5..2a1c83a2fa062e88fd2955aa2c69fbcca9a100b0 100644 (file)
@@ -289,6 +289,43 @@ void Vector<Number>::reinit (const size_type n,
 
 
 
+template <typename Number>
+inline
+void Vector<Number>::reinit_preserve (const size_type n)
+{
+  if (n==0)
+    {
+      values.reset();
+      max_vec_size = vec_size = 0;
+      thread_loop_partitioner.reset(new parallel::internal::TBBPartitioner());
+      return;
+    }
+
+  const size_type s = std::min(vec_size,n);
+  if (n>max_vec_size)
+    {
+      max_vec_size = n;
+      allocate(s);
+    }
+
+  if (vec_size != n)
+    {
+      vec_size = n;
+
+      // only reset the partitioner if we actually expect a significant vector
+      // size
+      if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size)
+        thread_loop_partitioner.reset(new parallel::internal::TBBPartitioner());
+    }
+
+  // pad with zeroes
+  for (size_type i = s; i < vec_size; ++i)
+    values[i] = Number();
+
+}
+
+
+
 template <typename Number>
 template <typename Number2>
 void Vector<Number>::reinit (const Vector<Number2> &v,
@@ -986,11 +1023,14 @@ Vector<Number>::memory_consumption () const
 
 template <typename Number>
 void
-Vector<Number>::allocate()
+Vector<Number>::allocate(const size_type copy_n_el)
 {
   // allocate memory with the proper alignment requirements of 64 bytes
   Number *new_values;
   Utilities::System::posix_memalign ((void **)&new_values, 64, sizeof(Number)*max_vec_size);
+  // copy:
+  for (size_type i = 0; i < copy_n_el; ++i)
+    new_values[i] = values[i];
   values.reset (new_values);
 }
 
index a12e833979f44a50b3576f4f04f30fef4ad42489..0ac15627e7c605f23d4e9f4528dc727df0eebc5c 100644 (file)
@@ -68,6 +68,7 @@ LAPACKFullMatrix<number>::operator = (const LAPACKFullMatrix<number> &M)
 }
 
 
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::reinit (const size_type n)
@@ -77,6 +78,23 @@ LAPACKFullMatrix<number>::reinit (const size_type n)
 }
 
 
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::reinit_preserve (const size_type n)
+{
+  TransposeTable<number> copy(*this);
+  const size_type s = std::min(std::min(this->m(), n), this->n());
+  this->TransposeTable<number>::reinit (n, n);
+  for (unsigned int i = 0; i < s; ++i)
+    for (unsigned int j = 0; j < s; ++j)
+      (*this)(i,j) = copy(i,j);
+
+  state = LAPACKSupport::matrix;
+}
+
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::reinit (const size_type m,
diff --git a/tests/lapack/full_matrix_26.cc b/tests/lapack/full_matrix_26.cc
new file mode 100644 (file)
index 0000000..5bf7500
--- /dev/null
@@ -0,0 +1,79 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 2015-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.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests reinitialization of square and rectangle LAPACKFullMatrix
+
+#include "../tests.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+
+#include <iostream>
+
+
+void test (const unsigned int size)
+{
+  AssertThrow (size>2, ExcInternalError());
+  const unsigned int smaller = size - 2;
+  const unsigned int larger  = size + 3;
+
+  // initialise a first matrix with the standard constructor and fill
+  // it with some numbers
+  LAPACKFullMatrix<double> M (size, size);
+
+  for (unsigned int i=0; i<size; ++i)
+    for (unsigned int j=0; j<size; ++j)
+      M(i,j) = i+2.*j;
+
+  LAPACKFullMatrix<double> M1(M), M2(M);
+
+  M1.reinit_preserve(smaller);
+  for (unsigned int i = 0; i < smaller; ++i)
+    for (unsigned int j = 0; j < smaller; ++j)
+      AssertThrow(M1(i,j) == M(i,j),
+                  ExcInternalError());
+
+  M2.reinit_preserve(larger);
+  for (unsigned int i = 0; i < larger; ++i)
+    for (unsigned int j = 0; j < larger; ++j)
+      {
+        if (i < size && j < size)
+          {
+            AssertThrow(M2(i,j) == M(i,j),
+                        ExcInternalError());
+          }
+        else
+          {
+            AssertThrow(M2(i,j) == 0.,
+                        ExcInternalError());
+          }
+      }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  test (4);
+  test (7);
+  test (11);
+  test (31);
+}
diff --git a/tests/lapack/full_matrix_26.output b/tests/lapack/full_matrix_26.output
new file mode 100644 (file)
index 0000000..5f42bb2
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/vector/vector_58.cc b/tests/vector/vector_58.cc
new file mode 100644 (file)
index 0000000..76fbb63
--- /dev/null
@@ -0,0 +1,79 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check Vector<double>::size()
+
+#include "../tests.h"
+#include <deal.II/lac/vector.h>
+
+static unsigned int counter = 0;
+
+template <typename Number>
+void fill (Vector<Number> &v)
+{
+  v = 0;
+  for (unsigned int i = 0; i < v.size(); i++)
+    v(i) = counter + i*2;
+
+  ++counter;
+}
+
+
+template <typename Number>
+void test (const std::vector<unsigned int> &size_sequence)
+{
+  Vector<Number> v, v_old;
+  for (unsigned int j = 0; j < size_sequence.size(); ++j)
+    {
+      const unsigned int s = size_sequence[j];
+      if (v.size() == 0)
+        {
+          v.reinit(s);
+        }
+      else
+        {
+          const unsigned int check_s = (s > v.size() ? v.size() : s);
+          v.reinit_preserve(s);
+          for (unsigned int i = 0; i < check_s; ++i)
+            AssertThrow(v(i) == v_old(i),
+                        ExcMessage("s=" + std::to_string(s) + " i=" + std::to_string(i) + " " +
+                                   std::to_string(v(i)) + "!=" + std::to_string(v_old(i))));
+
+          for (unsigned int i = check_s; i < s; ++i)
+            AssertThrow(v(i) == 0.,
+                        ExcMessage("s=" + std::to_string(s) + " i=" + std::to_string(i) + " " +
+                                   std::to_string(v(i)) + "!=0"));
+
+        }
+
+      fill(v);
+      v_old.reinit(s);
+      v_old = v;
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  initlog();
+
+  const std::vector<unsigned int> size_sequence = {{1,2,3,5,4,7,9,7,11,15,220,19,18,17,16,40,35,129,300,287}};
+  test<double>(size_sequence);
+}
diff --git a/tests/vector/vector_58.output b/tests/vector/vector_58.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

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.