From: Martin Kronbichler <kronbichler@lnm.mw.tum.de> Date: Tue, 2 Feb 2016 19:20:48 +0000 (+0100) Subject: Fix bug in matrix multiplication X-Git-Tag: v8.4.0-rc2~41^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e5b630adbf5b8f50cc6a583ad07af8d2db2d650b;p=dealii.git Fix bug in matrix multiplication --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 2d9235768a..4ee4086d67 100644 --- 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 1fbd74528d..bc6ea88fea 100644 --- 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 index 0000000000..cf98582b6e --- /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 index 0000000000..8d3ce8fd61 --- /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 index 0000000000..314f6e7788 --- /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 index 0000000000..75bb1cfa81 --- /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 index 0000000000..4635662fb0 --- /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 index 0000000000..75bb1cfa81 --- /dev/null +++ b/tests/lac/full_matrix_11.output @@ -0,0 +1,4 @@ + +DEAL::Difference: 0 +DEAL::Difference: 0 +DEAL::Difference: 0