const LAPACKFullMatrix<number> &B,
const bool adding=false) const;
+ /**
+ * Matrix-matrix-multiplication using transpose of <tt>this</tt> and a
+ * diagonal vector @p V.
+ *
+ * If the <code>adding=false</code> 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<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const Vector<number> &V,
+ const bool adding=false) const;
+
/**
* Matrix-matrix-multiplication using transpose of <tt>B</tt>.
*
const LAPACKFullMatrix<number> &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<number> &V);
+
/**
* Compute the LU factorization of the matrix using LAPACK function Xgetrf.
*/
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);
DEAL_II_NAMESPACE_OPEN
+
+/// Template wrapper for LAPACK functions dsyrk and ssyrk
+template <typename number>
+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 <typename number>
inline void
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()));
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()));
void
LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number> &C,
const LAPACKFullMatrix<number> &B,
+ const Vector<number> &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<number> &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<static_cast<types::blas_int>(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 <typename number>
+void
+LAPACKFullMatrix<number>::scale_rows(const Vector<number> &V)
+{
+ LAPACKFullMatrix<number> &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 <typename number>
+void
+LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &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()));
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()));
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()));
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()));
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()));
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()));
for (unsigned int j = 0; j < A.n(); ++j)
A(i,j) = random_value<typename FullMatrix::value_type>();
}
+
+
+
+template <typename NumberType>
+void create_random (Vector<NumberType> &V)
+{
+ for (unsigned int i = 0; i < V.size(); ++i)
+ V(i) = random_value<NumberType>();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <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 m, const unsigned int n, const unsigned int k, const NumberType eps)
+{
+ deallog << m << " " << n << " " << k << " " << std::endl;
+ FullMatrix<NumberType> A(k, m), B(k, n), C(m, n), D(k,k);
+ LAPACKFullMatrix<NumberType> AL(k,m), BL(k, n), CL(m, n);
+ Vector<NumberType> DL(k);
+
+ create_random(AL);
+ create_random(BL);
+ create_random(DL);
+
+ A = AL;
+ B = BL;
+ D = NumberType();
+ for (unsigned int i=0; i<k; ++i)
+ D(i,i) = DL(i);
+
+ // non-symmetric
+ C.triple_product(D,A,B,true);
+ AL.Tmmult(CL,BL,DL);
+
+ for (unsigned int i=0; i<m; ++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 non-symmetric" << std::endl;
+
+ AL.Tmmult(CL,BL,DL,true);
+ for (unsigned int i=0; i<m; ++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 non-symmetric adding" << std::endl;
+
+ // symmetric (B==A)
+ if (m == n)
+ {
+ C = 0.;
+ C.triple_product(D,A,A,true,false,1.);
+
+ AL.Tmmult(CL,AL,DL);
+
+ for (unsigned int i=0; i<m; ++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 symmetric" << std::endl;
+
+ AL.Tmmult(CL,AL,DL,true);
+
+ for (unsigned int i=0; i<m; ++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 symmetric 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,3>> 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<double>(el[0],el[1],el[2],1e-13);
+ deallog.pop();
+
+ deallog.push("float");
+ for (auto el : sizes)
+ test<float>(el[0],el[1],el[2],1e-5);
+ deallog.pop();
+
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <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(k, n), D(k,k);
+ LAPACKFullMatrix<NumberType> AL(k,n);
+ Vector<NumberType> DL(k);
+
+ create_random(AL);
+ create_random(DL);
+
+ A = AL;
+ D = NumberType();
+ for (unsigned int i=0; i<k; ++i)
+ D(i,i) = DL(i);
+
+ // C = D*A
+ D.mmult(C,A);
+
+ // A = D * A
+ AL.scale_rows(DL);
+
+ for (unsigned int i=0; i<k; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ AssertThrow(std::abs(C(i,j)-AL(i,j)) < eps * std::abs(AL(i,j)),
+ ExcEl(i,j,C(i,j),AL(i,j),C(i,j)-AL(i,j))
+ );
+
+ deallog << "OK" << 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}, {5,17}, {320,121}
+ };
+
+ 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();
+
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff