]> https://gitweb.dealii.org/ - dealii.git/commitdiff
implement rank-1 update for Cholesky factorization
authorDenis Davydov <davydden@gmail.com>
Mon, 18 Dec 2017 14:20:21 +0000 (15:20 +0100)
committerDenis Davydov <davydden@gmail.com>
Tue, 19 Dec 2017 12:27:14 +0000 (13:27 +0100)
source/lac/lapack_full_matrix.cc
tests/lapack/full_matrix_28.cc [new file with mode: 0644]
tests/lapack/full_matrix_28.output [new file with mode: 0644]

index 3e8246ce8a9a7f9b48491a224860c7b5928a50dc..c2007fd7cb08ba4c3b8151b36ded79b27f3eba3f 100644 (file)
@@ -22,6 +22,8 @@
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/block_vector.h>
 
+#include <deal.II/lac/utilities.h>
+
 #include <iostream>
 #include <iomanip>
 
@@ -211,13 +213,69 @@ LAPACKFullMatrix<number>::add (const number              a,
 
 
 
+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));
 
@@ -226,19 +284,28 @@ LAPACKFullMatrix<number>::add(const number a,
 
   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));
 }
 
 
@@ -1143,7 +1210,7 @@ LAPACKFullMatrix<number>::compute_eigenvalues(const bool right,
   if (info != 0)
     std::cerr << "LAPACK error in geev" << std::endl;
 
-  state = LAPACKSupport::State(eigenvalues | unusable);
+  state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable);
 }
 
 
@@ -1396,7 +1463,7 @@ LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric (
           eigenvectors[i](j) = values_A[col_begin+j];
         }
     }
-  state = LAPACKSupport::State(eigenvalues | unusable);
+  state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable);
 }
 
 
diff --git a/tests/lapack/full_matrix_28.cc b/tests/lapack/full_matrix_28.cc
new file mode 100644 (file)
index 0000000..4ddba3c
--- /dev/null
@@ -0,0 +1,110 @@
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
diff --git a/tests/lapack/full_matrix_28.output b/tests/lapack/full_matrix_28.output
new file mode 100644 (file)
index 0000000..6d36c7d
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::Cholesky:
+1.00000e+00  
+1.00000e+00  1.00000e+00  
+1.00000e+00  2.00000e+00  1.00000e+00  
+1.00000e+00  3.00000e+00  3.00000e+00  1.00000e+00  
+DEAL::Cholesky updated:
+3.16228e+00  
+2.21359e+00  1.04881e+00  
+1.26491e+00  2.09762e+00  1.00000e+00  
+5.05964e+00  2.66970e+00  3.00000e+00  1.80907e+00  

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.