From b05f32fe1151b8c62cf94d06415420cd58b2526e Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 16 Feb 2018 18:03:54 +0100 Subject: [PATCH] add LAPACK::Tmmult() to do a tripple product with a diagonal matrix --- include/deal.II/lac/lapack_full_matrix.h | 25 +++++ include/deal.II/lac/lapack_templates.h | 43 ++++++++ source/lac/lapack_full_matrix.cc | 91 ++++++++++++++--- tests/lapack/create_matrix.h | 9 ++ tests/lapack/full_matrix_31.cc | 125 +++++++++++++++++++++++ tests/lapack/full_matrix_31.output | 73 +++++++++++++ tests/lapack/full_matrix_32.cc | 86 ++++++++++++++++ tests/lapack/full_matrix_32.output | 33 ++++++ 8 files changed, 473 insertions(+), 12 deletions(-) create mode 100644 tests/lapack/full_matrix_31.cc create mode 100644 tests/lapack/full_matrix_31.output create mode 100644 tests/lapack/full_matrix_32.cc create mode 100644 tests/lapack/full_matrix_32.output diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index a629d5d9a3..9608a08415 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -420,6 +420,25 @@ public: const LAPACKFullMatrix &B, const bool adding=false) const; + /** + * Matrix-matrix-multiplication using transpose of this and a + * diagonal vector @p V. + * + * If the adding=false then the result is stored in the matrix + * $C = A^T \rm{diag}(V) B$ + * otherwise it is added $C \mathrel{+}= A^T \rm{diag}(V) B$. + * + * @note It is assumed that @p A, @p B and @p V have compatible sizes and that + * @p C already has the right size. + * + * @note This function is not provided by LAPACK. The function first forms $BV$ product and + * then uses Xgemm function. + */ + void Tmmult (LAPACKFullMatrix &C, + const LAPACKFullMatrix &B, + const Vector &V, + const bool adding=false) const; + /** * Matrix-matrix-multiplication using transpose of B. * @@ -475,6 +494,12 @@ public: const LAPACKFullMatrix &B, const bool adding=false) const; + /** + * Scale rows of this matrix by @p V . This is equivalent to premultiplication + * with a diagonal matrix $A\leftarrow {\rm diag}(V)A$. + */ + void scale_rows(const Vector &V); + /** * Compute the LU factorization of the matrix using LAPACK function Xgetrf. */ diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index 4777872444..45932cdb4a 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -54,6 +54,15 @@ extern "C" const float *alpha, const float *A, const dealii::types::blas_int *lda, const float *B, const dealii::types::blas_int *ldb, const float *beta, float *C, const dealii::types::blas_int *ldc); +// Symmetric rank-k update + void dsyrk_ (const char *uplo, const char *trans, + const dealii::types::blas_int *n, const dealii::types::blas_int *k, + const double *alpha, const double *A, const dealii::types::blas_int *lda, + const double *beta, double *C, const dealii::types::blas_int *ldc); + void ssyrk_ (const char *uplo, const char *trans, + const dealii::types::blas_int *n, const dealii::types::blas_int *k, + const float *alpha, const float *A, const dealii::types::blas_int *lda, + const float *beta, float *C, const dealii::types::blas_int *ldc); // Compute LU factorization void dgetrf_ (const dealii::types::blas_int *m, const dealii::types::blas_int *n, double *A, const dealii::types::blas_int *lda, dealii::types::blas_int *ipiv, dealii::types::blas_int *info); @@ -311,6 +320,40 @@ extern "C" DEAL_II_NAMESPACE_OPEN + +/// Template wrapper for LAPACK functions dsyrk and ssyrk +template +inline void +syrk(const char *, const char *, + const types::blas_int *, const types::blas_int *, + const number *, const number *, const types::blas_int *, + const number *, number *, const types::blas_int *) +{ + Assert (false, ExcNotImplemented()); +} + +#ifdef DEAL_II_WITH_LAPACK +inline void +syrk(const char *uplo, const char *trans, + const types::blas_int *n, const types::blas_int *k, + const double *alpha, const double *A, const types::blas_int *lda, + const double *beta, double *C, const types::blas_int *ldc) +{ + dsyrk_(uplo,trans,n,k,alpha,A,lda,beta,C,ldc); +} + +inline void +syrk(const char *uplo, const char *trans, + const types::blas_int *n, const types::blas_int *k, + const float *alpha, const float *A, const types::blas_int *lda, + const float *beta, float *C, const types::blas_int *ldc) +{ + ssyrk_(uplo,trans,n,k,alpha,A,lda,beta,C,ldc); +} +#endif + + + /// Template wrapper for LAPACK functions dsyr and ssyr template inline void diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index afb761b880..2073e5489d 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -578,8 +578,8 @@ LAPACKFullMatrix::mmult(LAPACKFullMatrix &C, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); - Assert(C.state == matrix || C.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); + Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state)); Assert (this->n() == B.m(), ExcDimensionMismatch(this->n(), B.m())); Assert (C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n())); Assert (C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m())); @@ -601,7 +601,7 @@ LAPACKFullMatrix::mmult(FullMatrix &C, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); Assert (this->n() == B.m(), ExcDimensionMismatch(this->n(), B.m())); Assert (C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n())); Assert (C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m())); @@ -623,11 +623,78 @@ template void LAPACKFullMatrix::Tmmult(LAPACKFullMatrix &C, const LAPACKFullMatrix &B, + const Vector &V, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); - Assert(C.state == matrix || C.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); + Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state)); + + const LAPACKFullMatrix &A = *this; + + Assert (A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m())); + Assert (C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n())); + Assert (C.m() == A.n(), ExcDimensionMismatch(A.n(), C.m())); + Assert (V.size() == A.m(), ExcDimensionMismatch(A.m(), V.size())); + + const types::blas_int mm = A.n(); + const types::blas_int nn = B.n(); + const types::blas_int kk = B.m(); + + // lapack does not have any tripple product routines, including the case of + // diagonal matrix in the middle, see + // https://stackoverflow.com/questions/3548069/multiplying-three-matrices-in-blas-with-the-middle-one-being-diagonal + // http://mathforum.org/kb/message.jspa?messageID=3546564 + + Threads::Mutex::ScopedLock lock (mutex); + // First, get V*B into "work" array + work.resize(kk*nn); + // following http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=768#p2577 + // do left-multiplication manually. Note that Xscal would require to first + // copy the input vector as multiplication is done inplace. + for (types::blas_int j = 0; j < nn; ++j) + for (types::blas_int i = 0; i < kk; ++i) + { + Assert(j*kk+i(work.size()), ExcInternalError()); + work[j*kk+i]=V(i)*B(i,j); + } + + // Now do the standard Tmmult: + const number alpha = 1.; + const number beta = (adding ? 1. : 0.); + + gemm("T", "N", &mm, &nn, &kk, &alpha, &this->values[0], &kk, &work[0], + &kk, &beta, &C.values[0], &mm); +} + + + +template +void +LAPACKFullMatrix::scale_rows(const Vector &V) +{ + LAPACKFullMatrix &A = *this; + Assert(state == matrix || state == inverse_matrix, ExcState(state)); + Assert (V.size() == A.m(), ExcDimensionMismatch(A.m(), V.size())); + + const types::blas_int nn = A.n(); + const types::blas_int kk = A.m(); + for (types::blas_int j = 0; j < nn; ++j) + for (types::blas_int i = 0; i < kk; ++i) + A(i,j)*=V(i); +} + + + +template +void +LAPACKFullMatrix::Tmmult(LAPACKFullMatrix &C, + const LAPACKFullMatrix &B, + const bool adding) const +{ + Assert(state == matrix || state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); + Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state)); Assert (this->m() == B.m(), ExcDimensionMismatch(this->m(), B.m())); Assert (C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n())); Assert (C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m())); @@ -649,7 +716,7 @@ LAPACKFullMatrix::Tmmult(FullMatrix &C, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); Assert (this->m() == B.m(), ExcDimensionMismatch(this->m(), B.m())); Assert (C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n())); Assert (C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m())); @@ -673,8 +740,8 @@ LAPACKFullMatrix::mTmult(LAPACKFullMatrix &C, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); - Assert(C.state == matrix || C.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); + Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state)); Assert (this->n() == B.n(), ExcDimensionMismatch(this->n(), B.n())); Assert (C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m())); Assert (C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m())); @@ -697,7 +764,7 @@ LAPACKFullMatrix::mTmult(FullMatrix &C, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); Assert (this->n() == B.n(), ExcDimensionMismatch(this->n(), B.n())); Assert (C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m())); Assert (C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m())); @@ -721,8 +788,8 @@ LAPACKFullMatrix::TmTmult(LAPACKFullMatrix &C, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); - Assert(C.state == matrix || C.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); + Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state)); Assert (this->m() == B.n(), ExcDimensionMismatch(this->m(), B.n())); Assert (C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m())); Assert (C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m())); @@ -744,7 +811,7 @@ LAPACKFullMatrix::TmTmult(FullMatrix &C, const bool adding) const { Assert(state == matrix || state == inverse_matrix, ExcState(state)); - Assert(B.state == matrix || B.state == inverse_matrix, ExcState(state)); + Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state)); Assert (this->m() == B.n(), ExcDimensionMismatch(this->m(), B.n())); Assert (C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m())); Assert (C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m())); diff --git a/tests/lapack/create_matrix.h b/tests/lapack/create_matrix.h index a06bd006bb..6e47841d03 100644 --- a/tests/lapack/create_matrix.h +++ b/tests/lapack/create_matrix.h @@ -50,3 +50,12 @@ void create_random (FullMatrix &A) for (unsigned int j = 0; j < A.n(); ++j) A(i,j) = random_value(); } + + + +template +void create_random (Vector &V) +{ + for (unsigned int i = 0; i < V.size(); ++i) + V(i) = random_value(); +} diff --git a/tests/lapack/full_matrix_31.cc b/tests/lapack/full_matrix_31.cc new file mode 100644 index 0000000000..9f2917e170 --- /dev/null +++ b/tests/lapack/full_matrix_31.cc @@ -0,0 +1,125 @@ +// --------------------------------------------------------------------- +// +// 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 additional vector +// (similar to _09) + +#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 m, const unsigned int n, const unsigned int k, const NumberType eps) +{ + deallog << m << " " << n << " " << k << " " << std::endl; + FullMatrix A(k, m), B(k, n), C(m, n), D(k,k); + LAPACKFullMatrix AL(k,m), BL(k, n), CL(m, n); + Vector DL(k); + + create_random(AL); + create_random(BL); + create_random(DL); + + A = AL; + B = BL; + D = NumberType(); + for (unsigned int i=0; i> sizes = + { + {3,3,3}, {7,7,7}, {51,51,51}, {320,320,320}, + {3,5,9}, {7,7,9}, {10,5,10}, {320,320,120} + }; + + deallog.push("double"); + for (auto el : sizes) + test(el[0],el[1],el[2],1e-13); + deallog.pop(); + + deallog.push("float"); + for (auto el : sizes) + test(el[0],el[1],el[2],1e-5); + deallog.pop(); + +} diff --git a/tests/lapack/full_matrix_31.output b/tests/lapack/full_matrix_31.output new file mode 100644 index 0000000000..2bb351db3e --- /dev/null +++ b/tests/lapack/full_matrix_31.output @@ -0,0 +1,73 @@ + +DEAL:double::3 3 3 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::OK symmetric +DEAL:double::OK symmetric adding +DEAL:double::7 7 7 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::OK symmetric +DEAL:double::OK symmetric adding +DEAL:double::51 51 51 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::OK symmetric +DEAL:double::OK symmetric adding +DEAL:double::320 320 320 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::OK symmetric +DEAL:double::OK symmetric adding +DEAL:double::3 5 9 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::7 7 9 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::OK symmetric +DEAL:double::OK symmetric adding +DEAL:double::10 5 10 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::320 320 120 +DEAL:double::OK non-symmetric +DEAL:double::OK non-symmetric adding +DEAL:double::OK symmetric +DEAL:double::OK symmetric adding +DEAL:float::3 3 3 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::OK symmetric +DEAL:float::OK symmetric adding +DEAL:float::7 7 7 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::OK symmetric +DEAL:float::OK symmetric adding +DEAL:float::51 51 51 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::OK symmetric +DEAL:float::OK symmetric adding +DEAL:float::320 320 320 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::OK symmetric +DEAL:float::OK symmetric adding +DEAL:float::3 5 9 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::7 7 9 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::OK symmetric +DEAL:float::OK symmetric adding +DEAL:float::10 5 10 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::320 320 120 +DEAL:float::OK non-symmetric +DEAL:float::OK non-symmetric adding +DEAL:float::OK symmetric +DEAL:float::OK symmetric adding diff --git a/tests/lapack/full_matrix_32.cc b/tests/lapack/full_matrix_32.cc new file mode 100644 index 0000000000..7c0cce2ac6 --- /dev/null +++ b/tests/lapack/full_matrix_32.cc @@ -0,0 +1,86 @@ +// --------------------------------------------------------------------- +// +// 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::scale_rows() + +#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(k, n), D(k,k); + LAPACKFullMatrix AL(k,n); + Vector DL(k); + + create_random(AL); + create_random(DL); + + A = AL; + D = NumberType(); + for (unsigned int i=0; i> sizes = + { + {3,3}, {7,7}, {51,51}, {320,320}, + {3,9}, {9,7}, {5,17}, {320,121} + }; + + 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_32.output b/tests/lapack/full_matrix_32.output new file mode 100644 index 0000000000..c91d76fb3f --- /dev/null +++ b/tests/lapack/full_matrix_32.output @@ -0,0 +1,33 @@ + +DEAL:double::3 3 +DEAL:double::OK +DEAL:double::7 7 +DEAL:double::OK +DEAL:double::51 51 +DEAL:double::OK +DEAL:double::320 320 +DEAL:double::OK +DEAL:double::3 9 +DEAL:double::OK +DEAL:double::9 7 +DEAL:double::OK +DEAL:double::5 17 +DEAL:double::OK +DEAL:double::320 121 +DEAL:double::OK +DEAL:float::3 3 +DEAL:float::OK +DEAL:float::7 7 +DEAL:float::OK +DEAL:float::51 51 +DEAL:float::OK +DEAL:float::320 320 +DEAL:float::OK +DEAL:float::3 9 +DEAL:float::OK +DEAL:float::9 7 +DEAL:float::OK +DEAL:float::5 17 +DEAL:float::OK +DEAL:float::320 121 +DEAL:float::OK -- 2.39.5