#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/utilities.h>
+
#include <iostream>
#include <iomanip>
+namespace
+{
+ template <typename number>
+ void cholesky_rank1(LAPACKFullMatrix<number> &A,
+ const number a,
+ const Vector<number> &v)
+ {
+ const unsigned int N = A.n();
+ Vector<number> z(v);
+ // Cholesky update / downdate, see
+ // 6.5.4 Cholesky Updating and Downdating, Golub 2013 Matrix computations
+ // Note that potrf() is called with LAPACKSupport::L , so the
+ // factorization is stored in lower triangular part.
+ if (a > 0.)
+ {
+ // simple update via a sequence of Givens rotations.
+ // Observe that
+ //
+ // | L^T |T | L^T |
+ // A <-- | | | | = L L^T + z z^T
+ // | z^T | | z^T |
+ //
+ // so we can get updated factor by doing a sequence of Givens
+ // rotations to make the matrix lower-triangular
+ // Also see LINPACK's dchud http://www.netlib.org/linpack/dchud.f
+ z *= std::sqrt(a);
+ for (unsigned int k = 0; k < N; ++k)
+ {
+ const std::array<number,3> csr = Utilities::LinearAlgebra::Givens_rotation(A(k,k),z(k));
+ A(k,k) = csr[2];
+ for (unsigned int i = k+1; i < N; ++i)
+ {
+ const number t = A(i,k);
+ A(i,k) = csr[0] * A(i,k) + csr[1] * z(i);
+ z(i) = -csr[1] * t + csr[0] * z(i);
+
+ //A(i,k) = ( t + csr[1] * z(i)) / csr[0];
+ //z(i) = - csr[1] * t + csr[0] * z(i) ;
+ }
+ }
+ }
+ else
+ {
+ AssertThrow(false, ExcNotImplemented());
+ }
+ }
+
+
+ template <typename number>
+ void cholesky_rank1(LAPACKFullMatrix<std::complex<number>> &A,
+ const std::complex<number> a,
+ const Vector<std::complex<number>> &v)
+ {
+ AssertThrow(false, ExcNotImplemented());
+ }
+}
+
+
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));
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);
+ if (state == LAPACKSupport::matrix)
+ {
+ 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);
+ }
+ else if (state == LAPACKSupport::cholesky)
+ {
+ cholesky_rank1(*this, a, v);
+ }
+ else
+ AssertThrow(false, ExcState(state));
}
if (info != 0)
std::cerr << "LAPACK error in geev" << std::endl;
- state = LAPACKSupport::State(eigenvalues | unusable);
+ state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable);
}
eigenvectors[i](j) = values_A[col_begin+j];
}
}
- state = LAPACKSupport::State(eigenvalues | unusable);
+ state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable);
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test LAPACKFullMatrix::add() for rank1 update of a Cholesky factorization
+
+/* MWE in Octave:
+A = pascal(4)
+A =
+
+ 1 1 1 1
+ 1 2 3 4
+ 1 3 6 10
+ 1 4 10 20
+
+x = [3 2 1 5]';
+
+R = chol(A)
+R =
+
+ 1 1 1 1
+ 0 1 2 3
+ 0 0 1 3
+ 0 0 0 1
+
+R1 = cholupdate(R,x)
+
+R1 =
+
+ 3.16228 2.21359 1.26491 5.05964
+ 0.00000 1.04881 2.09762 2.66970
+ 0.00000 0.00000 1.00000 3.00000
+ 0.00000 0.00000 0.00000 1.80907
+
+*/
+
+#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 = 4;
+ LAPACKFullMatrix<NumberType> A(size);
+ A.set_property(LAPACKSupport::symmetric);
+ Vector<NumberType> v(size);
+
+ A(0,0) = 1;
+ A(0,1) = 1;
+ A(0,2) = 1;
+ A(0,3) = 1;
+ A(1,0) = 1;
+ A(1,1) = 2;
+ A(1,2) = 3;
+ A(1,3) = 4;
+ A(2,0) = 1;
+ A(2,1) = 3;
+ A(2,2) = 6;
+ A(2,3) = 10;
+ A(3,0) = 1;
+ A(3,1) = 4;
+ A(3,2) =10;
+ A(3,3) = 20;
+
+ v(0) = 3;
+ v(1) = 2;
+ v(2) = 1;
+ v(3) = 5;
+
+ const NumberType a = 1.;
+
+ A.compute_cholesky_factorization();
+ deallog << "Cholesky:" << std::endl;
+ A.print_formatted(deallog.get_file_stream(),5);
+
+ A.add(a,v);
+
+ deallog << "Cholesky updated:" << std::endl;
+ A.print_formatted(deallog.get_file_stream(),5);
+}
+
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(5);
+ deallog.attach(logfile);
+
+ test<double>();
+}
projects / dealii.git / commitdiff