Fix bug in matrix multiplication

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 2d9235768ab0deadf6d9981de26ded5dead87d65..4ee4086d67556d86f51e6bcc0f0821460331f687 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -500,6 +500,13 @@ inconvenience this causes.
 
 
 <ol>
+  <li> Fixed: FullMatrix::TmTmult for matrix multiplication used to compute
+  wrong results for larger matrix sizes where external BLAS is called. This
+  has been fixed.
+  <br>
+  (Martin Kronbichler, 2016/02/02)
+  </li>
+
   <li> Fixed: parallel::distributed::Triangulation with periodic boundary
   conditions did not respect 2:1 balance over vertices on periodic
   boundaries. This lead to incomplete ghost layers on multigrid levels. This

diff --git a/include/deal.II/lac/full_matrix.templates.h b/include/deal.II/lac/full_matrix.templates.h
index 1fbd74528d63cdd60f093317f9e82bd1deaf0cf2..bc6ea88fea54334e4bf0ab505fa1f2a8d434c0f5 100644 (file)
--- a/include/deal.II/lac/full_matrix.templates.h
+++ b/include/deal.II/lac/full_matrix.templates.h
@@ -809,7 +809,7 @@ void FullMatrix<number>::TmTmult (FullMatrix<number2>       &dst,
         // transpose matrices, and read the result as if it were row-wise
         // again. In other words, we calculate (B A)^T, which is A^T B^T.
 
-        const int m = src.n();
+        const int m = src.m();
         const int n = this->n();
         const int k = this->m();
         const char *trans = "t";

diff --git a/tests/lac/full_matrix_09.cc b/tests/lac/full_matrix_09.cc
new file mode 100644 (file)
index 0000000..cf98582
--- /dev/null
+++ b/tests/lac/full_matrix_09.cc
@@ -0,0 +1,64 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+//check method TmTmult of FullMatrix on larger size than full_matrix_03 where
+//we interface to the external BLAS by comparison to mmult of the transpose
+//matrix
+
+#include "../tests.h"
+
+#include <deal.II/lac/full_matrix.h>
+
+
+template <typename Number>
+void test()
+{
+  FullMatrix<Number> A(2, 76), B(76, 3), C(2, 3), D(3, 2), E(2, 3);
+  for (unsigned int i=0; i<A.m(); ++i)
+    for (unsigned int j=0; j<A.n(); ++j)
+      A(i,j) = (double)Testing::rand()/RAND_MAX;
+  for (unsigned int i=0; i<B.m(); ++i)
+    for (unsigned int j=0; j<B.n(); ++j)
+      B(i,j) = (double)Testing::rand()/RAND_MAX;
+
+  A.mmult(C, B); // C = A * B
+  B.TmTmult(D, A); // D = B^T * A^T
+
+  E.copy_transposed(D); // E = D^T = C
+
+  C.add(-1., E);
+
+  deallog << "Difference 1: " << (double)C.l1_norm() << std::endl;
+
+  C = 0;
+  for (unsigned int i=0; i<A.m(); ++i)
+    for (unsigned int j=0; j<B.n(); ++j)
+      for (unsigned int k=0; k<A.n(); ++k)
+        C(i,j) += A(i,k) * B(k,j);
+  C.add(-1., E);
+
+  deallog << "Difference 2: " << (double)C.l1_norm() << std::endl;
+}
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1e-5);
+
+  test<double>();
+  test<float>();
+  test<long double>();
+}

diff --git a/tests/lac/full_matrix_09.output b/tests/lac/full_matrix_09.output
new file mode 100644 (file)
index 0000000..8d3ce8f
--- /dev/null
+++ b/tests/lac/full_matrix_09.output
@@ -0,0 +1,7 @@
+
+DEAL::Difference 1: 0
+DEAL::Difference 2: 0
+DEAL::Difference 1: 0
+DEAL::Difference 2: 0
+DEAL::Difference 1: 0
+DEAL::Difference 2: 0

