qmckl/org/qmckl_sherman_morrison_woodbury.org

Sherman-Morrison-Woodbury

• Headers
• Naïve Sherman-Morrison
  • qmckl_sherman_morrison_naive
    • Introduction
    • API
    • Requirements
    • Pedagogical kernel source (in Fortran)
      • C interface to the pedagogical kernel (not directly exposed)
    • C headers (exposed in qmckl.h)
    • C sources
    • Fortran interfaces (exposed in qmckl_f.F90)
    • Performance
    • Tests
• Sherman-Morrison with update splitting
  • qmckl_sherman_morrison_splitting
    • Introduction
    • API
    • Requirements
    • Pedagogical kernel source (in Fortran)
      • C interface to the pedagogical kernel (not directly exposed)
    • C headers (exposed in qmckl.h)
    • C source
    • Fortran interfaces (exposed in qmckl_f.F90)
    • Performance…
    • Test
• Slagel Splitting
  • qmckl_slagel_splitting
    • Introduction
    • API
    • Requirements
    • Pedagogical kernel source (in Fortran)
      • C interface to the pedagogical kernel (not directly exposed)
    • C headers (exposed in qmckl.h)
    • C sources
    • Fortran interfaces (exposed in qmckl_f.F90)
    • Performance
• End of files

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

Headers

#include "qmckl.h"
#include "assert.h"
#ifdef HAVE_CONFIG_H
  #include "config.h"
#endif
#include <math.h>

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

This is the range that determines the how many high performance kernel instantces will be generated, using the C-function templates defined in the sections below. If the name of the C-function template is called qmckl_kernel_{Dim}, then range(K, L+1) will results in kernel instances from qmckl_kernel_K to qmckl_kernel_L. #+NAME:kernel_generator_range

Naïve Sherman-Morrison

qmckl_sherman_morrison_naive

Introduction

This is the simplest of the available Sherman-Morrison-Woodbury kernels. 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 when an update is evaluated. It will exit with an error code of the denominator is too close to zero.

#+TODO Change the math notation so that the update vectors appear as row in the math so that it is consistent with the representation in C (memory)

The formula for any update $u_j$ (index $j$ is suppresed for clarity) that is applied is \[ (S + uv^T)^{-1} = S^{-1} - \frac{S^{-1} uv^T S^{-1}}{1 + v^T S^{-1} u} \]

where $S$ is the Slater-matrix, $u$ and $v^T$ are the column and row vectors containing the updates, $S^{-1}$ is the inverse of the Slater-matrix.

Even though the Slater-matrix $S$ with all updates applied at once is invertable, during the course of applying updates to the inverse Slater-matrix $S^{-1}$ one-by-one it can happen that one of the intermediate inverse matrices $S^{-1}$ becomes singular. Therefore a global threshold value $\epsilon$ is defined that is used to evaluate each individual update $u_j$ when it is applied.

This value sets the lower bound for which the denominator $1+v_j^TS^{-1}u_j$ is considered to be too small and will most probably result in a singular matrix $S$, or at least in an inverse of $S$ of very poor numerical quality. Therefore, when $1+v_j^TS^{-1}u_j \geq \epsilon$, the update is applied as usual and the kernel exits with return code \texttt{QMCKL_SUCCESS}. If $1+v_j^TS^{-1}u_j \leq \epsilon$ the update is rejected and the kernel exits with return code \texttt{QMCKL_FAILURE}.

If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting from applying the updates to the original matrix.

API

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
LDS	uint64_t	in	Leading dimension of Slater_inv
Dim	uint64_t	in	Dimension of Slater_inv
N_updates	uint64_t	in	Number of rank-1 updates to be applied to Slater_inv
Updates	double[N_updates*LDS]	in	Array containing the updates
Updates_index	uint64_t[N_updates]	in	Array containing the rank-1 updates
breakdown	double	in	Break-down parameter on which to fail or not
Slater_inv	double[Dim*LDS]	inout	Array containing the inverse of a Slater-matrix
determinant	double	inout	Determinant of the Slater-matrix

Requirements

• context is not QMCKL_NULL_CONTEXT
• LDS >= 2
• Dim >= 2
• N_updates >= 1
• Updates is allocated with $N_updates \times Dim$ elements
• Updates_index is allocated with $N_updates$ elements
• breakdown is a small number such that $0 < breakdown << 1$
• Slater_inv is allocated with $Dim \times Dim$ elements
• determinant > 0

Pedagogical kernel source (in Fortran)

The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore not be used in real workloads.

