From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Mon, 3 Jul 2017 08:12:30 +0000 (+0200)
Subject: Fix FullMatrix left_ and right_invert for square matrices
X-Git-Tag: v9.0.0-rc1~1451^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F4570%2Fhead;p=dealii.git

Fix FullMatrix left_ and right_invert for square matrices
---

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.
+<br>
+(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<number>::left_invert (const FullMatrix<number2> &A)
 {
   Assert (!A.empty(), ExcEmptyMatrix());
+
+  // If the matrix is square, simply do a
+  // standard inversion
+  if (A.m() == A.n())
+    {
+      FullMatrix<number2> 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<number>::right_invert (const FullMatrix<number2> &A)
 {
   Assert (!A.empty(), ExcEmptyMatrix());
+
+  // If the matrix is square, simply do a
+  // standard inversion
+  if (A.m() == A.n())
+    {
+      FullMatrix<number2> 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 <limits>
+
+
+using namespace dealii;
+
+
+template <typename number>
+void
+fill_matrix_invertible(FullMatrix<number> &A)
+{
+  for (unsigned int i=0; i<A.m(); i++)
+    for (unsigned int j=0; j<A.n(); j++)
+      {
+        A(i,j)=number(i*j);
+        if (i==j)
+          A(i,i)+=i+A.m();
+      }
+}
+
+
+template <typename number>
+bool
+calculate(const FullMatrix<number> A,
+          bool disp = true)
+{
+  bool retval = true;
+  // Different tolerance for different number types
+  const number tol = 1000*std::numeric_limits<number>::epsilon();
+
+  // Test left invert
+  if (A.m()>=A.n())
+    {
+      FullMatrix<number> A_l_inv(A.n(),A.m());
+      FullMatrix<number> identity(A.n(),A.n());
+      A_l_inv.left_invert(A);
+      A_l_inv.mmult(identity,A);
+
+      FullMatrix<double> 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<number> A_r_inv(A.n(),A.m());
+      FullMatrix<number> identity(A.m(),A.m());
+      A_r_inv.right_invert(A);
+      A.mmult(identity,A_r_inv);
+
+      FullMatrix<double> 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 <typename number>
+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<number> 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