From 0f996063ce14aa50bc18ce9656ab6bfff1dfd6b1 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Mon, 9 Oct 2017 11:06:34 +0200 Subject: [PATCH] add LAPACKFullMatrix::trace() --- doc/news/changes/minor/20171009DenisDavydov | 3 ++ include/deal.II/lac/lapack_full_matrix.h | 6 +++ source/lac/lapack_full_matrix.cc | 19 +++++++ tests/lapack/full_matrix_20.cc | 57 +++++++++++++++++++++ tests/lapack/full_matrix_20.output | 5 ++ 5 files changed, 90 insertions(+) create mode 100644 doc/news/changes/minor/20171009DenisDavydov create mode 100644 tests/lapack/full_matrix_20.cc create mode 100644 tests/lapack/full_matrix_20.output diff --git a/doc/news/changes/minor/20171009DenisDavydov b/doc/news/changes/minor/20171009DenisDavydov new file mode 100644 index 0000000000..10c2aa6080 --- /dev/null +++ b/doc/news/changes/minor/20171009DenisDavydov @@ -0,0 +1,3 @@ +New: add LAPACKFullMatrix::trace() to calculate trace of a matrix. +
+(Denis Davydov, 2017/10/09) diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index a1a2c413a5..ede5e3e48c 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -459,6 +459,12 @@ public: */ number frobenius_norm() const; + /** + * Compute trace of the matrix, i.e. the sum of the diagonal values. + * Obviously, the matrix needs to be quadratic for this function. + */ + number trace() const; + /** * Invert the matrix by first computing an LU/Cholesky factorization with the LAPACK * function Xgetrf/Xpotrf and then building the actual inverse using Xgetri/Xpotri. diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index dc78147416..e298960064 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -594,6 +594,25 @@ number LAPACKFullMatrix::norm(const char type) const } + +template +number LAPACKFullMatrix::trace() const +{ + Assert (state == LAPACKSupport::matrix || + state == LAPACKSupport::inverse_matrix, + ExcMessage("Trace can be called in matrix state only.")); + Assert (this->n_cols() == this->n_rows(), + ExcDimensionMismatch(this->n_cols(), this->n_rows())); + + number tr = 0; + for (size_type i=0; in_rows(); ++i) + tr += (*this)(i,i); + + return tr; +} + + + template void LAPACKFullMatrix::compute_cholesky_factorization() diff --git a/tests/lapack/full_matrix_20.cc b/tests/lapack/full_matrix_20.cc new file mode 100644 index 0000000000..363acef6fb --- /dev/null +++ b/tests/lapack/full_matrix_20.cc @@ -0,0 +1,57 @@ +// --------------------------------------------------------------------- +// +// 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::trace() by comparing to FullMatrix + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + // Full matrix: + FullMatrix F(size); + create_random(F); + + // Lapack: + LAPACKFullMatrix M(size); + M = F; + + deallog << "trace difference: " << (F.trace() - M.trace())<< 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,391}}; + for (const auto &s : sizes) + { + deallog << "size=" << s << std::endl; + test(s); + } +} diff --git a/tests/lapack/full_matrix_20.output b/tests/lapack/full_matrix_20.output new file mode 100644 index 0000000000..2ea8f8c44c --- /dev/null +++ b/tests/lapack/full_matrix_20.output @@ -0,0 +1,5 @@ + +DEAL::size=17 +DEAL::trace difference: 0.00000 +DEAL::size=391 +DEAL::trace difference: 0.00000 -- 2.39.5