integer function qmckl_sherman_morrison_naive_doc_f(context, &
    lds, dim, &
    nupdates, &
    upds, &
    updates_index, &
    breakdown, &
    s_inv, &
    determinant) result(info)

  use qmckl
  implicit none
  integer*8 , intent(in)    :: context
  integer*8 , intent(in)    :: lds, dim
  integer*8 , intent(in)    :: nupdates
  integer*8 , intent(in)    :: updates_index(nupdates)
  real*8    , intent(in)    :: upds(nupdates * lds)
  real*8    , intent(in)    :: breakdown
  real*8    , intent(inout) :: s_inv(dim * lds)
  real*8    , intent(inout) :: determinant

  real*8 , dimension(lds, nupdates) :: Updates
  real*8 , dimension(dim, lds)      :: Inverse
  real*8 , dimension(dim)           :: C
  real*8 , dimension(lds)           :: D
  real*8                            :: denominator, idenominator, update
  integer*8                         :: i, j, l, row

  info = QMCKL_FAILURE

  if (context == QMCKL_NULL_CONTEXT) then
    info = QMCKL_INVALID_CONTEXT
    return
  endif

  ! Convert 'upds' and 's_inv' into the more easily readable Fortran
  ! matrices 'Updates' and 'Inverse'.
  call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)

  l = 1;
  ! For each update do...
  do while (l < nupdates + 1)

    ! Compute C = S^{-1}U(l)
    do i = 1, dim
      C(i) = 0
      do j = 1, dim
        C(i) = C(i) + Inverse(i, j) * Updates(j, l)
      end do
    end do

    ! Compute denominator = 1 + V(l)^TC
    row = updates_index(l)
    denominator = 1 + C(row)

    ! Return early if denominator is too small
    if (abs(denominator) < breakdown) return
    idenominator = 1 / denominator

    ! Update det(S)
    determinant = determinant * denominator

    ! selecting column: v_l^T * S_inv
    D = Inverse(row, :)

    ! A^{-1} = A^{-1} - C x D / denominator
    do i = 1, dim
      do j = 1, dim
        update = C(i) * D(j) * idenominator
        Inverse(i, j) = Inverse(i, j) - update
      end do
    end do

    l = l + 1
  end do

  ! Copy updated inverse back to s_inv
  call copy_back(Inverse, s_inv, lds, dim)

  info = QMCKL_SUCCESS

end function qmckl_sherman_morrison_naive_doc_f

C interface to the pedagogical kernel (not directly exposed)

The following Fortran function qmckl_sherman_morrison_naive_doc makes sure that the pedagogical kernel qmckl_sherman_morrison_naive_doc_f, written in Fortran, can be called from C using the ISO_C_BINDING. The Fortran function qmckl_sherman_morrison_naive_doc will be exposed in the header file 'qmckl.h' for C users and in the module file 'qmckl_f.F90' for Fortran users.

C headers (exposed in qmckl.h)

qmckl_exit_code qmckl_sherman_morrison_naive (
      const qmckl_context context,
      const uint64_t LDS,
      const uint64_t Dim,
      const uint64_t N_updates,
      const double* Updates,
      const uint64_t* Updates_index,
      const double breakdown,
      double* Slater_inv,
      double* determinant );

qmckl_exit_code qmckl_sherman_morrison_naive_hpc (
      const qmckl_context context,
      const uint64_t LDS,
      const uint64_t Dim,
      const uint64_t N_updates,
      const double* Updates,
      const uint64_t* Updates_index,
      const double breakdown,
      double* Slater_inv,
      double* determinant );

qmckl_exit_code qmckl_sherman_morrison_naive_doc (
      const qmckl_context context,
      const uint64_t LDS,
      const uint64_t Dim,
      const uint64_t N_updates,
      const double* Updates,
      const uint64_t* Updates_index,
      const double breakdown,
      double* Slater_inv,
      double* determinant );

C sources

Common includes and macros used by all the Sherman-Morrison-Woodbury kernels.

#include <stdbool.h>
#include <math.h>
#include "qmckl.h"
#include "config.h"

// Order important because
// __GNUC__ also set in ICC, ICX and CLANG
// __clang__ also set in ICX
#if defined(__INTEL_COMPILER)
    #define IVDEP _Pragma("ivdep")
    #define ALIGNED _Pragma("vector aligned") 
#elif defined(__INTEL_LLVM_COMPILER)
    #define IVDEP _Pragma("ivdep")
    #define ALIGNED _Pragma("vector aligned") 
#elif defined(__clang__)
    #define IVDEP _Pragma("clang loop vectorize(enable)")
    #define ALIGNED
#elif defined(__GNUC__)
    #define IVDEP _Pragma("GCC ivdep")
    #define ALIGNED
#endif

