]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add rank-1 update of LAPACKFullMatrix
authorDenis Davydov <davydden@gmail.com>
Sun, 17 Dec 2017 16:06:15 +0000 (17:06 +0100)
committerDenis Davydov <davydden@gmail.com>
Tue, 19 Dec 2017 12:27:13 +0000 (13:27 +0100)
doc/news/changes/minor/20171217DenisDavydov [new file with mode: 0644]
include/deal.II/lac/lapack_full_matrix.h
include/deal.II/lac/lapack_templates.h
source/lac/lapack_full_matrix.cc
tests/lapack/full_matrix_27.cc [new file with mode: 0644]
tests/lapack/full_matrix_27.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20171217DenisDavydov b/doc/news/changes/minor/20171217DenisDavydov
new file mode 100644 (file)
index 0000000..b081117
--- /dev/null
@@ -0,0 +1,4 @@
+New: LAPACKFullMatrix::add(const number, const Vector<number> &) performs
+a rank-1 update of a matrix.
+<br>
+(Denis Davydov, 2017/12/17)
index 38ca6d253b1715a8ba5e2e0485088ce56b18bf31..f63a7dae6881bc734bdb0af1812a4fbba37ce5a9 100644 (file)
@@ -141,6 +141,13 @@ public:
   void add (const number                    a,
             const LAPACKFullMatrix<number> &A);
 
+  /**
+   * Perform a rank-1 update of a symmetric matrix
+   * $ A \leftarrow A + a \, \rm v \rm v^T $.
+   */
+  void add(const number a,
+           const Vector<number> &v);
+
   /**
    * Assignment from different matrix classes, performing the usual conversion
    * to the transposed format expected by LAPACK. This assignment operator
index c7cc243bc781396dee84e149c7c5f37e8c5cbf5b..dcda1559514a1cbbd63c3723005c90dcf223a02c 100644 (file)
@@ -297,11 +297,51 @@ extern "C"
                float *d, float *e, float *z,
                const int *ldz, float *work,
                int *info);
+// Rank-1 update for symmetric matrices
+  void dsyr_ (const char *uplo, const int *n,
+              const double *alpha, const double *x,
+              const int *incx,
+              double *A, const int *lda);
+  void ssyr_ (const char *uplo, const int *n,
+              const float *alpha, const float *x,
+              const int *incx,
+              float *A, const int *lda);
 
 }
 
 DEAL_II_NAMESPACE_OPEN
 
+/// Template wrapper for LAPACK functions dsyr and ssyr
+template <typename number>
+inline void
+syr(const char *, const int *,
+    const number *, const number *,
+    const int *,
+    number *, const int *)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+#ifdef DEAL_II_WITH_LAPACK
+inline void
+syr(const char *uplo, const int *n,
+    const double *alpha, const double *x,
+    const int *incx,
+    double *A, const int *lda)
+{
+  dsyr_(uplo,n,alpha,x,incx,A,lda);
+}
+inline void
+syr(const char *uplo, const int *n,
+    const float *alpha, const float *x,
+    const int *incx,
+    float *A, const int *lda)
+{
+  ssyr_(uplo,n,alpha,x,incx,A,lda);
+}
+#endif
+
+
 
 /// Template wrapper for LAPACK functions daxpy and saxpy
 template <typename number1, typename number2, typename number3>
index 0ac15627e7c605f23d4e9f4528dc727df0eebc5c..bed672572a4cf73e5207c901afcc07a2a146c093 100644 (file)
@@ -211,6 +211,39 @@ LAPACKFullMatrix<number>::add (const number              a,
 
 
 
+template <typename number>
+void
+LAPACKFullMatrix<number>::add(const number a,
+                              const Vector<number> &v)
+{
+  Assert(state == LAPACKSupport::matrix,
+         ExcState(state));
+  Assert(property == LAPACKSupport::symmetric,
+         ExcProperty(property));
+
+  Assert (n() == m(), ExcInternalError());
+  Assert (m() == v.size(), ExcDimensionMismatch(m(), v.size()));
+
+  AssertIsFinite(a);
+
+  const int N = this->n_rows();
+  const char uplo = LAPACKSupport::U;
+  const int lda = N;
+  const int incx=1;
+
+  syr(&uplo, &N, &a, v.begin(), &incx, this->values.begin(), &lda);
+
+  // FIXME: we should really only update upper or lower triangular parts
+  // of symmetric matrices and make sure the interface is consistent,
+  // for example operator(i,j) gives correct results regardless of storage.
+  for (unsigned int i = 0; i < N; ++i)
+    for (unsigned int j = 0; j < i; ++j)
+      (*this)(i,j) = (*this)(j,i);
+}
+
+
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::vmult (
diff --git a/tests/lapack/full_matrix_27.cc b/tests/lapack/full_matrix_27.cc
new file mode 100644 (file)
index 0000000..ce9d7aa
--- /dev/null
@@ -0,0 +1,80 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix::add() for rank1 update of a matrix
+
+#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>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+  // Lapack:
+  LAPACKFullMatrix<NumberType> A(size);
+  A.set_property(LAPACKSupport::symmetric);
+  Vector<NumberType> v(size);
+  for (unsigned int i = 0; i < size; ++i)
+    {
+      v(i) = random_value<NumberType>();
+      for (unsigned int j = i; j < size; ++j)
+        {
+          const NumberType val = random_value<NumberType>();
+          A(i,j) = val;
+          A(j,i) = val;
+        }
+    }
+
+
+  const NumberType a = random_value<NumberType>();
+
+  LAPACKFullMatrix<NumberType> B(size);
+  B = A;
+
+  B.add(a, v);
+
+  for (unsigned int i = 0; i < size; ++i)
+    for (unsigned int j = 0; j < size; ++j)
+      {
+        const NumberType diff = A(i,j) + a * v(i) * v(j) - B(i,j);
+        AssertThrow(std::abs(diff) < 1e-10 * std::abs(B(i,j)) ,
+                    ExcMessage("diff="+ std::to_string(diff)));
+      }
+
+  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<unsigned int> sizes = {{17,35,391}};
+  for (const auto &s : sizes)
+    {
+      deallog << "size=" << s << std::endl;
+      test<double>(s);
+    }
+}
diff --git a/tests/lapack/full_matrix_27.output b/tests/lapack/full_matrix_27.output
new file mode 100644 (file)
index 0000000..2f2ffd5
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::size=17
+DEAL::Ok
+DEAL::size=35
+DEAL::Ok
+DEAL::size=391
+DEAL::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.