qmckl/org/qmckl_sherman_morrison_woodbury.org

#+TITLE: Sherman-Morrison-Woodbury
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org

Low- and high-level functions that use the Sherman-Morrison and Woodbury matrix inversion formulas to update the
 inverse of a non-singualr matrix

* Headers
  #+begin_src elisp :noexport :results none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src

  #+begin_src c :comments link :tangle (eval c_test) :noweb yes
#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <cmath>
#ifndef THRESHOLD
#define THRESHOLD 1e-3
#endif

#include "qmckl.h"

int main() {
  qmckl_context context;
  context = qmckl_context_create();

  #+end_src
  
* Naïve   Sherman-Morrison
  
** ~qmckl_sherman_morrison~
   :PROPERTIES:
   :Name:     qmckl_sherman_morrison
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

   The Sherman-Morrison formula

   \begin{align}
	   S_k^{-1} &= (S_l + U_k)^-1 \\
			      &= S_l^{-1} - \frac{S_l^{-1}U_kS_l}{1+\underline{v}_k^tS_l^{-1}\underline{u}_k}
   \end{align}


   #+NAME: qmckl_sherman_morrison_args
   | qmckl_context | context                  | in    | Global state                                         |
   | uint64_t      | Dim                      | in    | Leading dimension of Slater_inv                      |
   | uint64_t      | N_updates                | in    | Number of rank-1 updates to be applied to Slater_inv |
   | double        | Updates[N_updates*Dim]   | in    | Array containing the updates                         |
   | uint64_t      | Updates_index[N_updates] | in    | Array containing the rank-1 updates                  |
   | double        | Slater_inv[Dim*Dim]      | inout | Array containing the inverse of a Slater-matrix      |

*** Requirements

      Add description of the input variables. (see for e.g. qmckl_distance.org)

*** C header

    #+CALL: generate_c_header(table=qmckl_sherman_morrison_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
// Sherman-Morrison-Woodbury break-down threshold
#ifndef THRESHOLD
#define THRESHOLD 1e-3
#endif
static double threshold();

    qmckl_exit_code qmckl_sherman_morrison_c (
	  const qmckl_context context,
	  const uint64_t Dim,
	  const uint64_t N_updates,
	  const double* Updates,
	  const uint64_t* Updates_index,
	  double* Slater_inv ); 
    #+end_src

*** Source Fortran

    #+begin_src f90 :tangle (eval f)
integer function qmckl_sherman_morrison_f(context, Slater_inv, Dim, N_updates,    &
         Updates, Updates_index) result(info) 
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  integer*8  , intent(in)  :: Dim, N_updates
  integer*8  , intent(in)  :: Updates_index(N_updates)
  real*8     , intent(in) :: Updates(N_updates*Dim)
  real*8     , intent(inout)  :: Slater_inv(Dim*Dim)
  !logical, external :: qmckl_sherman_morrison_f
  info = qmckl_sherman_morrison(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_f
    #+end_src

*** Source C

    #+begin_src c :tangle (eval c) :comments org
#include <stdbool.h>
#include "qmckl.h"

// Sherman-Morrison-Woodbury break-down threshold
static double threshold() {
  const double threshold = THRESHOLD;
#ifdef DEBUG
  std::cerr << "Break-down threshold set to: " << threshold << std::endl;
#endif
  return threshold;
}

qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, 
                                const uint64_t Dim,
                                const uint64_t N_updates,
                                const double* Updates,
                                const uint64_t* Updates_index,
                                double * Slater_inv) {
#ifdef DEBUG
  std::cerr << "Called qmckl_sherman_morrison with " << N_updates << " updates" << std::endl;
#endif

  double C[Dim];
  double D[Dim];

  unsigned int l = 0;
  // For each update
  while (l < N_updates) {
    // C = A^{-1} x U_l
    for (unsigned int i = 0; i < Dim; i++) {
      C[i] = 0;
      for (unsigned int j = 0; j < Dim; j++) {
        C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
      }
    }

    // Denominator
    double den = 1 + C[Updates_index[l] - 1];
    if (fabs(den) < threshold()) {
      return false;
    }
    double iden = 1 / den;

    // D = v^T x A^{-1}
    for (unsigned int j = 0; j < Dim; j++) {
      D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
    }

    // A^{-1} = A^{-1} - C x D / den
    for (unsigned int i = 0; i < Dim; i++) {
      for (unsigned int j = 0; j < Dim; j++) {
        double update = C[i] * D[j] * iden;
        Slater_inv[i * Dim + j] -= update;
      }
    }

    l += 1;
  }
  return QMCKL_SUCCESS;
}

    #+end_src



*** Performance

** C interface                                                     :noexport:

   #+CALL: generate_c_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval f) :comments org :exports none
   integer(c_int32_t) function qmckl_sherman_morrison &
       (context, Dim, N_updates, Updates, Updates_index, Slater_inv) &
       bind(C) result(info)

     use, intrinsic :: iso_c_binding
     implicit none

     integer (c_int64_t) , intent(in)  , value :: context
     integer (c_int64_t) , intent(in)  , value :: Dim
     integer (c_int64_t) , intent(in)  , value :: N_updates
     real    (c_double ) , intent(in)          :: Updates(N_updates*Dim)
     integer (c_int64_t) , intent(in)          :: Updates_index(N_updates)
     real    (c_double ) , intent(inout)        :: Slater_inv(Dim*Dim)

     integer(c_int32_t), external :: qmckl_sherman_morrison_c
     info = qmckl_sherman_morrison_c &
	    (context, Dim, N_updates, Updates, Updates_index, Slater_inv)

   end function qmckl_sherman_morrison
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

   #+RESULTS:
   #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
   interface
     integer(c_int32_t) function qmckl_sherman_morrison &
	 (context, Dim, N_updates, Updates, Updates_index, Slater_inv) &
	 bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

       integer (c_int64_t) , intent(in)  , value :: context
       integer (c_int64_t) , intent(in)  , value :: Dim
       integer (c_int64_t) , intent(in)  , value :: N_updates
       real    (c_double ) , intent(in)          :: Updates(N_updates*Dim)
       integer (c_int64_t) , intent(in)          :: Updates_index(N_updates)
       real    (c_double ) , intent(inout)        :: Slater_inv(Dim*Dim)

     end function qmckl_sherman_morrison
   end interface
   #+end_src

*** Test                                                           :noexport:

[TODO: FMJC] Write tests for the Sherman-Morrison part.

* End of files

   #+begin_src c :comments link :tangle (eval c_test)
  assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
  return 0;
}

   #+end_src


# -*- mode: org -*-
# vim: syntax=c