qmckl_sherman_morrison_naive_hpc is a high performance variation of qmckl_sherman_morrison_naive written in C. It is used in cases when Dim is smaller than the leading dimension LDS, irrespective of whetether LDS includes zero padding to benefit from SIMD instructions or not. Cases like this include situations where one wants to apply updates to a square submatrix of the full matrix. It takes advantage of memory aligned data and assumes no data dependencies inside the loops. The loops are fully vectorised whenever Dim is an integer multiple of SIMD_LEGTH.

qmckl_exit_code qmckl_sherman_morrison_naive_hpc(
    const qmckl_context context,
    const uint64_t LDS,
    const uint64_t Dim,
    const uint64_t N_updates,
    const double* __restrict Updates,
    const uint64_t* __restrict Updates_index,
    const double breakdown,
    double* __restrict Slater_inv,
    double* __restrict determinant) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
      return qmckl_failwith( context,
          QMCKL_NULL_CONTEXT,
          "qmckl_sherman_morrison_naive_hpc",
          NULL);
  }

  double __attribute__((aligned(8))) C[Dim];
  double __attribute__((aligned(8))) D[LDS];

  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.0f;
      IVDEP
      ALIGNED
      for (uint64_t j = 0; j < Dim; j++) {
        C[i] += Slater_inv[i * LDS + j] * Updates[l * LDS + j];
      }
    }

    // Denominator: v_l^T * C
    const int cui = Updates_index[l] - 1;
    double den = 1.0f + C[cui];

    if (fabs(den) < breakdown)
      return QMCKL_FAILURE;

    double iden = 1.0f / den;

    // Update det(A)
    if (determinant)
      *determinant *= den;

    // selecting column: v_l^T * S_inv
    IVDEP
    ALIGNED
    for (uint64_t j = 0; j < Dim; j++) {
      D[j] = Slater_inv[cui * LDS + j];
    }

    // A^{-1} = A^{-1} - C x D / den
    for (uint64_t i = 0; i < Dim; i++) {
      IVDEP
      ALIGNED
      for (uint64_t j = 0; j < Dim; j++) {
        const double update = C[i] * D[j] * iden;
        Slater_inv[i * LDS + j] -= update;
      }
    }
    l += 1;
  }
  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_sherman_morrison_naive_{Dim} is a C function-template that is used to genereate instances of C fucntions based on the range given above. The advantage of this method is that for each of these instances all the dimensions and loop-bounds are known at compile time, allowing the compiler to optimize more aggressively. #+NAME:naive_template_code

  static inline qmckl_exit_code qmckl_sherman_morrison_naive_{Dim}(
      const qmckl_context context,
      const uint64_t N_updates,
      const double* __restrict Updates,
      const uint64_t* __restrict Updates_index,
      const double breakdown,
      double* __restrict Slater_inv,
      double* __restrict determinant) {

    if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
      return qmckl_failwith(context,
          QMCKL_NULL_CONTEXT,
          "qmckl_sherman_morrison_naive_{Dim}",
          NULL);
    }

    #define D{Dim}_P ((1+({Dim}-1)/SIMD_LENGTH)*SIMD_LENGTH)

    double __attribute__((aligned(8))) C[{Dim}];
    double __attribute__((aligned(8))) D[D{Dim}_P];

    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;
        IVDEP
        ALIGNED
        for (uint64_t j = 0; j < D{Dim}_P; j++) {
          C[i] += Slater_inv[i * D{Dim}_P + j] * Updates[l * D{Dim}_P + j];
        }
      }

      // Denominator
      const int cui = Updates_index[l] - 1;
      double den = 1.0f + C[cui];

      if (fabs(den) < breakdown) {
        return QMCKL_FAILURE;
      }
      double iden = 1.0f / den;

      // Update det(A)
      if (determinant)
        *determinant *= den;

      // selecting column: D = v_l^T * S_inv
      IVDEP
      ALIGNED
      for (uint64_t j = 0; j < D{Dim}_P; j++) {
        D[j] = Slater_inv[cui * D{Dim}_P + j];
      }

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

      l += 1;
    }

    return QMCKL_SUCCESS;
  }

This is the kernel generator written in Python. It uses the kernel generator range and templates defined above to generate the C kernel instances. #+NAME:naive_kernel_generator

text="""
<<naive_template_code>>
"""
result = []
for Dim in <<kernel_generator_range>>:
    Dim=str(Dim)
    result.append(text.replace("{Dim}",Dim) )

return '\n'.join(result)

Python script that generated C switch cases that call individual kernel instances. #+NAME:naive_switch-case_generator

text="""
case {Dim}:  
  return qmckl_sherman_morrison_naive_{Dim}(context,
    N_updates,
    Updates,
    Updates_index,
    breakdown,
    Slater_inv,
    determinant);
"""
result = []
for Dim in <<kernel_generator_range>>:
    Dim=str(Dim)
    result.append(text.replace("{Dim}",Dim) )

return '\n'.join(result)

