From 6771a0994131946dff274dea0a50c7603b24fa32 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Thu, 15 Feb 2018 17:15:25 +0100
Subject: [PATCH] make Tmmult more clever if A==B

---
 source/lac/lapack_full_matrix.cc   | 18 +++++-
 tests/lapack/full_matrix_33.cc     | 88 ++++++++++++++++++++++++++++++
 tests/lapack/full_matrix_33.output | 49 +++++++++++++++++
 3 files changed, 153 insertions(+), 2 deletions(-)
 create mode 100644 tests/lapack/full_matrix_33.cc
 create mode 100644 tests/lapack/full_matrix_33.output

diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index 2073e5489d..83fc1c2f8b 100644
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -704,8 +704,22 @@ LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number>       &C,
   const number alpha = 1.;
   const number beta = (adding ? 1. : 0.);
 
-  gemm("T", "N", &mm, &nn, &kk, &alpha, &this->values[0], &kk, &B.values[0],
-       &kk, &beta, &C.values[0], &mm);
+  if (PointerComparison::equal (this, &B))
+    {
+      syrk(&LAPACKSupport::U,"T",&nn,&kk,&alpha,&this->values[0],&kk,&beta,&C.values[0],&nn);
+
+      // fill-in lower triangular part
+      for (types::blas_int j = 0; j < nn; ++j)
+        for (types::blas_int i = 0; i < j; ++i)
+          C(j,i) = C(i,j);
+
+      C.property = symmetric;
+    }
+  else
+    {
+      gemm("T", "N", &mm, &nn, &kk, &alpha, &this->values[0], &kk, &B.values[0],
+           &kk, &beta, &C.values[0], &mm);
+    }
 }
 
 
diff --git a/tests/lapack/full_matrix_33.cc b/tests/lapack/full_matrix_33.cc
new file mode 100644
index 0000000000..8f58e369f5
--- /dev/null
+++ b/tests/lapack/full_matrix_33.cc
@@ -0,0 +1,88 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests LAPACKFullMatrix::Tmmult with  A==B which shall use Xsyrk instead Xgemm
+
+#include "../tests.h"
+#include "create_matrix.h"
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+#include <tuple>
+
+DeclException5 (ExcEl, int, int, double, double,double,
+                << "Error in element (" << arg1 << "," << arg2 << "): "
+                << arg3 << " != " << arg4 << " delta=" << arg5);
+
+template <typename NumberType>
+void test(const unsigned int n, const unsigned int k, const NumberType eps)
+{
+  deallog << n << " " << k << " " << std::endl;
+  FullMatrix<NumberType> A(k, n), C(n, n);
+  LAPACKFullMatrix<NumberType> AL(k,n), CL(n, n);
+
+  create_random(AL);
+
+  A = AL;
+
+  A.Tmmult(C,A);
+  AL.Tmmult(CL,AL);
+
+  for (unsigned int i=0; i<n; ++i)
+    for (unsigned int j=0; j<n; ++j)
+      AssertThrow(std::abs(C(i,j)-CL(i,j)) < eps * std::abs(CL(i,j)),
+                  ExcEl(i,j,C(i,j),CL(i,j),C(i,j)-CL(i,j))
+                 );
+
+  deallog << "OK" << std::endl;
+
+  AL.Tmmult(CL,AL,true);
+
+  for (unsigned int i=0; i<n; ++i)
+    for (unsigned int j=0; j<n; ++j)
+      AssertThrow(std::abs(2.*C(i,j)-CL(i,j)) < eps * std::abs(CL(i,j)),
+                  ExcEl(i,j,2.*C(i,j),CL(i,j),2.*C(i,j)-CL(i,j))
+                 );
+
+  deallog << "OK adding" << std::endl;
+}
+
+int main()
+{
+  const std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  logfile.precision(3);
+  deallog.attach(logfile);
+
+  const std::vector<std::array<unsigned int,2>> sizes =
+  {
+    {3,3}, {7,7}, {51,51}, {320,320},
+    {3,9}, {9,7}, {10,5},  {320,120}
+  };
+
+  deallog.push("double");
+  for (auto el : sizes)
+    test<double>(el[0],el[1],1e-13);
+  deallog.pop();
+
+  deallog.push("float");
+  for (auto el : sizes)
+    test<float>(el[0],el[1],1e-5);
+  deallog.pop();
+
+}
diff --git a/tests/lapack/full_matrix_33.output b/tests/lapack/full_matrix_33.output
new file mode 100644
index 0000000000..8f1026ce09
--- /dev/null
+++ b/tests/lapack/full_matrix_33.output
@@ -0,0 +1,49 @@
+
+DEAL:double::3 3 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:double::7 7 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:double::51 51 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:double::320 320 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:double::3 9 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:double::9 7 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:double::10 5 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:double::320 120 
+DEAL:double::OK
+DEAL:double::OK adding
+DEAL:float::3 3 
+DEAL:float::OK
+DEAL:float::OK adding
+DEAL:float::7 7 
+DEAL:float::OK
+DEAL:float::OK adding
+DEAL:float::51 51 
+DEAL:float::OK
+DEAL:float::OK adding
+DEAL:float::320 320 
+DEAL:float::OK
+DEAL:float::OK adding
+DEAL:float::3 9 
+DEAL:float::OK
+DEAL:float::OK adding
+DEAL:float::9 7 
+DEAL:float::OK
+DEAL:float::OK adding
+DEAL:float::10 5 
+DEAL:float::OK
+DEAL:float::OK adding
+DEAL:float::320 120 
+DEAL:float::OK
+DEAL:float::OK adding
-- 
2.39.5