From fa28a4675b9281003bbd66445742ba31781ad5c6 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 6 Dec 2017 13:47:44 +0100 Subject: [PATCH] lapack: add reciprocal_condition_number() for triangular matrices --- include/deal.II/lac/lapack_full_matrix.h | 6 ++ include/deal.II/lac/lapack_templates.h | 71 +++++++++++++++++++++++ source/lac/lapack_full_matrix.cc | 32 +++++++++++ tests/lapack/full_matrix_22.cc | 73 ++++++++++++++++++++++++ tests/lapack/full_matrix_22.output | 7 +++ 5 files changed, 189 insertions(+) create mode 100644 tests/lapack/full_matrix_22.cc create mode 100644 tests/lapack/full_matrix_22.output diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index 7d75ba796b..83660f69b9 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -444,6 +444,12 @@ public: */ number reciprocal_condition_number(const number l1_norm) const; + /** + * Same as above but for triangular matrices. The matrix has to have + * the LAPACKSupport::Property set, see set_property(). + */ + number reciprocal_condition_number() const; + /** * Compute the determinant of a matrix. As it requires the LU factorization of * the matrix, this function can only be called after diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index 956e232e82..a48ee3739a 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -83,6 +83,16 @@ extern "C" void spocon_ (const char *uplo, const int *n, const float *A, const int *lda, const float *anorm, float *rcond, float *work, int *iwork, int *info); +// Estimate the reciprocal of the condition number of triangular matrices +// http://www.netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_gaff914510b1673e90752c095f5b9dcedf.html#gaff914510b1673e90752c095f5b9dcedf + void dtrcon_ (const char *norm, const char *uplo, const char *diag, + const int *n, const double *A, const int *lda, + double *rcond, + double *work, int *iwork, int *info); + void strcon_ (const char *norm, const char *uplo, const char *diag, + const int *n, const float *A, const int *lda, + float *rcond, + float *work, int *iwork, int *info); // Computes the inverse from Cholesky void dpotri_ (const char *uplo, const int *n, double *A, const int *lda, int *info); @@ -491,6 +501,67 @@ potrf (const char *uplo, const int *n, float *A, const int *lda, int *info) +/// Tempalte wrapper for trcon +template +inline void +trcon(const char *norm, const char *uplo, const char *diag, + const int *n, const number *A, const int *lda, + number *rcond, + number *work, int *iwork, int *info) +{ + Assert (false, ExcNotImplemented()); +} + +inline void +trcon (const char *norm, const char *uplo, const char *diag, + const int *n, const double *A, const int *lda, + double *rcond, + double *work, int *iwork, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + dtrcon_ (norm, uplo, diag, n, A, lda, rcond, work, iwork, info); +#else + (void) norm; + (void) uplo; + (void) diag; + (void) n; + (void) A; + (void) lda; + (void) rcond; + (void) work; + (void) iwork; + (void) info; + + Assert (false, LAPACKSupport::ExcMissing("dtrcon")); +#endif +} + +inline void +trcon (const char *norm, const char *uplo, const char *diag, + const int *n, const float *A, const int *lda, + float *rcond, + float *work, int *iwork, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + strcon_ (norm, uplo, diag, n, A, lda, rcond, work, iwork, info); +#else + (void) norm; + (void) uplo; + (void) diag; + (void) n; + (void) A; + (void) lda; + (void) rcond; + (void) work; + (void) iwork; + (void) info; + + Assert (false, LAPACKSupport::ExcMissing("dtrcon")); +#endif +} + + + /// Template wrapper for pocon template inline void diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index f53a0bb94d..344f09ec57 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -689,6 +689,38 @@ LAPACKFullMatrix::reciprocal_condition_number(const number a_norm) const +template +number +LAPACKFullMatrix::reciprocal_condition_number() const +{ + Threads::Mutex::ScopedLock lock (mutex); + Assert(property == upper_triangular || + property == lower_triangular, + ExcProperty(property)); + number rcond = 0.; + + const int N = this->n_rows(); + const number *values = &this->values[0]; + int info = 0; + const int lda = std::max(1,N); + work.resize(3*N); + iwork.resize(N); + + const char norm = '1'; + const char diag = 'N'; + const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L); + trcon(&norm, &uplo, &diag, + &N, values, &lda, + &rcond, + work.data(), iwork.data(), &info); + + Assert(info >= 0, ExcInternalError()); + + return rcond; +} + + + template void LAPACKFullMatrix::compute_svd() diff --git a/tests/lapack/full_matrix_22.cc b/tests/lapack/full_matrix_22.cc new file mode 100644 index 0000000000..8f94ff1585 --- /dev/null +++ b/tests/lapack/full_matrix_22.cc @@ -0,0 +1,73 @@ +// --------------------------------------------------------------------- +// +// 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::reciprocal_condition_number() for triangular matrices + +/* MWE for size=3 in Octave: +R = [1,2,3; 0, 5, 6; 0, 0, 9] +rcond(R) +ans = 0.055556 +*/ + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + LAPACKFullMatrix M(size); + M.set_property(LAPACKSupport::upper_triangular); + + M = 0.; + unsigned int counter = 1; + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j = 0; j < size; ++j) + { + if (j >= i) + M(i,j) = counter; + + counter++; + } + + // M.print_formatted(deallog.get_file_stream(), 3, false, 8); + const double rcond = M.reciprocal_condition_number(); + + deallog << rcond << 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 = {{1,3,11}}; + for (const auto &s : sizes) + { + deallog << "size=" << s << std::endl; + // test(s); + test(s); + } +} diff --git a/tests/lapack/full_matrix_22.output b/tests/lapack/full_matrix_22.output new file mode 100644 index 0000000000..1ea812f1cf --- /dev/null +++ b/tests/lapack/full_matrix_22.output @@ -0,0 +1,7 @@ + +DEAL::size=1 +DEAL::1.00000 +DEAL::size=3 +DEAL::0.0555556 +DEAL::size=11 +DEAL::0.00137741 -- 2.39.5