<<naive_kernel_generator()>>

qmckl_sherman_morrison_naive is a generic function that contains decision making logic that calls the proper kernel based on the used library configuration (--enable-doc and --enable-hpc) and the passed array dimensions LDS and Dim.

qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context,
                                const uint64_t LDS,
                                const uint64_t Dim,
                                const uint64_t N_updates,
                                const double* Updates,
                                const uint64_t* Updates_index,
                                const double breakdown,
                                double* Slater_inv,
                                double* determinant) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith(context,
                          QMCKL_NULL_CONTEXT,
                          "qmckl_sherman_morrison_naive",
                          NULL);
  }

  #ifdef HAVE_HPC
  if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases
    switch (Dim) {
      <<naive_switch-case_generator()>>
    }
  }
  else { // Updating smaller sub-matrix
    return qmckl_sherman_morrison_naive_hpc(context,
                                      LDS,
                                      Dim,
    	                              N_updates,
    	                              Updates,
    	                              Updates_index,
    	                              breakdown,
    	                              Slater_inv,
    	                              determinant);
  }
  #else
  return qmckl_sherman_morrison_naive_doc(context,
                                      LDS,
                                      Dim,
    	                              N_updates,
                                      Updates,
    	                              Updates_index,
    	                              breakdown,
    	                              Slater_inv,
    	                              determinant);
  #endif
  
  return QMCKL_FAILURE;
}

Fortran interfaces (exposed in qmckl_f.F90)

Performance

This function performs best when there is only 1 rank-1 update in the update cycle. It is not useful to use Sherman-Morrison with update splitting for these cycles since splitting can never resolve a situation where applying the update causes singular behaviour.

Tests

The tests for the kernels are executed on datasets that are extracted from a run of QMC=Chem on Benzene (21 spin-up/21 spin down electrons) using 329 unique alpha determinants. The tests are run such that the kernels reject the computed inverse whenever the computed intermediate determinants or denominators are smaller than 1e-3. This is the default value in QMC=Chem. The tests will return QMCKL_SUCCESS whenever all the elements of the final matrix $R=S.S^-1 - 1$ are smaller than the given tolerance value of 1e-3, and will return QMCKL_FAILURE if the values are larger than this tolerance value.

const uint64_t Dim = 21;
const uint64_t LDS = (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH;
const double breakdown = 1e-3;
const double tolerance = 1e-3;
double res[441];

#include "sm_test.h"

assert(Updates1 != NULL);
assert(Updates_index1 != NULL);
assert(Slater_inv1 != NULL);

// original determinant of Slater1 (before applying updates)
double det = 3.407025646103221e-10;
rc = qmckl_sherman_morrison_naive(context,
                            LDS,
                            Dim,
                            N_updates1,
                            Updates1,
                            Updates_index1,
                            breakdown,
                            Slater_inv1,
                            &det);

// Check that the determinant is updated properly
assert(fabs(det + 4.120398385068217e-10) < 1e-15);

for (unsigned int i = 0; i < Dim; i++) {
  for (unsigned int j = 0; j < Dim; j++) {
    res[i * Dim + j] = 0;
    for (unsigned int k = 0; k < Dim; k++) {
      res[i * Dim + j] += Slater1[i * Dim + k] * Slater_inv1[k * LDS + j];
    }
  }
}
rc = QMCKL_SUCCESS;
for (unsigned int i = 0; i < Dim; i++) {
  for (unsigned int j = 0; j < Dim; j++) {
    if (i == j && fabs(res[i * Dim + j] - 1) > tolerance) {
      rc = QMCKL_FAILURE;
    }
    if (i != j && fabs(res[i * Dim + j]) > tolerance) {
      rc = QMCKL_FAILURE;
    }
  }
}
assert(rc == QMCKL_SUCCESS);

Sherman-Morrison with update splitting

qmckl_sherman_morrison_splitting

Introduction

This is a variation on the 'Naive' Sherman-Morrison kernel. Whenever the denominator $1+v_j^T S^{-1} u_j$ in the Sherman-Morrison formula is deemed to be too close to zero, the update $u_j$ is split in half: $u_j \rightarrow \frac{1}{2} u_j$. One half is applied immediately –necessarily increasing the value of the denominator because of the split– while the other halve is put in a queue that will be applied when all the remaining updates have been treated.

The kernel is executed recursively until the queue is eiter empty and all updates are applied successfully, or the size of the queue equals the number of initial updates. In the last case the Slater-matrix that would have resulted from applying the updates is singular and therefore the kernel exits with an exit code.

If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting from applying the updates to the original matrix.

API

Variable	Type	In/Out	Description
context	qmckl_context	in	Global state
LDS	uint64_t	in	Leading dimension of Slater_inv
Dim	uint64_t	in	Dimension of Slater_inv
N_updates	uint64_t	in	Number of rank-1 updates to be applied to Slater_inv
Updates	double[N_updates*LDS]	in	Array containing the updates
Updates_index	uint64_t[N_updates]	in	Array containing the rank-1 updates
breakdown	double	in	Break-down parameter on which to fail or not
Slater_inv	double[Dim*LDS]	inout	Array containing the inverse of a Slater-matrix
determinant	double	inout	Determinant of the Slater-matrix

Requirements

• context is not QMCKL_NULL_CONTEXT
• LDS >= 2
• Dim >= 2
• N_updates >= 1
• Updates is allocated with $N_updates \times Dim$ elements
• Updates_index is allocated with $N_updates$ elements
• breakdown is a small number such that $0 < breakdown << 1$
• Slater_inv is allocated with $Dim \times Dim$ elements

Pedagogical kernel source (in Fortran)

The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore not be used in real workloads.

integer function qmckl_sherman_morrison_splitting_doc_f(context, &
    lds, dim, &
    nupdates, &
    upds, &
    updates_index, &
    breakdown, &
    s_inv, &
    determinant) result(info)

  use qmckl
  implicit none
  integer*8 , intent(in)    :: context
  integer*8 , intent(in)    :: lds, dim
  integer*8 , intent(in)    :: nupdates
  integer*8 , intent(in)    :: updates_index(nupdates)
  real*8    , intent(in)    :: upds(nupdates * lds)
  real*8    , intent(in)    :: breakdown
  real*8    , intent(inout) :: s_inv(dim * lds)
  real*8    , intent(inout) :: determinant

  real*8 , dimension(lds, nupdates) :: Updates
  real*8 , dimension(dim, lds)      :: Inverse

  info = QMCKL_FAILURE

  ! Convert 'upds' and 's_inv' into the more easily readable Fortran
  ! matrices 'Updates' and 'Inverse'.
  call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)

  ! YET TO BE IMPLEMENTED

  ! Copy updated inverse back to s_inv
  call copy_back(Inverse, s_inv, lds, dim)

  info = QMCKL_SUCCESS

