qmckl/org/qmckl_sherman_morrison_woodbury.org

Sherman-Morrison-Woodbury

• Headers
• Sherman-Morrison Helper Functions
  • qmckl_sherman_morrison_threshold
    • Requirements
    • C header
    • Source Fortran
    • Source C
    • Performance
• Naïve Sherman-Morrison
  • qmckl_sherman_morrison
    • Requirements
    • C header
    • Source Fortran
    • Source C
    • Performance
• Woodbury 2x2
  • qmckl_woodbury_2
    • Requirements
    • C header
    • Source Fortran
    • Source C
    • 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-singualr matrix

Headers

(org-babel-lob-ingest "../tools/lib.org")

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

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

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.

double

thresh

out

Threshold

Requirements

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

C header

// Sherman-Morrison-Woodbury break-down threshold
#ifndef THRESHOLD
#define THRESHOLD 1e-3
#endif

qmckl_exit_code qmckl_sherman_morrison_threshold_c (
double* const thresh );

Source Fortran

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

Source C

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

Performance

Naïve Sherman-Morrison

qmckl_sherman_morrison

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.

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

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

Source Fortran

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

Source C

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

Performance

Woodbury 2x2

[TODO: FMJC] Add main body intro.

qmckl_woodbury_2

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.

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

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

Source Fortran

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

Source C

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

Performance

End of files

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