TREX/qmckl

mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-12-22 20:36:01 +01:00

François 4a9cecff61 Removed unnecessary __restrict keywords from public exposed functoin headers.

2023-01-26 11:50:58 +01:00

53 KiB

Raw Blame History

Sherman-Morrison-Woodbury

Headers
Naïve Sherman-Morrison
- qmckl_sherman_morrison
Woodbury 2x2
- qmckl_woodbury_2
Woodbury 3x3
- qmckl_woodbury_3
Sherman-Morrison with update splitting
- qmckl_sherman_morrison_splitting
Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting
- qmckl_sherman_morrison_smw32s
Helper Functions
- qmckl_slagel_splitting
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;

Naïve Sherman-Morrison

`qmckl_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.

qmckl_context	context	in	Global state
uint64_t	LDS	in	Leading dimension of Slater_inv
uint64_t	Dim	in	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[LDS*Dim]	inout	Array containing the inverse of a Slater-matrix
double*	determinant	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

C header

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

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

static inline qmckl_exit_code qmckl_sherman_morrison_{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_{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-python

#+NAME:naive-python-switch

<<naive-python()>>

qmckl_exit_code qmckl_sherman_morrison(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",
                      NULL);
}

if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases
switch (Dim) {
  <<naive-python-switch()>>
}
}
else { // When SIMD_LENGTH > 1, called with LDS == Dim AND Dim != (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH)
return qmckl_sherman_morrison_hpc(context,
                                  LDS,
                                  Dim,
	                                N_updates,
	                                Updates,
	                                Updates_index,
	                                breakdown,
	                                Slater_inv,
	                                determinant);

}

return QMCKL_SUCCESS;
}

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.

Woodbury 2x2

`qmckl_woodbury_2`

The Woodbury 2x2 kernel. It is used to apply two rank-1 updates at once. The formula used in this algorithm is called the Woodbury Matrix Identity \[ (S + U V)^{-1} = S^{-1} - C B^{-1} D \] where $S$ is the Slater-matrix $U$ and $V$ are the matrices containing the updates and the canonical basis matrix $S^{-1}$ is the inverse of the Slater-matrix $C:= S^{-1}U$, a Dim $\times 2$ matrix $B := 1 + VC$, the $2 \times 2$ matrix that is going to be inverted $D := VS^{-1}$, a $2 \times Dim$ matrix

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.

qmckl_context	context	in	Global state
uint64_t	LDS	in	Leading dimension of Slater_inv
uint64_t	Dim	in	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[LDS*Dim]	inout	Array containing the inverse of a Slater-matrix
double*	determinant	inout	Determinant of Slater-matrix

Requirements

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

C header

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

C source

qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
                            const uint64_t LDS,
                            const uint64_t Dim,
                            const double* Updates,
                            const uint64_t* Updates_index,
                            const double breakdown,
                            double* Slater_inv,
                            double* determinant) {
/*
C := S^{-1} * U,    dim x 2
B := 1 + V * C,     2 x 2
D := V * S^{-1},    2 x dim
*/

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

const uint64_t row1 = (Updates_index[0] - 1);
const uint64_t row2 = (Updates_index[1] - 1);

// Compute C = (S^T)^{-1}U : Dim x 2
double __attribute__((aligned(8))) C[2 * Dim];
for (uint64_t i = 0; i < Dim; i++) {
  C[i * 2] = 0;
  C[i * 2 + 1] = 0;
  IVDEP
  ALIGNED
  for (uint64_t k = 0; k < LDS; k++) {
      C[i * 2] += Slater_inv[i * LDS + k] * Updates[k];
      C[i * 2 + 1] += Slater_inv[i * LDS + k] * Updates[LDS + k];
  }
}

// Compute B = 1 + VC : 2 x 2
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;
}

// Update det(S) when passed
if (determinant)
  *determinant *= det;

// Compute B^{-1} with explicit formula for 2 x 2 inversion
double __attribute__((aligned(8))) 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;

// tmp = B^{-1}D : 2 x LDS
double __attribute__((aligned(8))) tmp[2 * LDS];
double* r1dim = &(Slater_inv[row1 * LDS]);
double* r2dim = &(Slater_inv[row2 * LDS]);
IVDEP
ALIGNED
for (uint64_t j = 0; j < LDS; j++) {
  tmp[j] = Binv[0] * r1dim[j] + Binv[1] * r2dim[j];
  tmp[LDS + j] = Binv[2] * r1dim[j] + Binv[3] * r2dim[j];
}

// Compute (S^T)^{-1} - C * tmp : Dim x LDS
for (uint64_t i = 0; i < Dim; i++) {
  IVDEP
  ALIGNED
  for (uint64_t j = 0; j < LDS; j++) {
      Slater_inv[i * LDS + j] -= C[i * 2] * tmp[j];
      Slater_inv[i * LDS + j] -= C[i * 2 + 1] * tmp[LDS + j];
  }
}

return QMCKL_SUCCESS;
}

Performance

This function is most efficient when used in cases where there are only 2 rank-1 updates and it is sure they will not result in a singular matrix.

Woodbury 3x3

`qmckl_woodbury_3`

The 3x3 version of the Woodbury 2x2 kernel. It is used to apply three rank-1 updates at once. The formula used in this kernel is the same as for Woodbury 2x2, except for the sizes of the following matrices:

$C:= S^{-1}U$, a Dim $\times 3$ matrix $B := 1 + VC$, the $3 \times 3$ matrix that is going to be inverted $D := VS^{-1}$, a $3 \times Dim$ matrix

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. #pragma ivdep #pragma vector aligned

qmckl_context	context	in	Global state
uint64_t	LDS	in	Leading dimension of Slater_inv
uint64_t	Dim	in	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[LDS*Dim]	inout	Array containing the inverse of a Slater-matrix
double*	determinant	inout	Determinant of Slater-matrix

Requirements

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

C header

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

C source

qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
                            const uint64_t LDS,
                            const uint64_t Dim,
                            const double* Updates,
                            const uint64_t* Updates_index,
                            const double breakdown,
                            double* Slater_inv,
                            double* determinant) {
/*
C := S^{-1} * U,    dim x 3
B := 1 + V * C,     3 x 3
D := V * S^{-1},    3 x dim
,*/

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

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^T)^{-1}U : Dim x 3
double __attribute__((aligned(8))) C[3 * Dim];
for (uint64_t i = 0; i < Dim; i++) {
  C[i * 3] = 0;
  C[i * 3 + 1] = 0;
  C[i * 3 + 2] = 0;
  IVDEP
  ALIGNED
  for (uint64_t k = 0; k < LDS; k++) {
      C[i * 3] += Slater_inv[i * LDS + k] * Updates[k];
      C[i * 3 + 1] += Slater_inv[i * LDS + k] * Updates[LDS + k];
      C[i * 3 + 2] += Slater_inv[i * LDS + k] * Updates[2 * LDS + k];
  }
}

// Compute B = 1 + VC : 3 x 3
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;
}

// Update det(Slater) if passed
if (determinant)
  *determinant *= det;

// Compute B^{-1} with explicit formula for 3 x 3 inversion
double __attribute__((aligned(8))) 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;

// tmp = B^{-1}D : 3 x LDS
double __attribute__((aligned(8))) tmp[3 * LDS];
double* r1dim = &(Slater_inv[row1 * LDS]);
double* r2dim = &(Slater_inv[row2 * LDS]);
double* r3dim = &(Slater_inv[row3 * LDS]);
  IVDEP
  ALIGNED
for (uint64_t j = 0; j < LDS; j++) {
  tmp[j] = Binv[0] * r1dim[j] + Binv[1] * r2dim[j] + Binv[2] * r3dim[j];
  tmp[LDS + j] =
      Binv[3] * r1dim[j] + Binv[4] * r2dim[j] + Binv[5] * r3dim[j];
  tmp[2 * LDS + j] =
      Binv[6] * r1dim[j] + Binv[7] * r2dim[j] + Binv[8] * r3dim[j];
}

// Compute (S^T)^{-1} - C * tmp : Dim x LDS
for (uint64_t i = 0; i < Dim; i++) {
  IVDEP
  ALIGNED
  for (uint64_t j = 0; j < LDS; j++) {
      Slater_inv[i * LDS + j] -= C[i * 3] * tmp[j];
      Slater_inv[i * LDS + j] -= C[i * 3 + 1] * tmp[LDS + j];
      Slater_inv[i * LDS + j] -= C[i * 3 + 2] * tmp[2 * LDS + j];
  }
}

return QMCKL_SUCCESS;
}

Performance…

This function is most efficient when used in cases where there are only 3 rank-1 updates and it is sure they will not result in a singular matrix.

Sherman-Morrison with update splitting

`qmckl_sherman_morrison_splitting`

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.

qmckl_context	context	in	Global state
uint64_t	LDS	in	Leading dimension of Slater_inv
uint64_t	Dim	in	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[LDS*Dim]	inout	Array containing the inverse of a Slater-matrix
double*	determinant	inout	Determinant of the Slater-matrix

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.

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

C header

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

C source

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

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

(void) qmckl_slagel_splitting(
  LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv,
  later_updates, later_index, &later, determinant);

if (later > 0) {
  (void) qmckl_sherman_morrison_splitting(
      context, LDS, Dim, later, later_updates, later_index, breakdown,
      Slater_inv, determinant);
}

return QMCKL_SUCCESS;
}

Performance…

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

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

`qmckl_sherman_morrison_smw32s`

The Woodbury 3x3 and 2x2 kernel with Sherman-Morrison and update splitting combines the low-level Woodbury 3x3 kernel, the Woobury 2x2 kernel and Sherman-Morrison with update splitting. It works the almost the same as Woodbury 3x3 with Sherman-Morrison and update splitting, except that when there is a remainder of two rank-1 updates, it is first tried with Woodbury 2x2 instead of sending them all to Sherman-Morrison with update splitting. For example, in the case of 5 updates the updates are applied in 1 block of 3 updates end 1 block of 2 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.

qmckl_context	context	in	Global state
uint64_t	LDS	in	Leading dimension of Slater_inv
uint64_t	Dim	in	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[LDS*Dim]	inout	Array containing the inverse of a Slater-matrix
double*	determinant	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

C header

qmckl_exit_code qmckl_sherman_morrison_smw32s(
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_smw32s(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_smw32s",
                      NULL);
}

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

// Special case for 4 rank-1 updates: 2+2
if (N_updates == 4) {
qmckl_exit_code rc =
    qmckl_woodbury_2(context, LDS, Dim, Updates, Updates_index,
                       breakdown, Slater_inv, determinant);
  if (rc != QMCKL_SUCCESS) { // Send the entire block to slagel_splitting
    uint64_t l = 0;
    rc = qmckl_slagel_splitting(LDS, Dim, 2, Updates, Updates_index,
                                breakdown, Slater_inv,
                                later_updates + (LDS * later),
                                later_index + later, &l, determinant);
    later += l;
  }
  rc = qmckl_woodbury_2(context, LDS, Dim, &Updates[2 * LDS],
                          &Updates_index[2], breakdown, Slater_inv,
                          determinant);
  if (rc != QMCKL_SUCCESS) { // Send the entire block to slagel_splitting
    uint64_t l = 0;
    rc = qmckl_slagel_splitting(
          LDS, Dim, 2, &Updates[2 * LDS], &Updates_index[2], breakdown,
          Slater_inv, later_updates + (LDS * later), later_index + later,
          &l, determinant);
    later += l;
  }
  if (later > 0) {
    rc = qmckl_sherman_morrison_splitting(
          context, LDS, Dim, later, later_updates, later_index, breakdown,
          Slater_inv, determinant);
  }
  return QMCKL_SUCCESS;
}

// And for the other cases != 4
// Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates
// with Woodbury 3x3 kernel
uint64_t n_of_3blocks = N_updates / 3;
uint64_t remainder = N_updates % 3;
uint64_t length_3block = 3 * LDS;

if (n_of_3blocks > 0) {
for (uint64_t i = 0; i < n_of_3blocks; i++) {
  const double* Updates_3block = &Updates[i * length_3block];
  const uint64_t* Updates_index_3block = &Updates_index[i * 3];
  qmckl_exit_code rc = qmckl_woodbury_3(
                        context, LDS, Dim, Updates_3block, Updates_index_3block,
                        breakdown, Slater_inv, determinant);
  if (rc != QMCKL_SUCCESS) { // Send the entire block to slagel_splitting
    uint64_t l = 0;
    rc = qmckl_slagel_splitting(
          LDS, Dim, 3, Updates_3block, Updates_index_3block,
          breakdown, Slater_inv, later_updates + (LDS * later),
          later_index + later, &l, determinant);
    later += l;
  }
}
}

// Apply last remaining block of 2 updates with Woodbury 2x2 kernel
if (remainder == 2) {
const double* Updates_2block = &Updates[n_of_3blocks * length_3block];
const uint64_t* Updates_index_2block = &Updates_index[3 * n_of_3blocks];
qmckl_exit_code rc = qmckl_woodbury_2(
                      context, LDS, Dim, Updates_2block, Updates_index_2block,
                      breakdown, Slater_inv, determinant);
if (rc != QMCKL_SUCCESS) { // Send the entire block to slagel_splitting
  uint64_t l = 0;
  rc = qmckl_slagel_splitting(
        LDS, Dim, 2, Updates_2block, Updates_index_2block, breakdown,
        Slater_inv, later_updates + (LDS * later), later_index + later,
        &l, determinant);
  later += l;
}
}

// Apply last remaining update with slagel_splitting
if (remainder == 1) {
const double* Updates_1block = &Updates[n_of_3blocks * length_3block];
const uint64_t* Updates_index_1block = &Updates_index[3 * n_of_3blocks];
uint64_t l = 0;
(void) qmckl_slagel_splitting(
                      LDS, Dim, 1, Updates_1block, Updates_index_1block, breakdown,
                      Slater_inv, later_updates + (LDS * later), later_index + later, &l,
                      determinant);
later += l;
}

if (later > 0) {
(void) qmckl_sherman_morrison_splitting(
                      context, LDS, Dim, later, later_updates, later_index, breakdown,
                      Slater_inv, determinant);
}

return QMCKL_SUCCESS;
}

Performance…

This kernel performs best for update cycles with 2 or more rank-1 updates and the fail-rate is low.

Helper Functions

Private helper-functions that are used by the Sherman-Morrison-Woodbury kernels. These functions can only be used internally by the kernels in this module.

`qmckl_slagel_splitting`

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.

uint64_t	LDS	in	Leading dimension of Slater_inv
uint64_t	Dim	in	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[LDS*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
double*	determinant	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

C header

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 source

qmckl_exit_code qmckl_slagel_splitting(uint64_t LDS,
                                   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,
                                   double* 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;
}

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;

}

53 KiB Raw Blame History Unescape Escape

Sherman-Morrison-Woodbury

Headers

Naïve Sherman-Morrison

qmckl_sherman_morrison

Requirements

C header

Performance

Woodbury 2x2

qmckl_woodbury_2

Requirements

C header

C source

Performance

Woodbury 3x3

qmckl_woodbury_3

Requirements

C header

C source

Performance…

Sherman-Morrison with update splitting

qmckl_sherman_morrison_splitting

Requirements

C header

C source

Performance…

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

qmckl_sherman_morrison_smw32s

Requirements

C header

C source

Performance…

Helper Functions

qmckl_slagel_splitting

Requirements

C header

C source

Performance

End of files

53 KiB

Raw Blame History

`qmckl_sherman_morrison`

`qmckl_woodbury_2`

`qmckl_woodbury_3`

`qmckl_sherman_morrison_splitting`

`qmckl_sherman_morrison_smw32s`

`qmckl_slagel_splitting`