end function qmckl_sherman_morrison_splitting_doc_f

C interface to the pedagogical kernel (not directly exposed)

The following Fortran function qmckl_slagel_splitting_doc makes sure that the pedagogical kernel qmckl_slagel_splitting_doc_f, written in Fortran, can be called from C using the ISO_C_BINDING. The Fortran function qmckl_slagel_splitting_doc will be exposed in the header file 'qmckl.h' for C users and in the module file 'qmckl_f.F90' for Fortran users.

C headers (exposed in qmckl.h)

qmckl_exit_code qmckl_sherman_morrison_splitting (
      const qmckl_context context,
      const uint64_t LDS,
      const uint64_t Dim,
      const uint64_t N_updates,
      const double* Updates,
      const uint64_t* Updates_index,
      const double breakdown,
      double* Slater_inv,
      double* determinant );

qmckl_exit_code qmckl_sherman_morrison_splitting_hpc (
      const qmckl_context context,
      const uint64_t LDS,
      const uint64_t Dim,
      const uint64_t N_updates,
      const double* Updates,
      const uint64_t* Updates_index,
      const double breakdown,
      double* Slater_inv,
      double* determinant );

qmckl_exit_code qmckl_sherman_morrison_splitting_doc (
      const qmckl_context context,
      const uint64_t LDS,
      const uint64_t Dim,
      const uint64_t N_updates,
      const double* Updates,
      const uint64_t* Updates_index,
      const double breakdown,
      double* Slater_inv,
      double* determinant );

C source

qmckl_exit_code qmckl_sherman_morrison_splitting_hpc(const qmckl_context context,
                                const uint64_t LDS,
                                const uint64_t Dim,
                                const uint64_t N_updates,
                                const double* Updates,
                                const uint64_t* Updates_index,
                                const double breakdown,
                                double* Slater_inv,
                                double* determinant) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith(context,
                          QMCKL_NULL_CONTEXT,
                          "qmckl_sherman_morrison_splitting",
                          NULL);
  }

  double __attribute__((aligned(8))) later_updates[LDS * N_updates];
  uint64_t later_index[N_updates];
  uint64_t later = 0;

  qmckl_exit_code rc = qmckl_slagel_splitting(
    LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv,
    later_updates, later_index, &later, determinant);
  if (rc != QMCKL_SUCCESS) return QMCKL_FAILURE;

  if (later > 0) {
    qmckl_exit_code rc = qmckl_sherman_morrison_splitting(
      context, LDS, Dim, later, later_updates, later_index, breakdown,
      Slater_inv, determinant);
    if (rc != QMCKL_SUCCESS) return QMCKL_FAILURE;
  }

  return QMCKL_SUCCESS;
}

qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
                                const uint64_t LDS,
                                const uint64_t Dim,
                                const uint64_t N_updates,
                                const double* Updates,
                                const uint64_t* Updates_index,
                                const double breakdown,
                                double* Slater_inv,
                                double* determinant) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith(context,
                          QMCKL_NULL_CONTEXT,
                          "qmckl_sherman_morrison_splitting",
                          NULL);
  }

  #ifdef HAS_HPC
  return qmckl_sherman_morrison_splitting_hpc(context,
                                      LDS,
                                      Dim,
    	                              N_updates,
    	                              Updates,
    	                              Updates_index,
    	                              breakdown,
    	                              Slater_inv,
    	                              determinant);
  #else
  // return qmckl_sherman_morrison_splitting_doc(context,
  //                                  LDS,
  //                                  Dim,
  //   	                              N_updates,
  //   	                              Updates,
  //   	                              Updates_index,
  //   	                              breakdown,
  //   	                              Slater_inv,
  //   	                              determinant);
  return qmckl_sherman_morrison_splitting_hpc(context,
                                      LDS,
                                      Dim,
    	                              N_updates,
    	                              Updates,
    	                              Updates_index,
    	                              breakdown,
    	                              Slater_inv,
    	                              determinant);
  #endif
  
  return QMCKL_SUCCESS;
}

Fortran interfaces (exposed in qmckl_f.F90)

Performance…

This kernel performs best when there are 2 or more rank-1 update cycles and fail-rate is high.

Test

assert(Updates3 != NULL);
assert(Updates_index3 != NULL);
assert(Slater_inv3_2 != NULL);
det = -1.23743195512859e-09;
rc = qmckl_sherman_morrison_splitting(context, LDS, Dim, N_updates3, Updates3, Updates_index3, breakdown, Slater_inv3_2, &det);
assert(fabs(det - 1.602708950725074e-10) < 1e-15);
for (unsigned int i = 0; i < Dim; i++) {
  for (unsigned int j = 0; j < Dim; j++) {
    res[i * Dim + j] = 0;
    for (unsigned int k = 0; k < Dim; k++) {
      res[i * Dim + j] += Slater3[i * Dim + k] * Slater_inv3_2[k * LDS + j];
    }
  }
}
rc = QMCKL_SUCCESS;
for (unsigned int i = 0; i < Dim; i++) {
  for (unsigned int j = 0; j < Dim; j++) {
    if (i == j && fabs(res[i * Dim + j] - 1) > tolerance) {
      rc = QMCKL_FAILURE;
    }
    if (i != j && fabs(res[i * Dim + j]) > tolerance) {
      rc = QMCKL_FAILURE;
    }
  }
}
assert(rc == QMCKL_SUCCESS);

Slagel Splitting

qmckl_slagel_splitting

Introduction

qmckl_slagel_splitting is the non-recursive, inner part of the 'Sherman-Morrison with update splitting'-kernel. It is used internally to apply a collection of $N$ rank-1 updates to the inverse Slater-matrix $S^{-1}$ and splitting an update in two equal pieces if necessary. In case of a split, it applies the first half of the update, while putting the second half in a waiting queue to be applied at the end.

Therefore, when $1+v_j^TS^{-1}u_j \geq \epsilon$, the update is applied as usual. Otherwise, $u_j$ will be redefined as $\frac{1}{2}u_j$. One half is applied immediately, the other half will be applied at the end of the algorithm, using vectors $u_{j'}=\frac{1}{2}u_j$ and $v_{j'}^T=v_{j}^T$, which are stored in the array \texttt{later_updates}.

If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting from applying the updates to the original matrix.

API

Variable	Type	In/Out	Description
LDS	uint64_t	in	Leading dimension of Slater_inv
Dim	uint64_t	in	Dimension of Slater_inv
N_updates	uint64_t	in	Number of rank-1 updates to be applied to Slater_inv
Updates	double[N_updates*LDS]	in	Array containing the rank-1 updates
Updates_index	uint64_t[N_updates]	in	Array containing positions of the rank-1 updates
breakdown	double	in	Break-down parameter on which to fail or not
Slater_inv	double[Dim*LDS]	inout	Array containing the inverse Slater-matrix
later_updates	double[N_updates*LDS]	inout	Array containing the split updates for later
later_index	uint64_t[N_updates]	inout	Array containing the positions of the split updates for later
later	uint64_t	inout	Number of split updates for later
determinant	double	inout	Determinant of the Slater-matrix

Requirements

