From dbfb8a94f607c7b7f7354fd5c37030c42812ac15 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 18 Dec 2017 22:53:55 +0100 Subject: [PATCH] various cleanups and improvements --- doc/news/changes/minor/20171216DenisDavydov | 2 +- doc/news/changes/minor/20171217DenisDavydov | 2 +- include/deal.II/lac/lapack_full_matrix.h | 17 ++++++----- include/deal.II/lac/utilities.h | 33 +++++++++++++-------- include/deal.II/lac/vector.h | 5 ++-- include/deal.II/lac/vector.templates.h | 2 +- source/lac/lapack_full_matrix.cc | 10 +++---- tests/lac/utilities_04.cc | 25 +++++++++++----- tests/lac/utilities_04.output | 7 +++++ tests/lac/utilities_05.cc | 9 +++--- tests/lac/utilities_05.output | 1 + tests/lapack/full_matrix_26.cc | 6 ++-- tests/lapack/full_matrix_27.cc | 8 ++--- tests/lapack/full_matrix_27.output | 6 ++-- tests/lapack/full_matrix_28.cc | 4 +-- tests/lapack/full_matrix_29.cc | 6 ++-- tests/vector/vector_58.cc | 2 +- 17 files changed, 87 insertions(+), 58 deletions(-) diff --git a/doc/news/changes/minor/20171216DenisDavydov b/doc/news/changes/minor/20171216DenisDavydov index abfcd93b8e..7ef45c916e 100644 --- a/doc/news/changes/minor/20171216DenisDavydov +++ b/doc/news/changes/minor/20171216DenisDavydov @@ -1,4 +1,4 @@ -New: Add Vector::reinit_preserve() and LAPACKFullMatrix::reinit_preserve() +New: Add Vector::grow_or_shrink() and LAPACKFullMatrix::grow_or_shrink() to (partly) keep the previous values upon resizing.
(Denis Davydov, 2017/12/16) diff --git a/doc/news/changes/minor/20171217DenisDavydov b/doc/news/changes/minor/20171217DenisDavydov index b08111702d..b1b8ab4c53 100644 --- a/doc/news/changes/minor/20171217DenisDavydov +++ b/doc/news/changes/minor/20171217DenisDavydov @@ -1,4 +1,4 @@ -New: LAPACKFullMatrix::add(const number, const Vector &) performs +New: LAPACKFullMatrix::rank1_update(const number, const Vector &) performs a rank-1 update of a matrix.
(Denis Davydov, 2017/12/17) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index 0479b3ba0c..32cff48404 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -145,13 +145,16 @@ public: * Perform a rank-1 update of a symmetric matrix * $ A \leftarrow A + a \, \rm v \rm v^T $. * - * This function also works for Cholesky factorization. Updating ($a>0$) is + * This function also works for Cholesky factorization. + * In that case, updating ($a>0$) is * performed via Givens rotations, whereas downdating ($a<0$) via hyperbolic - * rotations. Note that the later case might lead to a negative definite - * matrix in which case the error will be thrown. + * rotations. Note that the latter case might lead to a negative definite + * matrix in which case the error will be thrown (because Cholesky + * factorizations are only valid for symmetric and positive definite + * matrices). */ - void add(const number a, - const Vector &v); + void rank1_update(const number a, + const Vector &v); /** * Assignment from different matrix classes, performing the usual conversion @@ -185,16 +188,16 @@ public: * 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} + * \rightarrow \mathbf A_{11} * \right) * \f] */ - void reinit_preserve (const size_type size); + void grow_or_shrink (const size_type size); /** * Regenerate the current matrix by one that has the same properties as if diff --git a/include/deal.II/lac/utilities.h b/include/deal.II/lac/utilities.h index b1c56f5bf4..352d00ad72 100644 --- a/include/deal.II/lac/utilities.h +++ b/include/deal.II/lac/utilities.h @@ -37,7 +37,8 @@ namespace Utilities { /** - * Return elements of continuous Givens rotation matrix and norm of the input vector. + * Return the elements of a continuous Givens rotation matrix and + * the norm of the input vector. * * That is for a given * pair @p x and @p y, return $c$ , $s$ and $\sqrt{x^2+y^2}$ such that @@ -66,7 +67,7 @@ namespace Utilities const NumberType &y); /** - * Return elements of hyperbolic rotation matrix. + * Return the elements of a hyperbolic rotation matrix. * * That is for a given * pair @p x and @p y, return $c$ , $s$ and $r$ such that @@ -209,6 +210,17 @@ namespace Utilities namespace LinearAlgebra { + template + std::array,3> hyperbolic_rotation(const std::complex &f, + const std::complex &g) + { + AssertThrow(false, ExcNotImplemented()); + std::array res; + return res; + } + + + template std::array hyperbolic_rotation(const NumberType &f, const NumberType &g) @@ -229,14 +241,6 @@ namespace Utilities return csr; } - template - std::array,3> hyperbolic_rotation(const std::complex &f, - const std::complex &g) - { - AssertThrow(false, ExcNotImplemented()); - std::array res; - return res; - } template @@ -248,19 +252,24 @@ namespace Utilities return res; } + + template std::array Givens_rotation(const NumberType &f, const NumberType &g) { std::array res; - // naieve calculation for "r" may overflow or underflow: + // naive calculation for "r" may overflow or underflow: // c = x / \sqrt{x^2+y^2} // s = -y / \sqrt{x^2+y^2} // See Golub 2013, Matrix computations, Chapter 5.1.8 // Algorithm 5.1.3 // and - // Anderson (2000), Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem. LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, December 4, 2000. + // Anderson (2000), + // Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem. + // LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, + // December 4, 2000. // Algorithm 4 // We implement the latter below: if (g == NumberType()) diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index eba1f004d8..23172d6067 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -272,17 +272,16 @@ public: * \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) + * \rightarrow \mathbf V_1 * \f] */ - void reinit_preserve (const size_type N); - + void grow_or_shrink (const size_type N); /** * Change the dimension to that of the vector @p V. The same applies as for diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index 2a1c83a2fa..04b03b18d9 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -291,7 +291,7 @@ void Vector::reinit (const size_type n, template inline -void Vector::reinit_preserve (const size_type n) +void Vector::grow_or_shrink (const size_type n) { if (n==0) { diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 5645f64f4e..f94cbbeee2 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -83,7 +83,7 @@ LAPACKFullMatrix::reinit (const size_type n) template void -LAPACKFullMatrix::reinit_preserve (const size_type n) +LAPACKFullMatrix::grow_or_shrink (const size_type n) { TransposeTable copy(*this); const size_type s = std::min(std::min(this->m(), n), this->n()); @@ -91,8 +91,6 @@ LAPACKFullMatrix::reinit_preserve (const size_type 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; } @@ -305,10 +303,11 @@ namespace } + template void -LAPACKFullMatrix::add(const number a, - const Vector &v) +LAPACKFullMatrix::rank1_update(const number a, + const Vector &v) { Assert(property == LAPACKSupport::symmetric, ExcProperty(property)); @@ -344,7 +343,6 @@ LAPACKFullMatrix::add(const number a, - template void LAPACKFullMatrix::vmult ( diff --git a/tests/lac/utilities_04.cc b/tests/lac/utilities_04.cc index bdff2d07fa..078658c9c4 100644 --- a/tests/lac/utilities_04.cc +++ b/tests/lac/utilities_04.cc @@ -46,7 +46,7 @@ void test (const NumberType a, const NumberType b) const double norm = res.l2_norm(); deallog << norm << std::endl; - if (norm > 1e-12) + if (norm > 1e-12 || csr[2] < 0.) { deallog << "x:" << std::endl; x.print(deallog.get_file_stream()); @@ -67,10 +67,21 @@ int main() deallog << std::setprecision(6); deallog.attach(logfile); - test(1., 0.); - test(0., 1.); - test(1., -2.); - test(-1., 2.); - test(-15.,3.); - test(15.,-3.); + deallog << "Residuals:"<< std::endl; + // g==0 + test( 3., 0.); + test(-2., 0.); + // f==0 + test( 0., 2.); + test( 0., -5.); + // |f| > |g| + test( 15., 3.); + test( 15.,-4.); + test(-17., 2.); + test(-18.,-5.); + // |f| < |g| + test( 2., -4.); + test(-2., 3.); + test(-3., -7.); + test( 5., 9.); } diff --git a/tests/lac/utilities_04.output b/tests/lac/utilities_04.output index 687ddbebbc..2a9b0035fe 100644 --- a/tests/lac/utilities_04.output +++ b/tests/lac/utilities_04.output @@ -1,4 +1,11 @@ +DEAL::Residuals: +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 +DEAL::0.00000 DEAL::0.00000 DEAL::0.00000 DEAL::0.00000 diff --git a/tests/lac/utilities_05.cc b/tests/lac/utilities_05.cc index 92839492af..88a31ebb6d 100644 --- a/tests/lac/utilities_05.cc +++ b/tests/lac/utilities_05.cc @@ -46,7 +46,7 @@ void test (const NumberType a, const NumberType b) const NumberType norm = res.l2_norm(); deallog << norm << std::endl; - if (norm > 1e-12) + if (norm > 1e-12 || csr[2] < 0.) { deallog << "x:" << std::endl; x.print(deallog.get_file_stream()); @@ -67,9 +67,10 @@ int main() deallog << std::setprecision(6); deallog.attach(logfile); - // check all combinations with real solutions: - test( 1., 0.); // g == 0 - test( 2., 1.); // both positive + deallog << "Residuals:"<< std::endl; + // check all combinations with real solutions: |f| > |g| + test( 3., 0.); // g == 0 + test( 2., 1.5); // both positive test( 3., -0.5); // g negative test(-4., -2.4); // both negative test(-5., 2); // f negative diff --git a/tests/lac/utilities_05.output b/tests/lac/utilities_05.output index 1f17d1a21b..35a82f964e 100644 --- a/tests/lac/utilities_05.output +++ b/tests/lac/utilities_05.output @@ -1,4 +1,5 @@ +DEAL::Residuals: DEAL::0.00000 DEAL::0.00000 DEAL::0.00000 diff --git a/tests/lapack/full_matrix_26.cc b/tests/lapack/full_matrix_26.cc index 5bf75008f5..62067264d4 100644 --- a/tests/lapack/full_matrix_26.cc +++ b/tests/lapack/full_matrix_26.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// Tests reinitialization of square and rectangle LAPACKFullMatrix +// Tests grow_or_shrink() of square and rectangle LAPACKFullMatrix #include "../tests.h" #include @@ -38,13 +38,13 @@ void test (const unsigned int size) LAPACKFullMatrix M1(M), M2(M); - M1.reinit_preserve(smaller); + M1.grow_or_shrink(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); + M2.grow_or_shrink(larger); for (unsigned int i = 0; i < larger; ++i) for (unsigned int j = 0; j < larger; ++j) { diff --git a/tests/lapack/full_matrix_27.cc b/tests/lapack/full_matrix_27.cc index ce9d7aab7b..64e62d469b 100644 --- a/tests/lapack/full_matrix_27.cc +++ b/tests/lapack/full_matrix_27.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2005 - 2015 by the deal.II authors +// Copyright (C) 2017 by the deal.II authors // // This file is part of the deal.II library. // @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// test LAPACKFullMatrix::add() for rank1 update of a matrix +// test LAPACKFullMatrix::rank1_update() for rank1 update of a matrix #include "../tests.h" #include "create_matrix.h" @@ -50,7 +50,7 @@ test(const unsigned int size) LAPACKFullMatrix B(size); B = A; - B.add(a, v); + B.rank1_update(a, v); for (unsigned int i = 0; i < size; ++i) for (unsigned int j = 0; j < size; ++j) @@ -60,7 +60,7 @@ test(const unsigned int size) ExcMessage("diff="+ std::to_string(diff))); } - deallog << "Ok" << std::endl; + deallog << "OK" << std::endl; } diff --git a/tests/lapack/full_matrix_27.output b/tests/lapack/full_matrix_27.output index 2f2ffd5b67..cee1d11396 100644 --- a/tests/lapack/full_matrix_27.output +++ b/tests/lapack/full_matrix_27.output @@ -1,7 +1,7 @@ DEAL::size=17 -DEAL::Ok +DEAL::OK DEAL::size=35 -DEAL::Ok +DEAL::OK DEAL::size=391 -DEAL::Ok +DEAL::OK diff --git a/tests/lapack/full_matrix_28.cc b/tests/lapack/full_matrix_28.cc index 4ddba3ce37..643872b3c2 100644 --- a/tests/lapack/full_matrix_28.cc +++ b/tests/lapack/full_matrix_28.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// test LAPACKFullMatrix::add() for rank1 update of a Cholesky factorization +// test LAPACKFullMatrix::rank1_update() for rank1 update of a Cholesky factorization /* MWE in Octave: A = pascal(4) @@ -92,7 +92,7 @@ test() deallog << "Cholesky:" << std::endl; A.print_formatted(deallog.get_file_stream(),5); - A.add(a,v); + A.rank1_update(a,v); deallog << "Cholesky updated:" << std::endl; A.print_formatted(deallog.get_file_stream(),5); diff --git a/tests/lapack/full_matrix_29.cc b/tests/lapack/full_matrix_29.cc index f7569b6947..912cea1031 100644 --- a/tests/lapack/full_matrix_29.cc +++ b/tests/lapack/full_matrix_29.cc @@ -14,7 +14,7 @@ // --------------------------------------------------------------------- -// test LAPACKFullMatrix::add() for rank1 downdate of a Cholesky factorization +// test LAPACKFullMatrix::rank1_update() for rank1 downdate of a Cholesky factorization /* MWE in Octave: A = pascal(4) @@ -103,14 +103,14 @@ test() deallog << "Cholesky:" << std::endl; A.print_formatted(deallog.get_file_stream(),5); - A.add(a,v); + A.rank1_update(a,v); deallog << "Cholesky updated:" << std::endl; A.print_formatted(deallog.get_file_stream(),5); v /= 2.; - A.add(-1.,v); + A.rank1_update(-1.,v); deallog << "Cholesky downdated:" << std::endl; A.print_formatted(deallog.get_file_stream(),5); diff --git a/tests/vector/vector_58.cc b/tests/vector/vector_58.cc index 76fbb63ae9..7967399d81 100644 --- a/tests/vector/vector_58.cc +++ b/tests/vector/vector_58.cc @@ -47,7 +47,7 @@ void test (const std::vector &size_sequence) else { const unsigned int check_s = (s > v.size() ? v.size() : s); - v.reinit_preserve(s); + v.grow_or_shrink(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) + " " + -- 2.39.5