-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.
<br>
(Denis Davydov, 2017/12/16)
-New: LAPACKFullMatrix::add(const number, const Vector<number> &) performs
+New: LAPACKFullMatrix::rank1_update(const number, const Vector<number> &) performs
a rank-1 update of a matrix.
<br>
(Denis Davydov, 2017/12/17)
* 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<number> &v);
+ void rank1_update(const number a,
+ const Vector<number> &v);
/**
* Assignment from different matrix classes, performing the usual conversion
* 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
{
/**
- * 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
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
namespace LinearAlgebra
{
+ template<typename NumberType>
+ std::array<std::complex<NumberType>,3> hyperbolic_rotation(const std::complex<NumberType> &f,
+ const std::complex<NumberType> &g)
+ {
+ AssertThrow(false, ExcNotImplemented());
+ std::array<NumberType,3> res;
+ return res;
+ }
+
+
+
template<typename NumberType>
std::array<NumberType,3> hyperbolic_rotation(const NumberType &f,
const NumberType &g)
return csr;
}
- template<typename NumberType>
- std::array<std::complex<NumberType>,3> hyperbolic_rotation(const std::complex<NumberType> &f,
- const std::complex<NumberType> &g)
- {
- AssertThrow(false, ExcNotImplemented());
- std::array<NumberType,3> res;
- return res;
- }
template<typename NumberType>
return res;
}
+
+
template<typename NumberType>
std::array<NumberType,3> Givens_rotation(const NumberType &f,
const NumberType &g)
{
std::array<NumberType,3> 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())
* \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
template <typename Number>
inline
-void Vector<Number>::reinit_preserve (const size_type n)
+void Vector<Number>::grow_or_shrink (const size_type n)
{
if (n==0)
{
template <typename number>
void
-LAPACKFullMatrix<number>::reinit_preserve (const size_type n)
+LAPACKFullMatrix<number>::grow_or_shrink (const size_type n)
{
TransposeTable<number> copy(*this);
const size_type s = std::min(std::min(this->m(), n), this->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>::add(const number a,
- const Vector<number> &v)
+LAPACKFullMatrix<number>::rank1_update(const number a,
+ const Vector<number> &v)
{
Assert(property == LAPACKSupport::symmetric,
ExcProperty(property));
-
template <typename number>
void
LAPACKFullMatrix<number>::vmult (
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());
deallog << std::setprecision(6);
deallog.attach(logfile);
- test<double>(1., 0.);
- test<double>(0., 1.);
- test<double>(1., -2.);
- test<double>(-1., 2.);
- test<double>(-15.,3.);
- test<double>(15.,-3.);
+ deallog << "Residuals:"<< std::endl;
+ // g==0
+ test<double>( 3., 0.);
+ test<double>(-2., 0.);
+ // f==0
+ test<double>( 0., 2.);
+ test<double>( 0., -5.);
+ // |f| > |g|
+ test<double>( 15., 3.);
+ test<double>( 15.,-4.);
+ test<double>(-17., 2.);
+ test<double>(-18.,-5.);
+ // |f| < |g|
+ test<double>( 2., -4.);
+ test<double>(-2., 3.);
+ test<double>(-3., -7.);
+ test<double>( 5., 9.);
}
+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
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());
deallog << std::setprecision(6);
deallog.attach(logfile);
- // check all combinations with real solutions:
- test<double>( 1., 0.); // g == 0
- test<double>( 2., 1.); // both positive
+ deallog << "Residuals:"<< std::endl;
+ // check all combinations with real solutions: |f| > |g|
+ test<double>( 3., 0.); // g == 0
+ test<double>( 2., 1.5); // both positive
test<double>( 3., -0.5); // g negative
test<double>(-4., -2.4); // both negative
test<double>(-5., 2); // f negative
+DEAL::Residuals:
DEAL::0.00000
DEAL::0.00000
DEAL::0.00000
// ---------------------------------------------------------------------
-// Tests reinitialization of square and rectangle LAPACKFullMatrix
+// Tests grow_or_shrink() of square and rectangle LAPACKFullMatrix
#include "../tests.h"
#include <deal.II/lac/lapack_full_matrix.h>
LAPACKFullMatrix<double> 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)
{
// ---------------------------------------------------------------------
//
-// 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.
//
// ---------------------------------------------------------------------
-// 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"
LAPACKFullMatrix<NumberType> 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)
ExcMessage("diff="+ std::to_string(diff)));
}
- deallog << "Ok" << std::endl;
+ deallog << "OK" << std::endl;
}
DEAL::size=17
-DEAL::Ok
+DEAL::OK
DEAL::size=35
-DEAL::Ok
+DEAL::OK
DEAL::size=391
-DEAL::Ok
+DEAL::OK
// ---------------------------------------------------------------------
-// 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)
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);
// ---------------------------------------------------------------------
-// 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)
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);
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) + " " +
projects / dealii.git / commitdiff