• LDS >= 2
• Dim >= 2
• N_updates >= 1
• Updates is allocated with $N_updates \times Dim$ elements
• Updates_index is allocated with $N_updates$ elements
• breakdown is a small number such that $0 < breakdown << 1$
• Slater_inv is allocated with $Dim \times Dim$ elements
• later_updates is allocated with $later \times Dim$ elements
• later_index is allocated with $N_updates$ elements
• later >= 0

Pedagogical kernel source (in Fortran)

The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore not be used in real workloads.

integer function qmckl_slagel_splitting_doc_f( &
    lds, dim, &
    nupdates, &
    upds, &
    updates_index, &
    breakdown, &
    s_inv, &
    later_updates, &
    later_index, &
    later, &
    determinant) result(info)

  use qmckl
  implicit none
  integer*8 , intent(in)    :: lds, dim
  integer*8 , intent(in)    :: nupdates
  integer*8 , intent(in)    :: updates_index(nupdates)
  real*8    , intent(in)    :: upds(nupdates * lds)
  real*8    , intent(in)    :: breakdown
  real*8    , intent(inout) :: s_inv(dim * lds)
  real*8    , intent(inout) :: later_updates(nupdates * lds)
  integer*8 , intent(inout) :: later_index(nupdates)
  integer*8 , intent(inout) :: later
  real*8    , intent(inout) :: determinant

  real*8 , dimension(lds, nupdates) :: Updates
  real*8 , dimension(dim, lds)      :: Inverse

  info = QMCKL_FAILURE

  ! Convert 'upds' and 's_inv' into the more easily readable Fortran
  ! matrices 'Updates' and 'Inverse'.
  call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)

  ! YET TO BE IMPLEMENTED

  ! Copy updated inverse back to s_inv
  call copy_back(Inverse, s_inv, lds, dim)

  info = QMCKL_SUCCESS

end function qmckl_slagel_splitting_doc_f

C interface to the pedagogical kernel (not directly exposed)

The following Fortran function qmckl_slagel_splitting_doc makes sure that the pedagogical kernel qmckl_slagel_splitting_doc_f, written in Fortran, can be called from C using the ISO_C_BINDING. The Fortran function qmckl_slagel_splitting_doc will be exposed in the header file 'qmckl.h' for C users and in the module file 'qmckl_f.F90' for Fortran users.

C headers (exposed in qmckl.h)

qmckl_exit_code qmckl_slagel_splitting (
      const uint64_t LDS,
      const uint64_t Dim,
      const uint64_t N_updates,
      const double* Updates,
      const uint64_t* Updates_index,
      const double breakdown,
      double* Slater_inv,
      double* later_updates,
      uint64_t* later_index,
      uint64_t* later,
      double* determinant );

C sources

qmckl_exit_code qmckl_slagel_splitting_hpc(
                  uint64_t LDS,
                  uint64_t Dim,
                  uint64_t N_updates,
                  const double* __restrict Updates,
                  const uint64_t* __restrict Updates_index,
                  const double breakdown,
                  double* __restrict Slater_inv,
                  double* __restrict later_updates,
                  uint64_t* __restrict later_index,
                  uint64_t* __restrict later,
                  double* __restrict determinant) {

  double __attribute__((aligned(8))) C[LDS];
  double __attribute__((aligned(8))) D[LDS];

  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.0f;
      IVDEP
      ALIGNED
      for (uint64_t j = 0; j < LDS; j++) {
        C[i] += Slater_inv[i * LDS + j] * Updates[l * LDS + j];
      }
    }

    // Denominator
    const int cui = Updates_index[l] - 1;
    double den = 1.0f + C[cui];
    if (fabs(den) < breakdown) {
      // U_l = U_l / 2: split the update in 2 equal halves and save the
      // second halve in later_updates
      IVDEP
      ALIGNED
      for (uint64_t i = 0; i < LDS; i++) {
        later_updates[*later * LDS + i] = Updates[l * LDS + i] * 0.5f;
        C[i] *= 0.5f;
      }
      later_index[*later] = Updates_index[l];
      (*later)++;

      den = 1.0f + C[cui];
    } // From here onwards we continue with applying the first halve of the
      // update to Slater_inv
    double iden = 1.0f / den;

    if (determinant)
      *determinant *= den;

    // D = v^T x S^{-1} : 1 x LDS
    IVDEP
    ALIGNED
    for (uint64_t j = 0; j < LDS; j++) {
      D[j] = Slater_inv[cui * LDS + j];
    }

    // S^{-1} = S^{-1} - C x D / den
    for (uint64_t i = 0; i < Dim; i++) {
      IVDEP
      ALIGNED
      for (uint64_t j = 0; j < LDS; j++) {
        const double update = C[i] * D[j] * iden;
        Slater_inv[i * LDS + j] -= update;
      }
    }
    l += 1;
  }

  return QMCKL_SUCCESS;
}

#+NAME:slagel_splitting_template_code

