--- /dev/null
+New: LAPACKFullMatrix::add(const number, const Vector<number> &) performs
+a rank-1 update of a matrix.
+<br>
+(Denis Davydov, 2017/12/17)
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
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>
+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 (
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+ }
+}
--- /dev/null
+
+DEAL::size=17
+DEAL::Ok
+DEAL::size=35
+DEAL::Ok
+DEAL::size=391
+DEAL::Ok