TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-04 02:16:15 +02:00
Code Issues Projects Releases Wiki Activity
181f662c68
qmckl/org/qmckl_sherman_morrison_woodbury.org
Francois Coppens 181f662c68 Added macro HPC/DOC switch
2023-02-15 19:03:11 +01:00

1284 lines
44 KiB
Org Mode
Raw Blame History

#+TITLE: Sherman-Morrison-Woodbury
#+SETUPFILE: ../tools/theme.setup
#+INCLUDE: ../tools/lib.org
#+STARTUP: content
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
#+begin_src elisp :noexport :results none :exports none
(org-babel-lob-ingest "../tools/lib.org")
#+end_src
#+begin_src c :comments link :tangle (eval c_test) :noweb yes
#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;
#+end_src
This is the range that determines the how many high performance kernel instantces will be generated, using the C-function templates defined in the sections below. If the name of the C-function template is called ~qmckl_kernel_{Dim}~, then ~range(K, L+1)~ will results in kernel instances from ~qmckl_kernel_K~ to ~qmckl_kernel_L~.
#+NAME:kernel_generator_range
#+begin_src python :noweb yes :exports none
range(2, 22)
#+end_src
* Naïve Sherman-Morrison
** ~qmckl_sherman_morrison_naive~
:PROPERTIES:
:Name: qmckl_sherman_morrison_naive
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
*** Introduction
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.
#+TODO
Change the math notation so that the update vectors appear as row in the math
so that it is consistent with the representation in C (memory)
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.
*** API
#+NAME: qmckl_sherman_morrison_naive_args
| Variable | Type | In/Out | Description |
|-----------------+-------------------------+--------+------------------------------------------------------|
| ~context~ | ~qmckl_context~ | in | Global state |
| ~LDS~ | ~uint64_t~ | in | Leading dimension of Slater_inv |
| ~Dim~ | ~uint64_t~ | in | Dimension of Slater_inv |
| ~N_updates~ | ~uint64_t~ | in | Number of rank-1 updates to be applied to Slater_inv |
| ~Updates~ | ~double[N_updates*LDS]~ | in | Array containing the updates |
| ~Updates_index~ | ~uint64_t[N_updates]~ | in | Array containing the rank-1 updates |
| ~breakdown~ | ~double~ | in | Break-down parameter on which to fail or not |
| ~Slater_inv~ | ~double[Dim*LDS]~ | inout | Array containing the inverse of a Slater-matrix |
| ~determinant~ | ~double~ | inout | Determinant of the Slater-matrix |
*** Pedagogical kernel source (in Fortran)
The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is
able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore
not be used in real workloads.
#+begin_src f90 :tangle (eval f) :comment org :exports none
subroutine convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)
implicit none
integer*8 , intent(in) :: lds, dim, nupdates
real*8 , intent(in) :: upds(nupdates * lds)
real*8 , intent(in) :: s_inv(dim * lds)
real*8 , intent(out) , dimension(lds, nupdates) :: Updates
real*8 , intent(out) , dimension(dim, lds) :: Inverse
integer*8 :: i, j
! Construct Updates: lds x nupdates
do i = 1, nupdates
do j = 1, lds
Updates(j, i) = upds((i - 1) * lds + j)
end do
end do
! Construct Inverse: dim x lds
do i = 1, dim
do j = 1, lds
Inverse(i, j) = s_inv((i - 1) * lds + j)
end do
end do
end subroutine convert
#+end_src
#+begin_src f90 :tangle (eval f) :comment org :exports none
subroutine copy_back(Inverse, s_inv, lds, dim)
implicit none
integer*8 , intent(in) :: lds, dim
real*8 , intent(in) , dimension(dim, lds) :: Inverse
real*8 , intent(out) :: s_inv(dim * lds)
integer*8 :: i, j
! Copy updated inverse back to s_inv
do i = 1, dim
do j = 1, lds
s_inv((i - 1) * lds + j) = Inverse(i, j)
end do
end do
end subroutine copy_back
#+end_src
#+begin_src f90 :tangle (eval f)
integer function qmckl_sherman_morrison_naive_doc_f(context, &
lds, dim, &
nupdates, &
upds, &
updates_index, &
breakdown, &
s_inv, &
determinant) result(info)
use qmckl
implicit none
integer*8 , intent(in) :: context
integer*8 , intent(in) :: lds, dim
integer*8 , intent(in) :: nupdates
integer*8 , intent(in) :: updates_index(nupdates)
real*8 , intent(in) :: upds(nupdates * lds)
real*8 , intent(in) :: breakdown
real*8 , intent(inout) :: s_inv(dim * lds)
real*8 , intent(inout) :: determinant
real*8 , dimension(lds, nupdates) :: Updates
real*8 , dimension(dim, lds) :: Inverse
real*8 , dimension(dim) :: C
real*8 , dimension(lds) :: D
real*8 :: denominator, idenominator, update
integer*8 :: i, j, l, row
info = QMCKL_FAILURE
if (context == QMCKL_NULL_CONTEXT) then
info = QMCKL_INVALID_CONTEXT
return
endif
! Convert 'upds' and 's_inv' into the more easily readable Fortran
! matrices 'Updates' and 'Inverse'.
call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)
l = 1;
! For each update do...
do while (l < nupdates + 1)
! Compute C = S^{-1}U(l)
do i = 1, dim
C(i) = 0
do j = 1, dim
C(i) = C(i) + Inverse(i, j) * Updates(j, l)
end do
end do
! Compute denominator = 1 + V(l)^TC
row = updates_index(l)
denominator = 1 + C(row)
! Return early if denominator is too small
if (abs(denominator) < breakdown) return
idenominator = 1 / denominator
! Update det(S)
determinant = determinant * denominator
! selecting column: v_l^T * S_inv
D = Inverse(row, :)
! A^{-1} = A^{-1} - C x D / denominator
do i = 1, dim
do j = 1, dim
update = C(i) * D(j) * idenominator
Inverse(i, j) = Inverse(i, j) - update
end do
end do
l = l + 1
end do
! Copy updated inverse back to s_inv
call copy_back(Inverse, s_inv, lds, dim)
info = QMCKL_SUCCESS
end function qmckl_sherman_morrison_naive_doc_f
#+end_src
**** C interface to the pedagogical kernel (not directly exposed)
The following Fortran function ~qmckl_sherman_morrison_naive_doc~ makes sure
that the pedagogical kernel ~qmckl_sherman_morrison_naive_doc_f~, written in
Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function ~qmckl_sherman_morrison_naive_doc~ will be exposed in the header file 'qmckl.h'
for C users and in the module file 'qmckl_f.F90' for Fortran users.
#+CALL: generate_c_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_doc")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_sherman_morrison_naive_doc &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: determinant
integer(c_int32_t), external :: qmckl_sherman_morrison_naive_doc_f
info = qmckl_sherman_morrison_naive_doc_f &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant)
end function qmckl_sherman_morrison_naive_doc
#+end_src
*** Requirements
* ~context~ is not ~QMCKL_NULL_CONTEXT~
* ~LDS >= 2~
* ~Dim >= 2~
* ~N_updates >= 1~
* ~Updates~ is allocated with $N_updates \times Dim$ elements
* ~Updates_index~ is allocated with $N_updates$ elements
* ~breakdown~ is a small number such that $0 < breakdown << 1$
* ~Slater_inv~ is allocated with $Dim \times Dim$ elements
* ~determinant > 0~
*** C headers (exposed in qmckl.h)
#+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_sherman_morrison_naive (
const qmckl_context context,
const uint64_t LDS,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv,
double* determinant );
#+end_src
#+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_hpc")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_sherman_morrison_naive_hpc (
const qmckl_context context,
const uint64_t LDS,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv,
double* determinant );
#+end_src
#+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_doc")
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_sherman_morrison_naive_doc (
const qmckl_context context,
const uint64_t LDS,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv,
double* determinant );
#+end_src
*** C sources
Common includes and macros used by all the Sherman-Morrison-Woodbury kernels.
#+begin_src c :tangle (eval c) :comments org
#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
#+end_src
~qmckl_sherman_morrison_naive_hpc~ is a high performance variation of
~qmckl_sherman_morrison_naive~ written in C. It is used in cases when ~Dim~ is
smaller than the leading dimension ~LDS~, irrespective of whetether ~LDS~
includes zero padding to benefit from SIMD instructions or not. Cases like this
include situations where one wants to apply updates to a square submatrix of the
full matrix.
It takes advantage of memory aligned data and assumes no data dependencies
inside the loops. The loops are fully vectorised whenever ~Dim~ is an integer
multiple of ~SIMD_LEGTH~.
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code qmckl_sherman_morrison_naive_hpc(
const qmckl_context context,
const uint64_t LDS,
const uint64_t Dim,
const uint64_t N_updates,
const double* __restrict Updates,
const uint64_t* __restrict Updates_index,
const double breakdown,
double* __restrict Slater_inv,
double* __restrict determinant) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( context,
QMCKL_NULL_CONTEXT,
"qmckl_sherman_morrison_naive_hpc",
NULL);
}
double __attribute__((aligned(8))) C[Dim];
double __attribute__((aligned(8))) D[LDS];
uint64_t l = 0;
// For each update
while (l < N_updates) {
// C = S^{-1} x u_l
for (uint64_t i = 0; i < Dim; i++) {
C[i] = 0.0f;
IVDEP
ALIGNED
for (uint64_t j = 0; j < Dim; 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 < Dim; 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 < Dim; j++) {
const double update = C[i] * D[j] * iden;
Slater_inv[i * LDS + j] -= update;
}
}
l += 1;
}
return QMCKL_SUCCESS;
}
#+end_src
~qmckl_exit_code qmckl_sherman_morrison_naive_{Dim}~ is a C function-template that is used to genereate instances of C fucntions based on the range given above. The advantage of this method is that for each of these instances all the dimensions and loop-bounds are known at compile time, allowing the compiler to optimize more aggressively.
#+NAME:naive_template_code
#+begin_src c
static inline qmckl_exit_code qmckl_sherman_morrison_naive_{Dim}(
const qmckl_context context,
const uint64_t N_updates,
const double* __restrict Updates,
const uint64_t* __restrict Updates_index,
const double breakdown,
double* __restrict Slater_inv,
double* __restrict determinant) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith(context,
QMCKL_NULL_CONTEXT,
"qmckl_sherman_morrison_naive_{Dim}",
NULL);
}
#define D{Dim}_P ((1+({Dim}-1)/SIMD_LENGTH)*SIMD_LENGTH)
double __attribute__((aligned(8))) C[{Dim}];
double __attribute__((aligned(8))) D[D{Dim}_P];
uint64_t l = 0;
// For each update
while (l < N_updates) {
// C = A^{-1} x U_l
for (uint64_t i = 0; i < {Dim}; i++) {
C[i] = 0;
IVDEP
ALIGNED
for (uint64_t j = 0; j < D{Dim}_P; j++) {
C[i] += Slater_inv[i * D{Dim}_P + j] * Updates[l * D{Dim}_P + j];
}
}
// Denominator
const int cui = Updates_index[l] - 1;
double den = 1.0f + C[cui];
if (fabs(den) < breakdown) {
return QMCKL_FAILURE;
}
double iden = 1.0f / den;
// Update det(A)
if (determinant)
*determinant *= den;
// selecting column: D = v_l^T * S_inv
IVDEP
ALIGNED
for (uint64_t j = 0; j < D{Dim}_P; j++) {
D[j] = Slater_inv[cui * D{Dim}_P + j];
}
// A^{-1} = A^{-1} - C x D / den
for (uint64_t i = 0; i < {Dim}; i++) {
IVDEP
ALIGNED
for (uint64_t j = 0; j < D{Dim}_P; j++) {
double update = C[i] * D[j] * iden;
Slater_inv[i * D{Dim}_P + j] -= update;
}
}
l += 1;
}
return QMCKL_SUCCESS;
}
#+end_src
This is the kernel generator written in Python. It uses the kernel generator range and templates defined above to generate the C kernel instances.
#+NAME:naive_kernel_generator
#+begin_src python :noweb yes
text="""
<<naive_template_code>>
"""
result = []
for Dim in <<kernel_generator_range>>:
Dim=str(Dim)
result.append(text.replace("{Dim}",Dim) )
return '\n'.join(result)
#+end_src
Python script that generated C switch cases that call individual kernel instances.
#+NAME:naive_switch-case_generator
#+begin_src python :noweb yes
text="""
case {Dim}:
return qmckl_sherman_morrison_naive_{Dim}(context,
N_updates,
Updates,
Updates_index,
breakdown,
Slater_inv,
determinant);
"""
result = []
for Dim in <<kernel_generator_range>>:
Dim=str(Dim)
result.append(text.replace("{Dim}",Dim) )
return '\n'.join(result)
#+end_src
#+begin_src c :tangle (eval c) :comments org :noweb yes
<<naive_kernel_generator()>>
#+end_src
~qmckl_sherman_morrison_naive~ is a generic function that contains decision making logic that calls the proper kernel based on the used library configuration (~--enable-doc~ and ~--enable-hpc~) and the passed array dimensions ~LDS~ and ~Dim~.
#+begin_src c :tangle (eval c) :comments org :noweb yes
qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context,
const uint64_t LDS,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv,
double* determinant) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith(context,
QMCKL_NULL_CONTEXT,
"qmckl_sherman_morrison_naive",
NULL);
}
#ifdef HAVE_HPC
if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases
switch (Dim) {
<<naive_switch-case_generator()>>
}
}
else { // Updating smaller sub-matrix
return qmckl_sherman_morrison_naive_hpc(context,
LDS,
Dim,
N_updates,
Updates,
Updates_index,
breakdown,
Slater_inv,
determinant);
}
#else
return qmckl_sherman_morrison_naive_doc(context,
LDS,
Dim,
N_updates,
Updates,
Updates_index,
breakdown,
Slater_inv,
determinant);
#endif
return QMCKL_FAILURE;
}
#+end_src
*** Fortran interfaces (exposed in qmckl_f.F90)
:PROPERTIES:
:Name: qmckl_sherman_morrison_naive
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive_hpc")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison_naive_hpc &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: determinant
end function qmckl_sherman_morrison_naive_hpc
end interface
#+end_src
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive_doc")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison_naive_doc &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: determinant
end function qmckl_sherman_morrison_naive_doc
end interface
#+end_src
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison_naive &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: determinant
end function qmckl_sherman_morrison_naive
end interface
#+end_src
*** 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.
*** Tests
The tests for the kernels are executed on datasets that are extracted from a run of
QMC=Chem on Benzene (21 spin-up/21 spin down electrons) using 329 unique alpha determinants.
The tests are run such that the kernels reject the computed inverse whenever the computed
intermediate determinants or denominators are smaller than 1e-3. This is the default value in
QMC=Chem. The tests will return QMCKL_SUCCESS whenever all the elements of the final matrix
$R=S.S^-1 - 1$ are smaller than the given tolerance value of 1e-3, and will return
QMCKL_FAILURE if the values are larger than this tolerance value.
#+begin_src c :tangle (eval c_test)
const uint64_t Dim = 21;
const uint64_t LDS = (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH;
const double breakdown = 1e-3;
const double tolerance = 1e-3;
double res[441];
#include "sm_test.h"
assert(Updates1 != NULL);
assert(Updates_index1 != NULL);
assert(Slater_inv1 != NULL);
// original determinant of Slater1 (before applying updates)
double det = 3.407025646103221e-10;
rc = qmckl_sherman_morrison_naive(context,
LDS,
Dim,
N_updates1,
Updates1,
Updates_index1,
breakdown,
Slater_inv1,
&det);
// Check that the determinant is updated properly
assert(fabs(det + 4.120398385068217e-10) < 1e-15);
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
res[i * Dim + j] = 0;
for (unsigned int k = 0; k < Dim; k++) {
res[i * Dim + j] += Slater1[i * Dim + k] * Slater_inv1[k * LDS + j];
}
}
}
rc = QMCKL_SUCCESS;
for (unsigned int i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) {
if (i == j && fabs(res[i * Dim + j] - 1) > tolerance) {
rc = QMCKL_FAILURE;
}
if (i != j && fabs(res[i * Dim + j]) > tolerance) {
rc = QMCKL_FAILURE;
}
}
}
assert(rc == QMCKL_SUCCESS);
#+end_src
* 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~
:PROPERTIES:
:Name: qmckl_slagel_splitting
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
*** Introduction ~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.
*** API
#+NAME: qmckl_slagel_splitting_args
| Variable | Type | In/Out | Description |
|-----------------+-------------------------+--------+---------------------------------------------------------------|
| ~LDS~ | ~uint64_t~ | in | Leading dimension of Slater_inv |
| ~Dim~ | ~uint64_t~ | in | Dimension of Slater_inv |
| ~N_updates~ | ~uint64_t~ | in | Number of rank-1 updates to be applied to Slater_inv |
| ~Updates~ | ~double[N_updates*LDS]~ | in | Array containing the rank-1 updates |
| ~Updates_index~ | ~uint64_t[N_updates]~ | in | Array containing positions of the rank-1 updates |
| ~breakdown~ | ~double~ | in | Break-down parameter on which to fail or not |
| ~Slater_inv~ | ~double[Dim*LDS]~ | inout | Array containing the inverse Slater-matrix |
| ~later_updates~ | ~double[N_updates*LDS]~ | inout | Array containing the split updates for later |
| ~later_index~ | ~uint64_t[N_updates]~ | inout | Array containing the positions of the split updates for later |
| ~later~ | ~uint64_t~ | inout | Number of split updates for later |
| ~determinant~ | ~double~ | inout | Determinant of the Slater-matrix |
*** Pedagogical kernel source (in Fortran)
The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is
able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore
not be used in real workloads.
#+begin_src f90 :tangle (eval f)
integer function qmckl_slagel_splitting_doc_f( &
lds, dim, &
nupdates, &
upds, &
updates_index, &
breakdown, &
s_inv, &
later_updates, &
later_index, &
later, &
determinant) result(info)
use qmckl
implicit none
integer*8 , intent(in) :: lds, dim
integer*8 , intent(in) :: nupdates
integer*8 , intent(in) :: updates_index(nupdates)
real*8 , intent(in) :: upds(nupdates * lds)
real*8 , intent(in) :: breakdown
real*8 , intent(inout) :: s_inv(dim * lds)
real*8 , intent(inout) :: later_updates(nupdates * lds)
integer*8 , intent(inout) :: later_index(nupdates)
integer*8 , intent(inout) :: later
real*8 , intent(inout) :: determinant
real*8 , dimension(lds, nupdates) :: Updates
real*8 , dimension(dim, lds) :: Inverse
info = QMCKL_FAILURE
! Convert 'upds' and 's_inv' into the more easily readable Fortran
! matrices 'Updates' and 'Inverse'.
call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)
! YET TO BE IMPLEMENTED
! Copy updated inverse back to s_inv
call copy_back(Inverse, s_inv, lds, dim)
info = QMCKL_SUCCESS
end function qmckl_slagel_splitting_doc_f
#+end_src
**** C interface to the pedagogical kernel (not directly exposed)
The following Fortran function ~qmckl_slagel_splitting_doc~ makes sure
that the pedagogical kernel ~qmckl_slagel_splitting_doc_f~, written in
Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function
~qmckl_slagel_splitting_doc~ will be exposed in the header file 'qmckl.h'
for C users and in the module file 'qmckl_f.F90' for Fortran users.
#+CALL: generate_c_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_slagel_splitting_doc")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_slagel_splitting_doc &
(LDS, &
Dim, &
N_updates, &
Updates, &
Updates_index, &
breakdown, &
Slater_inv, &
later_updates, &
later_index, &
later, &
determinant) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: later_updates(N_updates*LDS)
integer (c_int64_t) , intent(inout) :: later_index(N_updates)
integer (c_int64_t) , intent(inout) :: later
real (c_double ) , intent(inout) :: determinant
integer(c_int32_t), external :: qmckl_slagel_splitting_doc_f
info = qmckl_slagel_splitting_doc_f &
(LDS, &
Dim, &
N_updates, &
Updates, &
Updates_index, &
breakdown, &
Slater_inv, &
later_updates, &
later_index, &
later, &
determinant)
end function qmckl_slagel_splitting_doc
#+end_src
*** 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 headers (exposed in qmckl.h)
#+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
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 );
#+end_src
*** C sources
#+begin_src c :tangle (eval c) :comments org
qmckl_exit_code qmckl_slagel_splitting_hpc(
uint64_t LDS,
uint64_t Dim,
uint64_t N_updates,
const double* __restrict Updates,
const uint64_t* __restrict Updates_index,
const double breakdown,
double* __restrict Slater_inv,
double* __restrict later_updates,
uint64_t* __restrict later_index,
uint64_t* __restrict later,
double* __restrict 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;
}
#+end_src
#+NAME:slagel_splitting_template_code
#+begin_src c
static inline qmckl_exit_code qmckl_slagel_splitting_{Dim}(
uint64_t N_updates,
const double* __restrict Updates,
const uint64_t* __restrict Updates_index,
const double breakdown,
double* __restrict Slater_inv,
double* __restrict later_updates,
uint64_t* __restrict later_index,
uint64_t* __restrict later,
double* __restrict determinant) {
double __attribute__((aligned(8))) C[D{Dim}_P];
double __attribute__((aligned(8))) D[D{Dim}_P];
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 < 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) {
// 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 < D{Dim}_P; i++) {
later_updates[*later * D{Dim}_P + i] = Updates[l * D{Dim}_P + 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 D{Dim}_P
IVDEP
ALIGNED
for (uint64_t j = 0; j < D{Dim}_P; j++) {
D[j] = Slater_inv[cui * D{Dim}_P + 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 < D{Dim}_P; j++) {
const double update = C[i] * D[j] * iden;
Slater_inv[i * D{Dim}_P + j] -= update;
}
}
l += 1;
}
return QMCKL_SUCCESS;
}
#+end_src
#+NAME:slagel_splitting_kernel_generator
#+begin_src python :noweb yes
text="""
<<slagel_splitting_template_code>>
"""
result = []
for Dim in <<kernel_generator_range>>:
Dim=str(Dim)
result.append(text.replace("{Dim}",Dim) )
return '\n'.join(result)
#+end_src
#+NAME:slagel_splitting_switch-case_generator
#+begin_src python :noweb yes
text="""
case {Dim}:
return qmckl_slagel_splitting_{Dim}(
N_updates,
Updates,
Updates_index,
breakdown,
Slater_inv,
later_updates,
later_index,
later,
determinant);
"""
result = []
for Dim in <<kernel_generator_range>>:
Dim=str(Dim)
result.append(text.replace("{Dim}",Dim) )
return '\n'.join(result)
#+end_src
#+begin_src c :tangle (eval c) :comments org :noweb yes
<<slagel_splitting_kernel_generator()>>
#+end_src
#+begin_src c :tangle (eval c) :comments org :noweb yes
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) {
#ifdef HAVE_HPC
if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases
switch (Dim) {
<<slagel_splitting_switch-case_generator()>>
}
}
else { // Updating smaller sub-matrix
return qmckl_slagel_splitting_hpc(
LDS,
Dim,
N_updates,
Updates,
Updates_index,
breakdown,
Slater_inv,
later_updates,
later_index,
later,
determinant);
}
#else
// return qmckl_slagel_splitting_doc(
// LDS,
// Dim,
// N_updates,
// Updates,
// Updates_index,
// breakdown,
// Slater_inv,
// later_updates,
// later_index,
// later,
// determinant);
return qmckl_slagel_splitting_hpc(
LDS,
Dim,
N_updates,
Updates,
Updates_index,
breakdown,
Slater_inv,
later_updates,
later_index,
later,
determinant);
#endif
return QMCKL_FAILURE;
}
#+end_src
*** Fortran interfaces (exposed in qmckl_f.F90)
:PROPERTIES:
:Name: qmckl_slagel_splitting
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_hpc")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_slagel_splitting_hpc &
(LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: later_updates(N_updates*LDS)
integer (c_int64_t) , intent(inout) :: later_index(N_updates)
integer (c_int64_t) , intent(inout) :: later
real (c_double ) , intent(inout) :: determinant
end function qmckl_slagel_splitting_hpc
end interface
#+end_src
#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_doc")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_slagel_splitting_doc &
(LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: later_updates(N_updates*LDS)
integer (c_int64_t) , intent(inout) :: later_index(N_updates)
integer (c_int64_t) , intent(inout) :: later
real (c_double ) , intent(inout) :: determinant
end function qmckl_slagel_splitting_doc
end interface
#+end_src
#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_slagel_splitting &
(LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: later_updates(N_updates*LDS)
integer (c_int64_t) , intent(inout) :: later_index(N_updates)
integer (c_int64_t) , intent(inout) :: later
real (c_double ) , intent(inout) :: determinant
end function qmckl_slagel_splitting
end interface
#+end_src
*** 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
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;
}
#+end_src
# -*- mode: org -*-
# vim: syntax=c
Reference in New Issue View Git Blame Copy Permalink