From: Denis Davydov Date: Fri, 15 Dec 2017 15:18:58 +0000 (+0100) Subject: add LAPACKFullMatrix::reinit_preserve() and Vector::reinit_preserve() X-Git-Tag: v9.0.0-rc1~638^2~9 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=93fdd7a8575712b2c80a9bbd332b0ca7a27453d5;p=dealii.git add LAPACKFullMatrix::reinit_preserve() and Vector::reinit_preserve() to (partly) keep the old values. --- diff --git a/doc/news/changes/minor/20171216DenisDavydov b/doc/news/changes/minor/20171216DenisDavydov new file mode 100644 index 0000000000..abfcd93b8e --- /dev/null +++ b/doc/news/changes/minor/20171216DenisDavydov @@ -0,0 +1,4 @@ +New: Add Vector::reinit_preserve() and LAPACKFullMatrix::reinit_preserve() +to (partly) keep the previous values upon resizing. +
+(Denis Davydov, 2017/12/16) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index a1be34d185..38ca6d253b 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -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 diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index acccbd1458..eba1f004d8 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -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); }; /*@}*/ diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index 19eaafcf05..2a1c83a2fa 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -289,6 +289,43 @@ void Vector::reinit (const size_type n, +template +inline +void Vector::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 template void Vector::reinit (const Vector &v, @@ -986,11 +1023,14 @@ Vector::memory_consumption () const template void -Vector::allocate() +Vector::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); } diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index a12e833979..0ac15627e7 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -68,6 +68,7 @@ LAPACKFullMatrix::operator = (const LAPACKFullMatrix &M) } + template void LAPACKFullMatrix::reinit (const size_type n) @@ -77,6 +78,23 @@ LAPACKFullMatrix::reinit (const size_type n) } + +template +void +LAPACKFullMatrix::reinit_preserve (const size_type n) +{ + TransposeTable copy(*this); + const size_type s = std::min(std::min(this->m(), n), this->n()); + this->TransposeTable::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 void LAPACKFullMatrix::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 index 0000000000..5bf75008f5 --- /dev/null +++ b/tests/lapack/full_matrix_26.cc @@ -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 + +#include + + +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 M (size, size); + + for (unsigned int i=0; i 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 index 0000000000..5f42bb230f --- /dev/null +++ b/tests/lapack/full_matrix_26.output @@ -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 index 0000000000..76fbb63ae9 --- /dev/null +++ b/tests/vector/vector_58.cc @@ -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::size() + +#include "../tests.h" +#include + +static unsigned int counter = 0; + +template +void fill (Vector &v) +{ + v = 0; + for (unsigned int i = 0; i < v.size(); i++) + v(i) = counter + i*2; + + ++counter; +} + + +template +void test (const std::vector &size_sequence) +{ + Vector 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 size_sequence = {{1,2,3,5,4,7,9,7,11,15,220,19,18,17,16,40,35,129,300,287}}; + test(size_sequence); +} diff --git a/tests/vector/vector_58.output b/tests/vector/vector_58.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/vector/vector_58.output @@ -0,0 +1,2 @@ + +DEAL::OK