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 <math.h>
#ifndef THRESHOLD
#define THRESHOLD 1e-3
#endif

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

  qmckl_exit_code rc;
  #+end_src

* Sherman-Morrison Helper Functions

Helper functions that are used by the Sherman-Morrison-Woodbury kernels. These functions should only be used in the context of these kernels.

** ~qmckl_sherman_morrison_threshold~
   :PROPERTIES:
   :Name:     qmckl_sherman_morrison_threshold
   :CRetType: double
   :FRetType: double precision
   :END:

This function is used to set the threshold value that is used in the kernels to determine if a matrix is invertable or not. In the Sherman-Morrison kernels this is determined by comparing the denominator in the Sherman-Morrison formula to the value set in threshold. If the value is smaller than the threshold value it means the matrix is not invertable. In the Woodbury kernels the threshold value is compared with the value of the determinant of the update matrix.

   #+NAME: qmckl_sherman_morrison_threshold_args
   | double | thresh | out | Threshold |

*** Requirements

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

*** C header

    #+CALL: generate_c_header(table=qmckl_sherman_morrison_threshold_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

    qmckl_exit_code qmckl_sherman_morrison_threshold_c (
	  double* const thresh ); 
    #+end_src


*** Source Fortran

    #+begin_src f90 :tangle (eval f)
integer function qmckl_sherman_morrison_threshold_f(thresh) result(info) 
  use qmckl
  implicit none
  real*8     , intent(inout)  :: thresh
  !logical, external :: qmckl_sherman_morrison_f
  info = qmckl_sherman_morrison_threshold(thresh)
end function qmckl_sherman_morrison_threshold_f
    #+end_src

*** Source C

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

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


*** Performance

** C interface                                                     :noexport:

   #+CALL: generate_c_interface(table=qmckl_sherman_morrison_threshold_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_threshold &
       (thresh) &
       bind(C) result(info)

     use, intrinsic :: iso_c_binding
     implicit none

     real    (c_double ) , intent(out)         :: thresh

     integer(c_int32_t), external :: qmckl_sherman_morrison_threshold_c
     info = qmckl_sherman_morrison_threshold_c &
	    (thresh)

   end function qmckl_sherman_morrison_threshold
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_sherman_morrison_threshold_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_threshold &
	 (thresh) &
	 bind(C)
       use, intrinsic :: iso_c_binding
       import
       implicit none

       real    (c_double ) , intent(out)         :: thresh

     end function qmckl_sherman_morrison_threshold
   end interface
   #+end_src

*** Test                                                           :noexport:

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



  
* 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
    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), value  :: 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"

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];

  double threshold = 0.0;
  qmckl_exit_code rc = qmckl_sherman_morrison_threshold_c(&threshold);

  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];
    double thresh = 0.0;
    qmckl_exit_code rc = qmckl_sherman_morrison_threshold_c(&thresh);
    if (fabs(den) < thresh) {
      return QMCKL_FAILURE;
    }
    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.


     #+begin_src c :tangle (eval c_test)
const uint64_t Dim = 2;
const uint64_t N_updates = 2;
const uint64_t Updates_index[2] = {0, 0};
const double Updates[4] = {0.0, 0.0, 0.0, 0.0};
double Slater_inv[4] = {0.0, 0.0, 0.0, 0.0};

rc = qmckl_sherman_morrison_c(context, Dim, N_updates, Updates, Updates_index, Slater_inv);
assert(rc == QMCKL_SUCCESS);
     #+end_src

* 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