qmckl/qmckl_sherman_morrison_woodbury.org at 0e4ea1f81b406b329279f4a7eedb8a585799dbeb

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;

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_slagel_splitting`

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.

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	breakdown	in	Break-down parameter on which to fail or not
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

qmckl_exit_code qmckl_slagel_splitting_c (
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 );

Source Fortran

Source C

#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,
                                          const double breakdown,
                                          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];
if (fabs(den) < breakdown) {

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

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	breakdown	in	Break-down parameter on which to fail or not
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,
const double breakdown,
double* Slater_inv );

Source Fortran

integer function qmckl_sherman_morrison_f(context, Dim, N_updates,    &
     Updates, Updates_index, breakdown, Slater_inv) 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(in) :: breakdown
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, breakdown, 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,
                            const double breakdown,
                            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];
if (fabs(den) < breakdown) {
  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	breakdown	in	Break-down parameter on which to fail or not
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,
const double breakdown,
double* Slater_inv );

Source Fortran

integer function qmckl_woodbury_2_f(context, Dim,    &
     Updates, Updates_index, breakdown, Slater_inv) 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(in) :: breakdown
real*8     , intent(inout)  :: Slater_inv(Dim*Dim)
!logical, external :: qmckl_woodbury_2_f
info = qmckl_woodbury_2(context, Dim, Updates, Updates_index, breakdown, 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,
                            const double breakdown,
                            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;
if (fabs(det) < breakdown) {
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	breakdown	in	Break-down parameter on which to fail or not
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,
const double breakdown,
double* Slater_inv );

Source Fortran

integer function qmckl_woodbury_3_f(context, Dim,    &
     Updates, Updates_index, breakdown, Slater_inv) 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(in) :: breakdown
real*8     , intent(inout)  :: Slater_inv(Dim*Dim)
!logical, external :: qmckl_woodbury_3_f
info = qmckl_woodbury_3(context, Dim, Updates, Updates_index, breakdown, 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,
                            const double breakdown,
                            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);
if (fabs(det) < breakdown) {
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…

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`

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	breakdown	in	Break-down parameter on which to fail or not
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_splitting_c (
const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv );

Source Fortran

integer function qmckl_sherman_morrison_splitting_f(context, Dim, N_updates,    &
     Updates, Updates_index, breakdown, Slater_inv) 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(in) :: breakdown
real*8     , intent(inout)  :: Slater_inv(Dim*Dim)
info = qmckl_sherman_morrison_splitting(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv)
end function qmckl_sherman_morrison_splitting_f

Source C

#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,
                            const double breakdown,
                            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

qmckl_context local_context;
local_context = qmckl_context_create();
qmckl_exit_code rc;

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

rc = qmckl_slagel_splitting_c(Dim, N_updates, Updates, Updates_index,
                        breakdown, Slater_inv, later_updates, later_index, &later);

if (later > 0) {
rc = qmckl_sherman_morrison_splitting_c(local_context, Dim, later, 
                        later_updates, later_index, breakdown, Slater_inv);
}

return QMCKL_SUCCESS;
}

Performance…

Woodbury 2x2 with Sherman-Morrison and 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_smw2s`

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	breakdown	in	Break-down parameter on which to fail or not
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_smw2s_c (
const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv );

Source Fortran

integer function qmckl_sherman_morrison_smw2s_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_smw2s (context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_smw2s_f

Source C

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

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

qmckl_context local_context;
local_context = qmckl_context_create();
qmckl_exit_code rc;

uint64_t n_of_2blocks = N_updates / 2;
uint64_t remainder = N_updates % 2;
uint64_t length_2block = 2 * Dim;

// Apply first 2*n_of_2blocks updates in n_of_2blocks blocks of 2 updates with
// Woodbury 2x2 kernel
double later_updates[Dim * N_updates];
uint64_t later_index[N_updates];
uint64_t later = 0;
if (n_of_2blocks > 0) {
for (uint64_t i = 0; i < n_of_2blocks; i++) {
  double *Updates_2block = &Updates[i * length_2block];
  uint64_t *Updates_index_2block = &Updates_index[i * 2];
  rc = qmckl_woodbury_2_c(local_context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv);
  if (rc != 0) { // Send the entire block to slagel_splitting
    uint64_t l = 0;
    rc = qmckl_slagel_splitting_c(Dim, 2, Updates_2block, Updates_index_2block,
            breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
    later = later + l;
  }
}
}

if (remainder == 1) { // Apply last remaining update with slagel_splitting
double *Updates_1block = &Updates[n_of_2blocks * length_2block];
uint64_t *Updates_index_1block = &Updates_index[2 * n_of_2blocks];
uint64_t l = 0;
rc = qmckl_slagel_splitting_c(Dim, 1, Updates_1block, Updates_index_1block,
        breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
later = later + l;
}

if (later > 0) {
rc = qmckl_sherman_morrison_splitting_c(local_context, Dim, later, later_updates, later_index, breakdown, Slater_inv);
}
return QMCKL_SUCCESS;
}

Performance…

Woodbury 3x3 with Sherman-Morrison and 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_smw3s`

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	breakdown	in	Break-down parameter on which to fail or not
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_smw3s_c (
const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv );

Source Fortran

integer function qmckl_sherman_morrison_smw3s_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_smw3s(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_smw3s_f

Source C

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

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

qmckl_context local_context;
local_context = qmckl_context_create();
qmckl_exit_code rc;

uint64_t n_of_3blocks = N_updates / 3;
uint64_t remainder = N_updates % 3;
uint64_t length_3block = 3 * Dim;

// Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with
// Woodbury 3x3 kernel
double later_updates[Dim * N_updates];
uint64_t later_index[N_updates];
uint64_t later = 0;
if (n_of_3blocks > 0) {
for (uint64_t i = 0; i < n_of_3blocks; i++) {
  double *Updates_3block = &Updates[i * length_3block];
  uint64_t *Updates_index_3block = &Updates_index[i * 3];
  rc = qmckl_woodbury_3_c(local_context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv);
  if (rc != 0) { // Send the entire block to slagel_splitting
    uint64_t l = 0;
    rc = qmckl_slagel_splitting_c(Dim, 3, Updates_3block, Updates_index_3block,
            breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
    later = later + l;
  }
}
}

if (remainder != 0) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
double *Updates_remainder_block = &Updates[n_of_3blocks * length_3block];
uint64_t *Updates_index_remainder_block = &Updates_index[3 * n_of_3blocks];
uint64_t l = 0;
rc = qmckl_slagel_splitting_c(Dim, remainder, Updates_remainder_block, Updates_index_remainder_block,
        breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
later = later + l;
}

if (later > 0) {
rc = qmckl_sherman_morrison_splitting_c(local_context, Dim, later, later_updates, later_index, breakdown, Slater_inv);
}
return QMCKL_SUCCESS;
}

Performance…

Woodbury 3x3 and 2x2 with Sherman-Morrison and 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_smw32s`

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	breakdown	in	Break-down parameter on which to fail or not
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_smw32s_c (
const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv );

Source Fortran

integer function qmckl_sherman_morrison_smw32s_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_smw32s(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_smw32s_f

Source C

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

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

qmckl_context local_context;
local_context = qmckl_context_create();
qmckl_exit_code rc;

uint64_t n_of_3blocks = N_updates / 3;
uint64_t remainder = N_updates % 3;
uint64_t length_3block = 3 * Dim;

// Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with
// Woodbury 3x3 kernel
double later_updates[Dim * N_updates];
uint64_t later_index[N_updates];
uint64_t later = 0;
if (n_of_3blocks > 0) {
for (uint64_t i = 0; i < n_of_3blocks; i++) {
  double *Updates_3block = &Updates[i * length_3block];
  uint64_t *Updates_index_3block = &Updates_index[i * 3];
  rc = qmckl_woodbury_3_c(local_context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv);
  if (rc != 0) { // Send the entire block to slagel_splitting
    uint64_t l = 0;
    rc = qmckl_slagel_splitting_c(Dim, 3, Updates_3block, Updates_index_3block,
            breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
    later = later + l;
  }
}
}

if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
double *Updates_2block = &Updates[n_of_3blocks * length_3block];
uint64_t *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
rc = qmckl_woodbury_2_c(local_context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv);
if (rc != 0) { // Send the entire block to slagel_splitting
  uint64_t l = 0;
  rc = qmckl_slagel_splitting_c(Dim, 2, Updates_2block, Updates_index_2block,
          breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
  later = later + l;
}
}
else if (remainder == 1) { // Apply last remaining update with slagel_splitting
double *Updates_1block = &Updates[n_of_3blocks * length_3block];
uint64_t *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
uint64_t l = 0;
rc = qmckl_slagel_splitting_c(Dim, 1, Updates_1block, Updates_index_1block,
        breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l);
later = later + l;
}

if (later > 0) {
rc = qmckl_sherman_morrison_splitting_c(local_context, Dim, later, later_updates, later_index, breakdown, Slater_inv);
}
return QMCKL_SUCCESS;
}

Performance…

End of files

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

62 KiB Raw Blame History

Sherman-Morrison-Woodbury

Headers

Helper Functions

qmckl_slagel_splitting

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…

Sherman-Morrison with update splitting

qmckl_sherman_morrison_splitting

Requirements

C header

Source Fortran

Source C

Performance…

Woodbury 2x2 with Sherman-Morrison and update splitting

qmckl_sherman_morrison_smw2s

Requirements

C header

Source Fortran

Source C

Performance…

Woodbury 3x3 with Sherman-Morrison and update splitting

qmckl_sherman_morrison_smw3s

Requirements

C header

Source Fortran

Source C

Performance…

Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting

qmckl_sherman_morrison_smw32s

Requirements

C header

Source Fortran

Source C

Performance…

End of files

62 KiB

Raw Blame History

`qmckl_slagel_splitting`

`qmckl_sherman_morrison`

`qmckl_woodbury_2`

`qmckl_woodbury_3`

`qmckl_sherman_morrison_splitting`

`qmckl_sherman_morrison_smw2s`

`qmckl_sherman_morrison_smw3s`

`qmckl_sherman_morrison_smw32s`