static inline qmckl_exit_code qmckl_slagel_splitting_{Dim}(
                                uint64_t N_updates,
                                const double* __restrict Updates,
                                const uint64_t* __restrict Updates_index,
                                const double breakdown,
                                double* __restrict Slater_inv,
                                double* __restrict later_updates,
                                uint64_t* __restrict later_index,
                                uint64_t* __restrict later,
                                double* __restrict determinant) {

  double __attribute__((aligned(8))) C[D{Dim}_P];
  double __attribute__((aligned(8))) D[D{Dim}_P];

  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.0f;
      IVDEP
      ALIGNED
      for (uint64_t j = 0; j < D{Dim}_P; j++) {
        C[i] += Slater_inv[i * D{Dim}_P + j] * Updates[l * D{Dim}_P + j];
      }
    }

    // Denominator
    const int cui = Updates_index[l] - 1;
    double den = 1.0f + C[cui];
    if (fabs(den) < breakdown) {
      // U_l = U_l / 2: split the update in 2 equal halves and save the
      // second halve in later_updates
      IVDEP
      ALIGNED
      for (uint64_t i = 0; i < D{Dim}_P; i++) {
        later_updates[*later * D{Dim}_P + i] = Updates[l * D{Dim}_P + i] * 0.5f;
        C[i] *= 0.5f;
      }
      later_index[*later] = Updates_index[l];
      (*later)++;

      den = 1.0f + C[cui];
    } // From here onwards we continue with applying the first halve of the
      // update to Slater_inv
    double iden = 1.0f / den;

    if (determinant)
      *determinant *= den;

    // D = v^T x S^{-1} : 1 x D{Dim}_P
    IVDEP
    ALIGNED
    for (uint64_t j = 0; j < D{Dim}_P; j++) {
      D[j] = Slater_inv[cui * D{Dim}_P + j];
    }

    // S^{-1} = S^{-1} - C x D / den
    for (uint64_t i = 0; i < {Dim}; i++) {
      IVDEP
      ALIGNED
      for (uint64_t j = 0; j < D{Dim}_P; j++) {
        const double update = C[i] * D[j] * iden;
        Slater_inv[i * D{Dim}_P + j] -= update;
      }
    }
    l += 1;
  }

  return QMCKL_SUCCESS;
}

#+NAME:slagel_splitting_kernel_generator

text="""
<<slagel_splitting_template_code>>
"""
result = []
for Dim in <<kernel_generator_range>>:
  Dim=str(Dim)
  result.append(text.replace("{Dim}",Dim) )

return '\n'.join(result)

#+NAME:slagel_splitting_switch-case_generator

text="""
case {Dim}:  
  return qmckl_slagel_splitting_{Dim}(
          N_updates,
          Updates,
          Updates_index,
          breakdown,
          Slater_inv,
          later_updates,
          later_index,
          later,
          determinant);
"""
result = []
for Dim in <<kernel_generator_range>>:
  Dim=str(Dim)
  result.append(text.replace("{Dim}",Dim) )

return '\n'.join(result)

<<slagel_splitting_kernel_generator()>>

qmckl_exit_code qmckl_slagel_splitting(
                      const uint64_t LDS,
                      const uint64_t Dim,
                      const uint64_t N_updates,
                      const double* Updates,
                      const uint64_t* Updates_index,
                      const double breakdown,
                      double* Slater_inv,
                      double* later_updates,
                      uint64_t* later_index,
                      uint64_t* later,
                      double* determinant) {

  #ifdef HAVE_HPC
  if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases
    switch (Dim) {
      <<slagel_splitting_switch-case_generator()>>
    }
  }
  else { // Updating smaller sub-matrix
    return qmckl_slagel_splitting_hpc(
                  LDS,
                  Dim,
                  N_updates,
                  Updates,
                  Updates_index,
                  breakdown,
                  Slater_inv,
                  later_updates,
                  later_index,
                  later,
                  determinant);
  }
  #else
  // return qmckl_slagel_splitting_doc(
  //               LDS,
  //               Dim,
  //               N_updates,
  //               Updates,
  //               Updates_index,
  //               breakdown,
  //               Slater_inv,
  //               later_updates,
  //               later_index,
  //               later,
  //               determinant);
  return qmckl_slagel_splitting_hpc(
                LDS,
                Dim,
                N_updates,
                Updates,
                Updates_index,
                breakdown,
                Slater_inv,
                later_updates,
                later_index,
                later,
                determinant);
  #endif

  return QMCKL_FAILURE;
}

Fortran interfaces (exposed in qmckl_f.F90)

Performance

This function cannot be used by itself and is used in Sherman-Morrison with update splitting and Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting. Please look at the performance reccomendations for those two kernels.

End of files

assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;

}