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 :exports 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_threshold~
   :PROPERTIES:
   :Name:     qmckl_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_threshold_args
   | double | thresh | inout | Threshold |

*** Requirements

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

*** C header

    #+CALL: generate_c_header(table=qmckl_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_threshold_c (
	  double* const thresh ); 
    #+end_src

*** Source Fortran

*** Source C

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

// Sherman-Morrison-Woodbury break-down threshold
qmckl_exit_code qmckl_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

** ~qmckl_slagel_splitting~
   :PROPERTIES:
   :Name:     qmckl_slagel_splitting
   :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_slagel_splitting_args
   | 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 rank-1 updates                           |
   | uint64_t | Updates_index[N_updates]       | in    | Array containing positions of the rank-1 updates              |
   | double   | Slater_inv[Dim*Dim]            | inout | Array containing the inverse Slater-matrix                    |
   | double   | later_updates[Dim * N_updates] | inout | Array containing the split updates for later                  |
   | uint64_t | later_index[N_updates]         | inout | Array containing the positions of the split updates for later |
   | uint64_t | later                          | inout | Number of split updates for later                             |
  

*** Requirements

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

*** C header

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

*** Source Fortran

*** Source C

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

qmckl_exit_code qmckl_slagel_splitting_c(uint64_t Dim,
                                              uint64_t N_updates,
                                              const double *Updates, 
                                              const uint64_t *Updates_index,
                                              double *Slater_inv,
                                              double *later_updates, 
                                              uint64_t *later_index,
                                              uint64_t *later) {
// #ifdef DEBUG //  Leave commented out since debugging information is not yet implemented in QMCkl.
//   std::cerr << "Called slagel_splitting with " << N_updates << " updates" << std::endl;
// #endif

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

  uint64_t l = 0;
  // For each update
  while (l < N_updates) {
    // C = S^{-1} x U_l
    for (uint64_t i = 0; i < Dim; i++) {
      C[i] = 0;
      for (uint64_t 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_threshold_c(&thresh);
    if (fabs(den) < thresh) {

      // U_l = U_l / 2 (do the split)
      for (uint64_t i = 0; i < Dim; i++) {
        later_updates[*later * Dim + i] = Updates[l * Dim + i] / 2.0;
        C[i] /= 2.0;
      }
      later_index[*later] = Updates_index[l];
      (*later)++;

      den = 1 + C[Updates_index[l] - 1];
    }
    double iden = 1 / den;

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

    // S^{-1} = S^{-1} - C x D / den
    for (uint64_t i = 0; i < Dim; i++) {
      for (uint64_t 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_f_interface(table=qmckl_threshold_args,rettyp=get_value("FRetType"),fname=get_value("Name"))

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

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

  uint64_t l = 0;
  // For each update
  while (l < N_updates) {
    // C = A^{-1} x U_l
    for (uint64_t i = 0; i < Dim; i++) {
      C[i] = 0;
      for (uint64_t 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_threshold_c(&thresh);
    if (fabs(den) < thresh) {
      return QMCKL_FAILURE;
    }
    double iden = 1 / den;

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

    // A^{-1} = A^{-1} - C x D / den
    for (uint64_t i = 0; i < Dim; i++) {
      for (uint64_t 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};

// [TODO : FMJC ] add realistic tests

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


* Woodbury 2x2

This is the Woodbury 3x3 kernel.
  
** ~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 uint64_t row1 = (Updates_index[0] - 1);
  const uint64_t row2 = (Updates_index[1] - 1);

  // Compute C = S_inv * U  !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE
  // OF LAYOUT OF 'Updates' !!
  double C[2 * Dim];
  for (uint64_t i = 0; i < Dim; i++) {
    for (uint64_t j = 0; j < 2; j++) {
      C[i * 2 + j] = 0;
      for (uint64_t 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_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 (uint64_t i = 0; i < 2; i++) {
    for (uint64_t 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 (uint64_t i = 0; i < Dim; i++) {
    for (uint64_t 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 woodbury_Dim = 2;
const uint64_t woodbury_Updates_index[2] = {1, 1};
const double woodbury_Updates[4] = {1.0, 1.0, 1.0, 1.0};
double woodbury_Slater_inv[4] = {1.0, 1.0, 1.0, 1.0};

// [TODO : FMJC ] add realistic tests

rc = qmckl_woodbury_2_c(context, woodbury_Dim, woodbury_Updates, woodbury_Updates_index, woodbury_Slater_inv);
assert(rc == QMCKL_SUCCESS);
     #+end_src


* Woodbury 3x3

This is the Woodbury 3x3 kernel.
  
** ~qmckl_woodbury_3~
   :PROPERTIES:
   :Name:     qmckl_woodbury_3
   :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_3_args
   | qmckl_context | context             | in    | Global state                                    |
   | uint64_t      | Dim                 | in    | Leading dimension of Slater_inv                 |
   | double        | Updates[3*Dim]      | in    | Array containing the updates                    |
   | uint64_t      | Updates_index[3]    | 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_3_args,rettyp=get_value("CRetType"),fname=get_value("Name"))

    #+RESULTS:
    #+begin_src c :tangle (eval h_func) :comments org
    qmckl_exit_code qmckl_woodbury_3_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_3_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(3)
  real*8     , intent(in) :: Updates(3*Dim)
  real*8     , intent(inout)  :: Slater_inv(Dim*Dim)
  !logical, external :: qmckl_woodbury_3_f
  info = qmckl_woodbury_2(context, Dim, Updates, Updates_index, Slater_inv)
end function qmckl_woodbury_3_f
    #+end_src

*** Source C

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

qmckl_exit_code qmckl_woodbury_3_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 3
    B := 1 + V * C,     3 x 3
    D := V * S^{-1},    3 x dim
,*/
// #ifdef DEBUG //  Leave commented out since debugging information is not yet implemented in QMCkl.
//   std::cerr << "Called Woodbury 3x3 kernel" << std::endl;
// #endif

  const uint64_t row1 = (Updates_index[0] - 1);
  const uint64_t row2 = (Updates_index[1] - 1);
  const uint64_t row3 = (Updates_index[2] - 1);

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

  // Compute B = 1 + V.C
  const double B0 = C[row1 * 3] + 1;
  const double B1 = C[row1 * 3 + 1];
  const double B2 = C[row1 * 3 + 2];
  const double B3 = C[row2 * 3];
  const double B4 = C[row2 * 3 + 1] + 1;
  const double B5 = C[row2 * 3 + 2];
  const double B6 = C[row3 * 3];
  const double B7 = C[row3 * 3 + 1];
  const double B8 = C[row3 * 3 + 2] + 1;

  // Check if determinant of B is not too close to zero
  double det;
  det = B0 * (B4 * B8 - B5 * B7) - B1 * (B3 * B8 - B5 * B6) +
        B2 * (B3 * B7 - B4 * B6);
  double thresh = 0.0;
  qmckl_exit_code rc = qmckl_threshold_c(&thresh);
  if (fabs(det) < thresh) {
    return QMCKL_FAILURE;
  }

  // Compute B^{-1} with explicit formula for 3x3 inversion
  double Binv[9], idet = 1.0 / det;
  Binv[0] = (B4 * B8 - B7 * B5) * idet;
  Binv[1] = -(B1 * B8 - B7 * B2) * idet;
  Binv[2] = (B1 * B5 - B4 * B2) * idet;
  Binv[3] = -(B3 * B8 - B6 * B5) * idet;
  Binv[4] = (B0 * B8 - B6 * B2) * idet;
  Binv[5] = -(B0 * B5 - B3 * B2) * idet;
  Binv[6] = (B3 * B7 - B6 * B4) * idet;
  Binv[7] = -(B0 * B7 - B6 * B1) * idet;
  Binv[8] = (B0 * B4 - B3 * B1) * idet;

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

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

  return QMCKL_SUCCESS;
}

    #+end_src

*** Performance...

** C interface                                                     :noexport:

   #+CALL: generate_c_interface(table=qmckl_woodbury_3_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_3 &
       (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(3*Dim)
     integer (c_int64_t) , intent(in)          :: Updates_index(3)
     real    (c_double ) , intent(inout)        :: Slater_inv(Dim*Dim)

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

   end function qmckl_woodbury_3
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_woodbury_3_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_3 &
	 (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(3*Dim)
       integer (c_int64_t) , intent(in)          :: Updates_index(3)
       real    (c_double ) , intent(inout)        :: Slater_inv(Dim*Dim)

     end function qmckl_woodbury_3
   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 woodbury3_Dim = 3;
const uint64_t woodbury3_Updates_index[3] = {1, 1, 1};
const double woodbury3_Updates[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
double woodbury3_Slater_inv[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};

// [TODO : FMJC ] add realistic tests

rc = qmckl_woodbury_3_c(context, woodbury3_Dim, woodbury3_Updates, woodbury3_Updates_index, woodbury3_Slater_inv);
assert(rc == QMCKL_SUCCESS);
     #+end_src


* Sherman-Morrison with update splitting

This is like naïve Sherman-Morrising, but whenever a denominator is
found that is too close to zero the update is split in half. Then one
half is applied immediately and the other have is ket for later. When
all the updates have been processed, the list of split updates that
have been kept for later are processed. If again applying an update
results in a denominator that is too close to zero, it is split in
half again. One half is applied immediately and one half is kept for
later. The algorithm is done when no more updates have been kept for
later. This recursion will always end in a finite number of steps,
unless the last original update causes a singular Slater-matrix.
  
** ~qmckl_sherman_morrison_splitting~
   :PROPERTIES:
   :Name:     qmckl_sherman_morrison_splitting
   :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_splitting_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_splitting_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_splitting_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_splitting_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)
  info = qmckl_sherman_morrison_splitting(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_splitting_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_splitting_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 //  Leave commented out since debugging information is not yet implemented in QMCkl.
//   std::cerr << "Called qmckl_sherman_morrison_splitting with " << N_updates << " updates" << std::endl;
// #endif

  double later_updates[Dim * N_updates];
  uint64_t later_index[N_updates];
  uint64_t later = 0;

  qmckl_exit_code smss = qmckl_slagel_splitting_c(Dim, N_updates, Updates, Updates_index,
                            Slater_inv, later_updates, later_index, &later);
  
  if (later > 0) {
    qmckl_context context = qmckl_context_create();
    qmckl_exit_code sms = qmckl_sherman_morrison_splitting_c(context, Dim, later, 
                            later_updates, later_index, Slater_inv);
  }

  return QMCKL_SUCCESS;
}

    #+end_src

*** Performance...

** C interface                                                     :noexport:

   #+CALL: generate_c_interface(table=qmckl_sherman_morrison_splitting_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_splitting &
       (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_splitting_c
     info = qmckl_sherman_morrison_splitting_c &
	    (context, Dim, N_updates, Updates, Updates_index, Slater_inv)

   end function qmckl_sherman_morrison_splitting
   #+end_src

   #+CALL: generate_f_interface(table=qmckl_sherman_morrison_splitting_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_splitting &
	 (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_splitting
   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 splitting_Dim = 3;
const uint64_t splitting_N_updates = 3;
const uint64_t splitting_Updates_index[3] = {1, 1, 1};
const double splitting_Updates[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
double splitting_Slater_inv[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};

// [TODO : FMJC ] add realistic tests

rc = qmckl_sherman_morrison_splitting_c(context, splitting_Dim, splitting_N_updates, splitting_Updates, splitting_Updates_index, splitting_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