to (partly) keep the old values.
--- /dev/null
+New: Add Vector::reinit_preserve() and LAPACKFullMatrix::reinit_preserve()
+to (partly) keep the previous values upon resizing.
+<br>
+(Denis Davydov, 2017/12/16)
*/
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
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.
/**
* 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);
};
/*@}*/
+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,
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);
}
}
+
template <typename number>
void
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,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+}
--- /dev/null
+
+DEAL::OK