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:

   This is the simplest of the available Sherman-Morrison-Woodbury kernels in QMCkl. It applies rank-1 updates one by one in the order that is given. It only checks if the denominator in the Sherman-Morrison formula is not too close to zero (and exit with an error if it does) during the application of an update. 

   #+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


* Woodbury 2x2
  
** ~qmckl_woodbury_2~
   :PROPERTIES:
   :Name:     qmckl_woodbury_2
   :CRetType: qmckl_exit_code
   :FRetType: qmckl_exit_code
   :END:

   This is the simplest of the available Sherman-Morrison-Woodbury kernels in QMCkl. It applies rank-1 updates one by one in the order that is given. It only checks if the denominator in the Sherman-Morrison formula is not too close to zero (and exit with an error if it does) during the application of an update. 

   #+NAME: qmckl_woodbury_2_args
   | qmckl_context | context                  | in    | Global state                                         |
   | uint64_t      | Dim                      | in    | Leading dimension of Slater_inv                      |
   | double        | Updates[2*Dim]           | in    | Array containing the updates                         |
   | uint64_t      | Updates_index[2]         | 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_woodbury_2_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_woodbury_2_c (
	  const qmckl_context context,
	  const uint64_t Dim,
	  const double* Updates,
	  const uint64_t* Updates_index,
	  double* Slater_inv ); 
    #+end_src

*** Source Fortran

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

*** Source C

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

qmckl_exit_code qmckl_woodbury_2_c(const qmckl_context context, 
                                const uint64_t Dim,
                                const double* Updates,
                                const uint64_t* Updates_index,
                                double * Slater_inv) {
/*
    C := S^{-1} * U,    dim x 2
    B := 1 + V * C,     2 x 2
    D := V * S^{-1},    2 x dim
*/
// #ifdef DEBUG //  Leave commented out since debugging information is not yet implemented in QMCkl.
//   std::cerr << "Called Woodbury 2x2 kernel" << std::endl;
// #endif

  const unsigned int row1 = (Updates_index[0] - 1);
  const unsigned int row2 = (Updates_index[1] - 1);

  // Compute C = S_inv * U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE
  // OF LAYOUT OF 'Updates' !!
  double C[2 * Dim];
  for (unsigned int i = 0; i < Dim; i++) {
    for (unsigned int j = 0; j < 2; j++) {
      C[i * 2 + j] = 0;
      for (unsigned int k = 0; k < Dim; k++) {
        C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k];
      }
    }
  }

  // Compute B = 1 + V * C
  const double B0 = C[row1 * 2] + 1;
  const double B1 = C[row1 * 2 + 1];
  const double B2 = C[row2 * 2];
  const double B3 = C[row2 * 2 + 1] + 1;

  // Check if determinant of inverted matrix is not zero
  double det = B0 * B3 - B1 * B2;
  double thresh = 0.0;
  qmckl_exit_code rc = qmckl_sherman_morrison_threshold_c(&thresh);
  if (fabs(det) < thresh) {
    return QMCKL_FAILURE;
  }

  // Compute B^{-1} with explicit formula for 2x2 inversion
  double Binv[4], idet = 1.0 / det;
  Binv[0] = idet * B3;
  Binv[1] = -1.0 * idet * B1;
  Binv[2] = -1.0 * idet * B2;
  Binv[3] = idet * B0;

  // Compute tmp = B^{-1} x (V.S^{-1})
  double tmp[2 * Dim];
  for (unsigned int i = 0; i < 2; i++) {
    for (unsigned int j = 0; j < Dim; j++) {
      tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1 * Dim + j];
      tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2 * Dim + j];
    }
  }

  // Compute (S + U V)^{-1} = S^{-1} - C x tmp
  for (unsigned int i = 0; i < Dim; i++) {
    for (unsigned int j = 0; j < Dim; j++) {
      Slater_inv[i * Dim + j] -= C[i * 2] * tmp[j];
      Slater_inv[i * Dim + j] -= C[i * 2 + 1] * tmp[Dim + j];
    }
  }

  return QMCKL_SUCCESS;
}

    #+end_src



*** Performance

** C interface                                                     :noexport:

   #+CALL: generate_c_interface(table=qmckl_woodbury_2_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_woodbury_2 &
       (context, Dim, 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
     real    (c_double ) , intent(in)          :: Updates(2*Dim)
     integer (c_int64_t) , intent(in)          :: Updates_index(2)
     real    (c_double ) , intent(inout)        :: Slater_inv(Dim*Dim)

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

   end function qmckl_woodbury_2
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_woodbury_2_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_woodbury_2 &
	 (context, Dim, 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
       real    (c_double ) , intent(in)          :: Updates(2*Dim)
       integer (c_int64_t) , intent(in)          :: Updates_index(2)
       real    (c_double ) , intent(inout)        :: Slater_inv(Dim*Dim)

     end function qmckl_woodbury_2
   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 Dim2 = 2;
const uint64_t Updates_index2[2] = {0, 0};
const double Updates2[4] = {0.0, 0.0, 0.0, 0.0};
double Slater_inv2[4] = {0.0, 0.0, 0.0, 0.0};

rc = qmckl_woodbury_2_c(context, Dim, 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