From 6771a0994131946dff274dea0a50c7603b24fa32 Mon Sep 17 00:00:00 2001 From: Denis Davydov 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::Tmmult(LAPACKFullMatrix &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 +#include +#include + +#include +#include + +DeclException5 (ExcEl, int, int, double, double,double, + << "Error in element (" << arg1 << "," << arg2 << "): " + << arg3 << " != " << arg4 << " delta=" << arg5); + +template +void test(const unsigned int n, const unsigned int k, const NumberType eps) +{ + deallog << n << " " << k << " " << std::endl; + FullMatrix A(k, n), C(n, n); + LAPACKFullMatrix 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> 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(el[0],el[1],1e-13); + deallog.pop(); + + deallog.push("float"); + for (auto el : sizes) + test(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