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
• Woodbury 3x3
  • qmckl_woodbury_3
    • 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

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

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

Performance

Woodbury 2x2

This is the Woodbury 3x3 kernel.

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 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_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 (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;
}

Performance

Woodbury 3x3

This is the Woodbury 3x3 kernel.

qmckl_woodbury_3

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[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

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

Source Fortran

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

Source C

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

Performance…

End of files

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