From 1d7f8f56659afdfa2e1d8978ca2663c274fc6df1 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 18 Dec 2017 15:20:21 +0100 Subject: [PATCH] implement rank-1 update for Cholesky factorization --- source/lac/lapack_full_matrix.cc | 101 +++++++++++++++++++++----- tests/lapack/full_matrix_28.cc | 110 +++++++++++++++++++++++++++++ tests/lapack/full_matrix_28.output | 11 +++ 3 files changed, 205 insertions(+), 17 deletions(-) create mode 100644 tests/lapack/full_matrix_28.cc create mode 100644 tests/lapack/full_matrix_28.output diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 3e8246ce8a..c2007fd7cb 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -22,6 +22,8 @@ #include #include +#include + #include #include @@ -211,13 +213,69 @@ LAPACKFullMatrix::add (const number a, +namespace +{ + template + void cholesky_rank1(LAPACKFullMatrix &A, + const number a, + const Vector &v) + { + const unsigned int N = A.n(); + Vector z(v); + // Cholesky update / downdate, see + // 6.5.4 Cholesky Updating and Downdating, Golub 2013 Matrix computations + // Note that potrf() is called with LAPACKSupport::L , so the + // factorization is stored in lower triangular part. + if (a > 0.) + { + // simple update via a sequence of Givens rotations. + // Observe that + // + // | L^T |T | L^T | + // A <-- | | | | = L L^T + z z^T + // | z^T | | z^T | + // + // so we can get updated factor by doing a sequence of Givens + // rotations to make the matrix lower-triangular + // Also see LINPACK's dchud http://www.netlib.org/linpack/dchud.f + z *= std::sqrt(a); + for (unsigned int k = 0; k < N; ++k) + { + const std::array csr = Utilities::LinearAlgebra::Givens_rotation(A(k,k),z(k)); + A(k,k) = csr[2]; + for (unsigned int i = k+1; i < N; ++i) + { + const number t = A(i,k); + A(i,k) = csr[0] * A(i,k) + csr[1] * z(i); + z(i) = -csr[1] * t + csr[0] * z(i); + + //A(i,k) = ( t + csr[1] * z(i)) / csr[0]; + //z(i) = - csr[1] * t + csr[0] * z(i) ; + } + } + } + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + + + template + void cholesky_rank1(LAPACKFullMatrix> &A, + const std::complex a, + const Vector> &v) + { + AssertThrow(false, ExcNotImplemented()); + } +} + + template void LAPACKFullMatrix::add(const number a, const Vector &v) { - Assert(state == LAPACKSupport::matrix, - ExcState(state)); Assert(property == LAPACKSupport::symmetric, ExcProperty(property)); @@ -226,19 +284,28 @@ LAPACKFullMatrix::add(const number a, AssertIsFinite(a); - const int N = this->n_rows(); - const char uplo = LAPACKSupport::U; - const int lda = N; - const int incx=1; - - syr(&uplo, &N, &a, v.begin(), &incx, this->values.begin(), &lda); - - // FIXME: we should really only update upper or lower triangular parts - // of symmetric matrices and make sure the interface is consistent, - // for example operator(i,j) gives correct results regardless of storage. - for (unsigned int i = 0; i < N; ++i) - for (unsigned int j = 0; j < i; ++j) - (*this)(i,j) = (*this)(j,i); + if (state == LAPACKSupport::matrix) + { + const int N = this->n_rows(); + const char uplo = LAPACKSupport::U; + const int lda = N; + const int incx=1; + + syr(&uplo, &N, &a, v.begin(), &incx, this->values.begin(), &lda); + + // FIXME: we should really only update upper or lower triangular parts + // of symmetric matrices and make sure the interface is consistent, + // for example operator(i,j) gives correct results regardless of storage. + for (unsigned int i = 0; i < N; ++i) + for (unsigned int j = 0; j < i; ++j) + (*this)(i,j) = (*this)(j,i); + } + else if (state == LAPACKSupport::cholesky) + { + cholesky_rank1(*this, a, v); + } + else + AssertThrow(false, ExcState(state)); } @@ -1143,7 +1210,7 @@ LAPACKFullMatrix::compute_eigenvalues(const bool right, if (info != 0) std::cerr << "LAPACK error in geev" << std::endl; - state = LAPACKSupport::State(eigenvalues | unusable); + state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable); } @@ -1396,7 +1463,7 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( eigenvectors[i](j) = values_A[col_begin+j]; } } - state = LAPACKSupport::State(eigenvalues | unusable); + state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable); } diff --git a/tests/lapack/full_matrix_28.cc b/tests/lapack/full_matrix_28.cc new file mode 100644 index 0000000000..4ddba3ce37 --- /dev/null +++ b/tests/lapack/full_matrix_28.cc @@ -0,0 +1,110 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::add() for rank1 update of a Cholesky factorization + +/* MWE in Octave: +A = pascal(4) +A = + + 1 1 1 1 + 1 2 3 4 + 1 3 6 10 + 1 4 10 20 + +x = [3 2 1 5]'; + +R = chol(A) +R = + + 1 1 1 1 + 0 1 2 3 + 0 0 1 3 + 0 0 0 1 + +R1 = cholupdate(R,x) + +R1 = + + 3.16228 2.21359 1.26491 5.05964 + 0.00000 1.04881 2.09762 2.66970 + 0.00000 0.00000 1.00000 3.00000 + 0.00000 0.00000 0.00000 1.80907 + +*/ + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test() +{ + const unsigned int size = 4; + LAPACKFullMatrix A(size); + A.set_property(LAPACKSupport::symmetric); + Vector v(size); + + A(0,0) = 1; + A(0,1) = 1; + A(0,2) = 1; + A(0,3) = 1; + A(1,0) = 1; + A(1,1) = 2; + A(1,2) = 3; + A(1,3) = 4; + A(2,0) = 1; + A(2,1) = 3; + A(2,2) = 6; + A(2,3) = 10; + A(3,0) = 1; + A(3,1) = 4; + A(3,2) =10; + A(3,3) = 20; + + v(0) = 3; + v(1) = 2; + v(2) = 1; + v(3) = 5; + + const NumberType a = 1.; + + A.compute_cholesky_factorization(); + deallog << "Cholesky:" << std::endl; + A.print_formatted(deallog.get_file_stream(),5); + + A.add(a,v); + + deallog << "Cholesky updated:" << std::endl; + A.print_formatted(deallog.get_file_stream(),5); +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(5); + deallog.attach(logfile); + + test(); +} diff --git a/tests/lapack/full_matrix_28.output b/tests/lapack/full_matrix_28.output new file mode 100644 index 0000000000..6d36c7d29d --- /dev/null +++ b/tests/lapack/full_matrix_28.output @@ -0,0 +1,11 @@ + +DEAL::Cholesky: +1.00000e+00 +1.00000e+00 1.00000e+00 +1.00000e+00 2.00000e+00 1.00000e+00 +1.00000e+00 3.00000e+00 3.00000e+00 1.00000e+00 +DEAL::Cholesky updated: +3.16228e+00 +2.21359e+00 1.04881e+00 +1.26491e+00 2.09762e+00 1.00000e+00 +5.05964e+00 2.66970e+00 3.00000e+00 1.80907e+00 -- 2.39.5