]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add LAPACK::Tmmult() to do a tripple product with a diagonal matrix
authorDenis Davydov <davydden@gmail.com>
Fri, 16 Feb 2018 17:03:54 +0000 (18:03 +0100)
committerDenis Davydov <davydden@gmail.com>
Fri, 16 Feb 2018 17:03:54 +0000 (18:03 +0100)
include/deal.II/lac/lapack_full_matrix.h
include/deal.II/lac/lapack_templates.h
source/lac/lapack_full_matrix.cc
tests/lapack/create_matrix.h
tests/lapack/full_matrix_31.cc [new file with mode: 0644]
tests/lapack/full_matrix_31.output [new file with mode: 0644]
tests/lapack/full_matrix_32.cc [new file with mode: 0644]
tests/lapack/full_matrix_32.output [new file with mode: 0644]

index a629d5d9a31a685a7aea9bd98a0ab3e1b43099f1..9608a08415578df2370e1b222183cb5947125293 100644 (file)
@@ -420,6 +420,25 @@ public:
                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>.
    *
@@ -475,6 +494,12 @@ public:
                 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.
    */
index 477787244453d91624b5a0894a713ecde70e10a5..45932cdb4a95dbc6cd7b0a02108f92d55548f1b0 100644 (file)
@@ -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 <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
index afb761b880d81ec8494be7a00b88e12787a8bf59..2073e5489dfe0d4732a763206c864696bda5bcad 100644 (file)
@@ -578,8 +578,8 @@ LAPACKFullMatrix<number>::mmult(LAPACKFullMatrix<number>       &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<number>::mmult(FullMatrix<number>             &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 <typename number>
 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()));
@@ -649,7 +716,7 @@ LAPACKFullMatrix<number>::Tmmult(FullMatrix<number>             &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<number>::mTmult(LAPACKFullMatrix<number>       &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<number>::mTmult(FullMatrix<number>             &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<number>::TmTmult(LAPACKFullMatrix<number>       &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<number>::TmTmult(FullMatrix<number>             &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()));
index a06bd006bbcceb7299d8cdd55ae02ba0aa7cd82d..6e47841d03aadcb7a6648ae7c7d9a7e08c58a2e6 100644 (file)
@@ -50,3 +50,12 @@ void create_random (FullMatrix &A)
     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>();
+}
diff --git a/tests/lapack/full_matrix_31.cc b/tests/lapack/full_matrix_31.cc
new file mode 100644 (file)
index 0000000..9f2917e
--- /dev/null
@@ -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 <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();
+
+}
diff --git a/tests/lapack/full_matrix_31.output b/tests/lapack/full_matrix_31.output
new file mode 100644 (file)
index 0000000..2bb351d
--- /dev/null
@@ -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 (file)
index 0000000..7c0cce2
--- /dev/null
@@ -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 <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();
+
+}
diff --git a/tests/lapack/full_matrix_32.output b/tests/lapack/full_matrix_32.output
new file mode 100644 (file)
index 0000000..c91d76f
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.