From: Denis Davydov Date: Sun, 17 Dec 2017 16:06:15 +0000 (+0100) Subject: add rank-1 update of LAPACKFullMatrix X-Git-Tag: v9.0.0-rc1~638^2~8 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9497504a19d0247446cc22dce5ddda26e72dec24;p=dealii.git add rank-1 update of LAPACKFullMatrix --- diff --git a/doc/news/changes/minor/20171217DenisDavydov b/doc/news/changes/minor/20171217DenisDavydov new file mode 100644 index 0000000000..b08111702d --- /dev/null +++ b/doc/news/changes/minor/20171217DenisDavydov @@ -0,0 +1,4 @@ +New: LAPACKFullMatrix::add(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 38ca6d253b..f63a7dae68 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -141,6 +141,13 @@ public: void add (const number a, const LAPACKFullMatrix &A); + /** + * Perform a rank-1 update of a symmetric matrix + * $ A \leftarrow A + a \, \rm v \rm v^T $. + */ + void add(const number a, + const Vector &v); + /** * Assignment from different matrix classes, performing the usual conversion * to the transposed format expected by LAPACK. This assignment operator diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index c7cc243bc7..dcda155951 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -297,11 +297,51 @@ extern "C" float *d, float *e, float *z, const int *ldz, float *work, int *info); +// Rank-1 update for symmetric matrices + void dsyr_ (const char *uplo, const int *n, + const double *alpha, const double *x, + const int *incx, + double *A, const int *lda); + void ssyr_ (const char *uplo, const int *n, + const float *alpha, const float *x, + const int *incx, + float *A, const int *lda); } DEAL_II_NAMESPACE_OPEN +/// Template wrapper for LAPACK functions dsyr and ssyr +template +inline void +syr(const char *, const int *, + const number *, const number *, + const int *, + number *, const int *) +{ + Assert (false, ExcNotImplemented()); +} + +#ifdef DEAL_II_WITH_LAPACK +inline void +syr(const char *uplo, const int *n, + const double *alpha, const double *x, + const int *incx, + double *A, const int *lda) +{ + dsyr_(uplo,n,alpha,x,incx,A,lda); +} +inline void +syr(const char *uplo, const int *n, + const float *alpha, const float *x, + const int *incx, + float *A, const int *lda) +{ + ssyr_(uplo,n,alpha,x,incx,A,lda); +} +#endif + + /// Template wrapper for LAPACK functions daxpy and saxpy template diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 0ac15627e7..bed672572a 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -211,6 +211,39 @@ LAPACKFullMatrix::add (const number a, +template +void +LAPACKFullMatrix::add(const number a, + const Vector &v) +{ + Assert(state == LAPACKSupport::matrix, + ExcState(state)); + Assert(property == LAPACKSupport::symmetric, + ExcProperty(property)); + + Assert (n() == m(), ExcInternalError()); + Assert (m() == v.size(), ExcDimensionMismatch(m(), v.size())); + + 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); +} + + + + template void LAPACKFullMatrix::vmult ( diff --git a/tests/lapack/full_matrix_27.cc b/tests/lapack/full_matrix_27.cc new file mode 100644 index 0000000000..ce9d7aab7b --- /dev/null +++ b/tests/lapack/full_matrix_27.cc @@ -0,0 +1,80 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 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. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::add() for rank1 update of a matrix + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + // Lapack: + LAPACKFullMatrix A(size); + A.set_property(LAPACKSupport::symmetric); + Vector v(size); + for (unsigned int i = 0; i < size; ++i) + { + v(i) = random_value(); + for (unsigned int j = i; j < size; ++j) + { + const NumberType val = random_value(); + A(i,j) = val; + A(j,i) = val; + } + } + + + const NumberType a = random_value(); + + LAPACKFullMatrix B(size); + B = A; + + B.add(a, v); + + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j = 0; j < size; ++j) + { + const NumberType diff = A(i,j) + a * v(i) * v(j) - B(i,j); + AssertThrow(std::abs(diff) < 1e-10 * std::abs(B(i,j)) , + ExcMessage("diff="+ std::to_string(diff))); + } + + deallog << "Ok" << std::endl; +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + const std::vector sizes = {{17,35,391}}; + for (const auto &s : sizes) + { + deallog << "size=" << s << std::endl; + test(s); + } +} diff --git a/tests/lapack/full_matrix_27.output b/tests/lapack/full_matrix_27.output new file mode 100644 index 0000000000..2f2ffd5b67 --- /dev/null +++ b/tests/lapack/full_matrix_27.output @@ -0,0 +1,7 @@ + +DEAL::size=17 +DEAL::Ok +DEAL::size=35 +DEAL::Ok +DEAL::size=391 +DEAL::Ok