qmckl/org/qmckl_sherman_morrison_woodbury.org

Sherman-Morrison-Woodbury

• Headers
• Naïve Sherman-Morrison
  • Pedagogical kernel source (in Fortran)
    • C interface to the pedagogical kernel
  • Requirements
  • C headers (exposed in qmckl.h)
  • C sources
  • Fortran interfaces (exposed in qmckl_f.F90)
  • Tests
• 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;

#+NAME:kernel_generator_range

Naïve Sherman-Morrison

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.

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.

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

Pedagogical kernel source (in Fortran)

integer function qmckl_sherman_morrison_naive_doc_f(context, &
    LDS, Dim, &
    N_updates, &
    Updates, &
    Updates_index, &
    breakdown, &
    Slater_inv, &
    determinant) result(info)

  use qmckl
  implicit none
  integer*8 , intent(in)    :: context
  integer*8 , intent(in)    :: LDS, Dim
  integer*8 , intent(in)    :: N_updates
  integer*8 , intent(in)    :: Updates_index(N_updates)
  real*8    , intent(in)    :: Updates(N_updates*LDS)
  real*8    , intent(in)    :: breakdown
  real*8    , intent(inout) :: Slater_inv(Dim*LDS)
  real*8    , intent(inout) :: determinant

  info = 0

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

  write(*,*) "Function 'qmckl_sherman_morrison_naive_doc_f' does nothing for now..."

  info = QMCKL_SUCCESS

end function qmckl_sherman_morrison_naive_doc_f

C interface to the pedagogical kernel

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

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

#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_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 < LDS; 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 < LDS; 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 < LDS; j++) {
        const double update = C[i] * D[j] * iden;
        Slater_inv[i * LDS + j] -= update;
      }
    }
    l += 1;
  }
  return QMCKL_SUCCESS;
}

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

#+NAME:naive_kernel_generator

#+NAME:naive_switch-case_generator

<<naive_kernel_generator()>>

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 { // When SIMD_LENGTH > 1, called with LDS == Dim AND Dim != (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH)
    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)

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

End of files

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

}