diff --git a/tests/lac/full_matrix_10.cc b/tests/lac/full_matrix_10.cc
new file mode 100644 (file)
index 0000000..314f6e7
--- /dev/null
+++ b/tests/lac/full_matrix_10.cc
@@ -0,0 +1,55 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+//check method mTmult of FullMatrix on larger size than full_matrix_02 where
+//we interface to the external BLAS
+
+#include "../tests.h"
+
+#include <deal.II/lac/full_matrix.h>
+
+
+template <typename Number>
+void test()
+{
+  FullMatrix<Number> A(2, 76), B(3, 76), C(2, 3), D(2, 3);
+  for (unsigned int i=0; i<A.m(); ++i)
+    for (unsigned int j=0; j<A.n(); ++j)
+      A(i,j) = (double)Testing::rand()/RAND_MAX;
+  for (unsigned int i=0; i<B.m(); ++i)
+    for (unsigned int j=0; j<B.n(); ++j)
+      B(i,j) = (double)Testing::rand()/RAND_MAX;
+
+  A.mTmult(C, B); // C = A * B^T
+
+  for (unsigned int i=0; i<A.m(); ++i)
+    for (unsigned int j=0; j<B.m(); ++j)
+      for (unsigned int k=0; k<A.n(); ++k)
+        D(i,j) += A(i,k) * B(j,k);
+  C.add(-1., D);
+
+  deallog << "Difference: " << (double)C.l1_norm() << std::endl;
+}
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1e-5);
+
+  test<double>();
+  test<float>();
+  test<long double>();
+}

diff --git a/tests/lac/full_matrix_10.output b/tests/lac/full_matrix_10.output
new file mode 100644 (file)
index 0000000..75bb1cf
--- /dev/null
+++ b/tests/lac/full_matrix_10.output
@@ -0,0 +1,4 @@
+
+DEAL::Difference: 0
+DEAL::Difference: 0
+DEAL::Difference: 0

diff --git a/tests/lac/full_matrix_11.cc b/tests/lac/full_matrix_11.cc
new file mode 100644 (file)
index 0000000..4635662
--- /dev/null
+++ b/tests/lac/full_matrix_11.cc
@@ -0,0 +1,55 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+//check method Tmmult of FullMatrix on larger size than full_matrix_02 where
+//we interface to the external BLAS
+
+#include "../tests.h"
+
+#include <deal.II/lac/full_matrix.h>
+
+
+template <typename Number>
+void test()
+{
+  FullMatrix<Number> A(76, 2), B(76, 3), C(2, 3), D(2, 3);
+  for (unsigned int i=0; i<A.m(); ++i)
+    for (unsigned int j=0; j<A.n(); ++j)
+      A(i,j) = (double)Testing::rand()/RAND_MAX;
+  for (unsigned int i=0; i<B.m(); ++i)
+    for (unsigned int j=0; j<B.n(); ++j)
+      B(i,j) = (double)Testing::rand()/RAND_MAX;
+
+  A.Tmmult(C, B); // C = A^T * B
+
+  for (unsigned int i=0; i<A.n(); ++i)
+    for (unsigned int j=0; j<B.n(); ++j)
+      for (unsigned int k=0; k<A.m(); ++k)
+        D(i,j) += A(k,i) * B(k,j);
+  C.add(-1., D);
+
+  deallog << "Difference: " << (double)C.l1_norm() << std::endl;
+}
+
+int main()
+{
+  initlog();
+  deallog.threshold_double(1e-5);
+
+  test<double>();
+  test<float>();
+  test<long double>();
+}

diff --git a/tests/lac/full_matrix_11.output b/tests/lac/full_matrix_11.output
new file mode 100644 (file)
index 0000000..75bb1cf
--- /dev/null
+++ b/tests/lac/full_matrix_11.output
@@ -0,0 +1,4 @@
+
+DEAL::Difference: 0
+DEAL::Difference: 0
+DEAL::Difference: 0