make mTmult more clever if A==B

diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index 83fc1c2f8bdab22e4eb98e37ceb7e819fa0734bf..95ea3d6a1d729d4906e4eab08847c807544ab1a9 100644 (file)
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -765,8 +765,22 @@ LAPACKFullMatrix<number>::mTmult(LAPACKFullMatrix<number>       &C,
   const number alpha = 1.;
   const number beta = (adding ? 1. : 0.);
 
-  gemm("N", "T", &mm, &nn, &kk, &alpha, &this->values[0], &mm, &B.values[0],
-       &nn, &beta, &C.values[0], &mm);
+  if (PointerComparison::equal (this, &B))
+    {
+      syrk(&LAPACKSupport::U,"N",&nn,&kk,&alpha,&this->values[0],&nn,&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("N", "T", &mm, &nn, &kk, &alpha, &this->values[0], &mm, &B.values[0],
+           &nn, &beta, &C.values[0], &mm);
+    }
 }
 
 

diff --git a/tests/lapack/full_matrix_34.cc b/tests/lapack/full_matrix_34.cc
new file mode 100644 (file)
index 0000000..c86f298
--- /dev/null
+++ b/tests/lapack/full_matrix_34.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(n,k), C(n, n);
+  LAPACKFullMatrix<NumberType> AL(n,k), CL(n, n);
+
+  create_random(AL);
+
+  A = AL;
+
+  A.mTmult(C,A);
+  AL.mTmult(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.mTmult(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_34.output b/tests/lapack/full_matrix_34.output
new file mode 100644 (file)
index 0000000..8f1026c
--- /dev/null
+++ b/tests/lapack/full_matrix_34.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