From 9b6dc3d72d8b55d64d855eeeedbcdce962f4096c Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 3 Jul 2017 10:12:30 +0200 Subject: [PATCH] Fix FullMatrix left_ and right_invert for square matrices --- ...170702AndreasKergassnerJean-PaulPelteret-2 | 4 + include/deal.II/lac/full_matrix.templates.h | 22 +++ tests/full_matrix/full_matrix_invert_01.cc | 127 ++++++++++++++++++ .../full_matrix_invert_01.with_lapack.output | 3 + 4 files changed, 156 insertions(+) create mode 100644 doc/news/changes/minor/20170702AndreasKergassnerJean-PaulPelteret-2 create mode 100644 tests/full_matrix/full_matrix_invert_01.cc create mode 100644 tests/full_matrix/full_matrix_invert_01.with_lapack.output diff --git a/doc/news/changes/minor/20170702AndreasKergassnerJean-PaulPelteret-2 b/doc/news/changes/minor/20170702AndreasKergassnerJean-PaulPelteret-2 new file mode 100644 index 0000000000..1ef69e6d4d --- /dev/null +++ b/doc/news/changes/minor/20170702AndreasKergassnerJean-PaulPelteret-2 @@ -0,0 +1,4 @@ +Fixed: FullMatrix::left_invert() and FullMatrix::right_invert() no longer +fail when the matrix is square. +
+(Andreas Kergassner, Jean-Paul Pelteret, 2017/07/02) diff --git a/include/deal.II/lac/full_matrix.templates.h b/include/deal.II/lac/full_matrix.templates.h index f086bae1c0..3582c72761 100644 --- a/include/deal.II/lac/full_matrix.templates.h +++ b/include/deal.II/lac/full_matrix.templates.h @@ -1541,6 +1541,17 @@ void FullMatrix::left_invert (const FullMatrix &A) { Assert (!A.empty(), ExcEmptyMatrix()); + + // If the matrix is square, simply do a + // standard inversion + if (A.m() == A.n()) + { + FullMatrix left_inv(A.n(),A.m()); + left_inv.invert(A); + *this = std::move(left_inv); + return; + } + Assert(A.m()>A.n(), ExcDimensionMismatch(A.m(), A.n())); Assert(this->m()==A.n(), ExcDimensionMismatch(this->m(), A.n())); Assert(this->n()==A.m(), ExcDimensionMismatch(this->n(), A.m())); @@ -1569,6 +1580,17 @@ void FullMatrix::right_invert (const FullMatrix &A) { Assert (!A.empty(), ExcEmptyMatrix()); + + // If the matrix is square, simply do a + // standard inversion + if (A.m() == A.n()) + { + FullMatrix right_inv(A.n(),A.m()); + right_inv.invert(A); + *this = std::move(right_inv); + return; + } + Assert(A.n()>A.m(), ExcDimensionMismatch(A.n(), A.m())); Assert(this->m()==A.n(), ExcDimensionMismatch(this->m(), A.n())); Assert(this->n()==A.m(), ExcDimensionMismatch(this->n(), A.m())); diff --git a/tests/full_matrix/full_matrix_invert_01.cc b/tests/full_matrix/full_matrix_invert_01.cc new file mode 100644 index 0000000000..03de111673 --- /dev/null +++ b/tests/full_matrix/full_matrix_invert_01.cc @@ -0,0 +1,127 @@ +// --------------------------------------------------------------------- +// +// 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 left and right inversion of FullMatrix + +#include "../tests.h" +#include "full_matrix_common.h" +#include + + +using namespace dealii; + + +template +void +fill_matrix_invertible(FullMatrix &A) +{ + for (unsigned int i=0; i +bool +calculate(const FullMatrix A, + bool disp = true) +{ + bool retval = true; + // Different tolerance for different number types + const number tol = 1000*std::numeric_limits::epsilon(); + + // Test left invert + if (A.m()>=A.n()) + { + FullMatrix A_l_inv(A.n(),A.m()); + FullMatrix identity(A.n(),A.n()); + A_l_inv.left_invert(A); + A_l_inv.mmult(identity,A); + + FullMatrix M(IdentityMatrix(identity.n())); + M.add(-1, identity); + if (disp || M.linfty_norm() > tol) + { + // deallog << "A matrix" << std::endl; + // display_matrix(A); + deallog << "Left inverse" << std::endl; + display_matrix(A_l_inv); + // deallog << "Identity = A_l_inv*A" << std::endl; + // display_matrix(identity); + retval = false; + } + } + + // Test right invert + if (A.m()<=A.n()) + { + FullMatrix A_r_inv(A.n(),A.m()); + FullMatrix identity(A.m(),A.m()); + A_r_inv.right_invert(A); + A.mmult(identity,A_r_inv); + + FullMatrix M(IdentityMatrix(identity.n())); + M.add(-1, identity); + if (disp || M.linfty_norm() > tol) + { + // deallog << "A matrix"<< std::endl; + // display_matrix(A); + deallog << "Right inverse"<< std::endl; + display_matrix(A_r_inv); + // deallog << "Identity = A*A_r_inv" << std::endl; + // display_matrix(identity); + // deallog << std::endl; + retval = false; + } + } + + return retval; +} + + +template +void +check () +{ + int nFails = 0; + int maxDim = 10; + for (unsigned int i=1; i <= maxDim; ++i) + { + for (unsigned int j=1; j <= maxDim; ++j) + { + FullMatrix A(i,j); + fill_matrix_invertible(A); + if (!calculate(A, false)) + { + nFails++; + } + } + } + + if (nFails > 0) + deallog + << nFails + << " out of " + << maxDim *maxDim + << " calls with wrong result" + << std::endl; + else + deallog << "OK" << std::endl; +} diff --git a/tests/full_matrix/full_matrix_invert_01.with_lapack.output b/tests/full_matrix/full_matrix_invert_01.with_lapack.output new file mode 100644 index 0000000000..8b3b075900 --- /dev/null +++ b/tests/full_matrix/full_matrix_invert_01.with_lapack.output @@ -0,0 +1,3 @@ + +DEAL::OK +DEAL::OK -- 2.39.5