2021-07-19 12:01:07 +02:00
#+TITLE : Sherman-Morrison-Woodbury
#+SETUPFILE : ../tools/theme.setup
#+INCLUDE : ../tools/lib.org
2021-09-14 09:55:55 +02:00
#+STARTUP : content
2021-07-19 12:01:07 +02:00
2021-09-07 12:22:39 +02:00
Low- and high-level functions that use the Sherman-Morrison and
2021-07-21 17:56:04 +02:00
Woodbury matrix inversion formulas to update the inverse of a
2021-09-07 09:27:22 +02:00
non-singular matrix
2021-09-07 11:22:54 +02:00
2021-07-19 12:01:07 +02:00
* Headers
2021-07-21 17:56:04 +02:00
#+begin_src elisp :noexport :results none :exports none
2021-07-19 12:01:07 +02:00
(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
2021-07-20 11:15:33 +02:00
#include <math.h >
2021-07-19 17:38:53 +02:00
2021-07-19 12:01:07 +02:00
int main() {
qmckl_context context;
context = qmckl_context_create();
2021-07-20 19:34:51 +02:00
qmckl_exit_code rc;
2021-07-19 12:01:07 +02:00
#+end_src
2021-10-06 23:44:06 +02:00
2021-07-21 17:30:12 +02:00
* Naïve Sherman-Morrison
2021-10-13 10:45:54 +02:00
2021-07-19 12:01:07 +02:00
** ~qmckl_sherman_morrison~
:PROPERTIES:
:Name: qmckl_sherman_morrison
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-09-07 11:22:54 +02:00
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
2021-09-01 16:06:51 +02:00
\[
2021-09-01 15:29:14 +02:00
(S + uv^T)^{-1} = S^{-1} - \frac{S^{-1} uv^T S^{-1}}{1 + v^T S^{-1} u}
2021-09-01 16:06:51 +02:00
\]
2021-09-01 15:29:14 +02:00
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.
2021-07-19 12:01:07 +02:00
2021-09-07 11:22:54 +02:00
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$,
2021-09-07 12:22:39 +02:00
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}.
2021-09-07 11:22:54 +02:00
2021-10-13 10:31:24 +02:00
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.
2021-07-19 12:01:07 +02:00
#+NAME: qmckl_sherman_morrison_args
2021-07-19 18:25:10 +02:00
| qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv |
2021-10-13 13:19:35 +02:00
| uint64_t | N_updates | in | Number of rank-1 updates to be applied to Slater_inv |
2021-07-19 18:25:10 +02:00
| double | Updates[N_updates*Dim] | in | Array containing the updates |
| uint64_t | Updates_index[N_updates] | in | Array containing the rank-1 updates |
2021-07-26 17:41:21 +02:00
| double | breakdown | in | Break-down parameter on which to fail or not |
2021-07-19 18:25:10 +02:00
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
2021-10-13 12:03:18 +02:00
| double* | determinant | inout | Determinant of the Slater-matrix |
2021-07-19 12:01:07 +02:00
*** Requirements
2021-10-06 23:44:06 +02:00
2021-09-07 11:22:54 +02:00
* ~context~ is not ~QMCKL_NULL_CONTEXT~
* ~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
2021-10-06 23:44:06 +02:00
2021-07-19 12:01:07 +02:00
*** C header
2021-07-19 14:00:10 +02:00
2021-07-19 14:04:47 +02:00
#+CALL: generate_c_header(table=qmckl_sherman_morrison_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
2021-07-19 14:00:10 +02:00
2021-07-19 17:38:53 +02:00
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
2021-09-16 13:44:12 +02:00
qmckl_exit_code qmckl_sherman_morrison(
2021-07-19 18:25:10 +02:00
const qmckl_context context,
2021-10-08 18:03:11 +02:00
const uint64_t Dim,
const uint64_t N_updates,
2021-07-19 18:25:10 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-08 18:03:11 +02:00
const double breakdown,
2021-10-13 10:31:24 +02:00
double* Slater_inv,
2021-10-13 12:03:18 +02:00
double* determinant);
2021-07-19 18:25:10 +02:00
#+end_src
2021-10-09 22:23:12 +02:00
*** C source
2021-07-19 18:25:10 +02:00
#+begin_src c :tangle (eval c) :comments org
#include <stdbool.h >
2021-09-10 17:25:56 +02:00
#include <math.h >
2021-07-19 18:25:10 +02:00
#include "qmckl.h"
2021-10-08 18:03:11 +02:00
qmckl_exit_code qmckl_sherman_morrison(const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
2021-07-19 18:25:10 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-08 18:03:11 +02:00
const double breakdown,
2021-10-13 10:31:24 +02:00
double* Slater_inv,
double* determinant) {
2021-10-06 23:44:06 +02:00
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
2021-07-19 17:38:53 +02:00
double C[Dim];
double D[Dim];
2021-07-22 09:59:02 +02:00
uint64_t l = 0;
2021-07-19 17:38:53 +02:00
// For each update
while (l < N_updates) {
// C = A^{-1} x U_l
2021-07-22 09:59:02 +02:00
for (uint64_t i = 0; i < Dim; i++) {
2021-07-19 17:38:53 +02:00
C[i] = 0;
2021-07-22 09:59:02 +02:00
for (uint64_t j = 0; j < Dim; j++) {
2021-07-19 17:38:53 +02:00
C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
}
}
// Denominator
double den = 1 + C[Updates_index[l] - 1];
2021-10-13 10:31:24 +02:00
2021-07-26 17:41:21 +02:00
if (fabs(den) < breakdown) {
2021-07-20 12:09:43 +02:00
return QMCKL_FAILURE;
2021-07-19 17:38:53 +02:00
}
double iden = 1 / den;
2021-10-13 10:31:24 +02:00
// Update det(A)
if (determinant != NULL)
*determinant *= den;
2021-07-19 17:38:53 +02:00
// D = v^T x A^{-1}
2021-07-22 09:59:02 +02:00
for (uint64_t j = 0; j < Dim; j++) {
2021-07-19 17:38:53 +02:00
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
}
// A^{-1} = A^{-1} - C x D / den
2021-07-22 09:59:02 +02:00
for (uint64_t i = 0; i < Dim; i++) {
for (uint64_t j = 0; j < Dim; j++) {
2021-07-19 17:38:53 +02:00
double update = C[i] * D[j] * iden;
Slater_inv[i * Dim + j] -= update;
}
}
l += 1;
}
2021-07-29 11:48:38 +02:00
2021-07-19 18:25:10 +02:00
return QMCKL_SUCCESS;
2021-07-19 17:38:53 +02:00
}
#+end_src
2021-07-19 12:01:07 +02:00
*** Performance
2021-09-07 12:22:39 +02:00
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.
2021-09-01 15:29:14 +02:00
2021-10-09 22:23:12 +02:00
** Fortran interface :noexport:
2021-09-14 12:13:17 +02:00
:PROPERTIES:
:Name: qmckl_sherman_morrison
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-10-06 23:44:06 +02:00
2021-07-19 14:04:47 +02:00
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
2021-09-15 10:52:00 +02:00
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison &
2021-10-13 10:31:24 +02:00
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
2021-09-15 10:52:00 +02:00
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
2021-09-18 18:07:05 +02:00
integer (c_int64_t) , intent(in) , value :: context
2021-10-08 18:03:11 +02:00
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
2021-09-18 18:07:05 +02:00
real (c_double ) , intent(in) :: Updates(N_updates*Dim)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
2021-10-08 18:03:11 +02:00
real (c_double ) , intent(in) , value :: breakdown
2021-09-18 18:07:05 +02:00
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
2021-10-13 10:31:24 +02:00
real (c_double ) , intent(inout) :: determinant
2021-09-15 10:52:00 +02:00
end function qmckl_sherman_morrison
end interface
#+end_src
2021-07-19 14:04:47 +02:00
*** Test :noexport:
2021-09-05 11:22:41 +02:00
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.
2021-07-20 19:34:51 +02:00
#+begin_src c :tangle (eval c_test)
2021-09-05 11:22:41 +02:00
const uint64_t Dim = 21;
2021-07-26 17:41:21 +02:00
const double breakdown = 1e-3;
2021-09-05 11:22:41 +02:00
const double tolerance = 1e-3;
double res[441];
const uint64_t N_updates1 = 1;
const uint64_t Updates_index1[1] = {21};
const double Updates1[21] = {-0.0068946573883295389, -0.024496730417013196, -0.2336846236139537, -0.011720726499334019, 0.094674163497984229, 0.049728685989975902, 0.025834936182945929, 0.0586009658873081, -0.06664357893168929, -0.060625394340604558, 0.1100576994940636, -0.077925791963934898, -0.07255193404853344, 0.04257450252771329, -0.079073280096053591, 0.053145386278629192, 0.014488439075648779, -0.070511359721422195, -0.13172587938606731, -0.01733743725344538, 0.26186770945787408};
const double Slater1[441] = {-2.8945870399475102, 3.5455725193023699, 2.0470812320709202, -3.5464441776275599, -2.0474903583526598, -2.89596366882324, -0.61329728364944502, 0.70991641283035301, 0.45664468407630898, 0.59523195028305098, 0.26079276204109197, -0.027727147564291999, -0.35093706846237199, -0.095610238611698206, -0.130077064037323, 0.10946778208017301, 0.021800471469759899, 0.048480678349733401, -0.092234551906585693, -0.0160505045205355, 0.0065372241660952603, -0.78021097183227495, -0.95558542013168302, -0.55174458026885997, -0.95579683780670199, -0.55179446935653698, 0.78044348955154397, -0.115299895405769, -0.11754634231329, -0.0318448171019554, -0.00082233443390577999, 0.064001239836216001, 0.12740932404995001, -0.049412809312343597, -0.074464745819568606, 0.19851659238338501, -0.088878795504569993, 0.135610401630402, -0.13736967742443101, -0.101879440248013, 0.102390937507153, 0.061368178576231003, -0.00527230883017182, 0.00029381641070358499, 0.0072919740341603799, -0.00026352208806201799, 0.0068838247098028703, 0.0047294595278799499, 0.11728889495134399, -0.062190085649490398, -0.15969239175319699, -0.10645916312933, 0.112943567335606, -0.0015487050404772199, -0.045737843960523598, 0.13196058571338701, 0.038660705089569099, -0.050266433507204097, -0.13745591044426, -0.037740666419267703, 0.11781705915927899, 0.19417031109332999, -0.0186148826032877, -0.30874741077423101, -0.37815779447555498, 0.218254208564758, 0.378139227628708, -0.21839827299118, -0.30879178643226601, 0.025539221242070202, 0.112043909728527, 0.037676706910133403, 0.025347150862216901, -0.19991625845432301, 0.14163509011268599, 0.13326919078826899, 0.213842943310738, 0.131471157073975, 0.14626120030879999, -0.0118067460134625, 0.093547157943248693, 0.22439678013324699, 0.0082425400614738499, -0.00068585784174501896, -1.52881526947021, -1.9163775505148799e-05, -2.1623263359069802, -1.5540548702119799e-05, 2.16284847259521, -1.52953028678894, -0.26402890682220498, 0.024944067001342801, -0.37140902876853898, -0.044071510434150703, -0.29843890666961698, -0.043834518641233403, -0.177285701036453, -0.058585006743669503, -0.0183692276477814, -0.026075478643178902, -0.15623773634433699, -0.011319605633616401, 0.054057534784078598, -0.23151709139347099, 0.10071966797113401, -0.067333526909351293, 0.00043256155913695698, 0.095035225152969402, -0.00043632232700474598, 0.094624765217304202, 0.066768139600753798, 0.14448715746402699, -0.12090456485748299, -0.203432962298393, -0.21472349762916601, 0.154288679361343, 0.017385326325893399, -0.071518480777740506, 0.27918133139610302, 0.080454014241695404, -0.0361245274543762, 0.0566458702087402, -0.030157707631588, 0.25150132179260298, -0.082445412874221802, -0.017144737765192999, -2.2779068946838401, -2.3393469746224599e-05, -3.2218937873840301, -1.7939200915861899e-05, 3.22278833389282, -2.2791447639465301, -0.47945570945739702, 0.045197278261184699, -0.64478874206543002, -0.080234304070472703, -0.47389149665832497, 0.024908870458602898, -0.23455648124218001, -0.10695217549800901, -0.033054649829864502, -0.195749551057816, -0.040076982229948002, 0.094596333801746396, 0.098401188850402804, -0.059563819319009802, 0.028169330209493599, -0.0025739683769643298, -0.0020174346864223502, -0.0024467820767313199, -0.00204527890309691, 0.00013988610589876801, 0.00073755031917244196, 0.17339251935482, 0.099962860345840496, 0.14622049033641801, -0.088966980576515198, 0.040044382214546197, -0.14437201619148299, -1.5416861060657501e-05, -0.081159926950931494, -0.11909016221761699, 0.170753449201584, 0.031124599277973199, -0.052109345793724102, 0.121445283293724, 0.0289684906601906, 0.0189813487231731, -0.0043673445470631097, -0.00014539182302542001, 0.0061219600029289696, 0.000189467871678062, 0.0060523007996380303, 0.00425094133242965, 0.0905451029539108, 0.058540407568216303, -0.15191346406936601, 0.11511342227459, 0.15016922354698201, 0.10024280846118901, -0.11066282540559801, -0.15089659392833699, -0.021995550021529201, 0.10021472722291901, 0.019932290539145501, 0.070281185209751101, -0.120184391736984, -0.030
double Slater_inv1[441] = {-0.0525082963130096, -192.960477122746, -167.343357271998, 1.69063377058866, -1378.4882402189, 239.360652372684, 874.40192183061697, -513.31846478968203, 807.35483045933495, -19.026422090297299, -17.824591748935202, 165.6645926386, -332.909289060628, -18.6107453916377, -11.6439930601979, 317.18559450822198, 214.12532640356201, -430.80943299473, -344.65005510575799, 740.81014093255203, -84.923325284754696, 0.064335245989650805, 234.71397330210399, 205.1186636226, -3.4153770027768502, 1680.13890084581, -293.40065007215702, -1065.8419078679401, 624.86259814866003, -989.76729197131897, 23.298249286417001, 21.573213707404399, -202.41174445438401, 405.86421806904502, 22.922945938042901, 14.236362125535999, -386.65825218682801, -261.06337565261498, 528.20337075570103, 421.94063881187702, -902.75507799084596, 105.03985525163, 0.044200668381308601, -138.85152101714399, -121.27002019292, 2.0061662585667599, -992.69022926063894, 173.24384564131401, 629.63518796071799, -369.68832373377103, 582.99742710801104, -13.474468953452901, -12.4187314531643, 119.002229516372, -239.76386110607999, -13.4968690647037, -8.1499663302329708, 228.499936475189, 154.09732894897101, -311.21651724011002, -248.26815844269399, 533.37480106787098, -61.991136485706697, -0.0765497914729725, 240.772681546397, 209.115172008988, -2.12401058528051, 1722.15420882834, -299.07379390801702, -1092.4670851015201, 641.09661184520598, -1007.44097957056, 23.333413806919001, 21.9897514440317, -206.53257016044699, 415.81861108697302, 23.3743606050752, 14.107639254523001, -396.52477956346598, -267.77854087146301, 537.68662188870303, 429.83533323670702, -925.65731644586003, 106.070284414301, -0.037109975520249097, -136.738457148886, -118.917063549582, 1.19704609156923, -977.22998951029899, 169.63000196416201, 620.01253859081601, -364.02043006288699, 570.85982682561803, -13.451266756380299, -12.146454471873501, 117.435790064483, -235.687342903548, -13.183975346097, -8.2288378293500593, 224.66974233421601, 151.73703725446899, -304.613662218336, -243.40661913323399, 525.55077671846004, -60.142862384234697, -0.062476778936209301, 195.74640208908301, 170.54051550313201, -2.8328765312753998, 1399.12576024585, -244.28628529648901, -887.63077081451104, 520.884480704167, -823.92607500410702, 19.0479002693959, 18.195743596123101, -168.030291787452, 337.586201042798, 18.991146825575299, 11.5224539224244, -321.78784999564101, -217.23371709456501, 439.82292557372398, 351.28300986264401, -752.091241883468, 87.457548717202599, 4.3572379528639997e-05, -23.790404094575401, -19.504220922406098, 0.26595531593048999, -161.34910319146201, 30.163156787264001, 102.337776116854, -58.672051676794297, 103.227708063905, -0.22802366652163999, -2.2306129106953398, 21.778747237132102, -40.783330246607299, -1.38732280869919, -0.75357165711863805, 39.269880423319599, 26.312297182321299, -55.061514171221297, -43.569942636118199, 85.973367077774896, -10.016061028885201, 4.0427159001010801e-05, 28.392374743152502, 26.373254860492299, -0.36980926031100297, 209.22398819265001, -37.9473984838014, -132.76472733653301, 77.720615205630693, -127.924723556082, 0.29745864731730998, 2.7793514223195599, -29.019189334306802, 53.826499495808498, 4.2582570650914597, 0.51526747227423098, -50.914050124737102, -31.724328055085401, 68.0683646567051, 55.266627392249198, -110.63482872154501, 14.1152956761065, 0.00055720168222872903, -2.6345228739497299, -3.4087974363069198, 0.048641237778327802, -18.3623844122964, 3.8195276251277801, 11.6411523542838, -5.8812824661877396, 15.735088038012799, -0.058858662594294599, -0.23845780430312999, 6.5124119812799099, -6.1494099868393697, 0.340564236694477, -0.880332042883611, 5.7252222161672197, 2.8545094725865199, -8.3447036743893506, -8.0497901808946395, 9.3503553819376606, -1.76541796469797, 0.00098753020434672305, -4.6360828892267198, -4.7638931470942598, 0.052453544062221497, -34.4916331053766, 7.7511206271403301, 21.997597399118899, -10.6934789694381, 32.345600850876203, -0.10706305963427901, -0.54284354835464799, 12.037079682564, -12.6961140374279, -1.7388
const uint64_t Updates_index2[2] = {20, 21};
const double Updates2[42] = {-0.012056605890393198, 0.072118259966373, -0.14986032247543299, -0.023892195895314251, 0.2306191368843431, 0.070976656861603302, 0.059228850266663378, -0.0081594046205282003, 0.03775668889284136, -0.32248186320066397, -0.30175055470317586, -0.048578753136098399, -0.050512738234829151, 0.074964269995690003, 0.022375471889972701, 0.20063611492514641, -0.013817949220538101, 0.034675425849854905, 0.14269321411848099, 0.16303040878847203, 0.20520395040512121, 0.027209745720028898, -0.019752152264117806, 0.24686626344919221, 0.021544795483350802, 0.0053255971870385145, -0.074287162162363515, 0.0020896954229101539, -0.0015901625156403004, 0.045653678011149175, -0.095019596163183437, 0.015724316006526361, 0.057632524985820083, 0.0043175870669074383, -0.075611688196659296, -0.2787729352712634, 0.071667610667645931, -0.0016040895134210986, -0.018613442778587338, 0.16365851461887351, 0.018654582090675831, -0.43227782845497098};
const double Slater2[441] = {-2.8945870399475102, 3.5455725193023699, 2.0470812320709202, -3.5464441776275599, -2.0474903583526598, -2.89596366882324, -0.61329728364944502, 0.70991641283035301, 0.45664468407630898, 0.59523195028305098, 0.26079276204109197, -0.027727147564291999, -0.35093706846237199, -0.095610238611698206, -0.130077064037323, 0.10946778208017301, 0.021800471469759899, 0.048480678349733401, -0.092234551906585693, -0.028107110410928698, 0.013431881554424799, -0.78021097183227495, -0.95558542013168302, -0.55174458026885997, -0.95579683780670199, -0.55179446935653698, 0.78044348955154397, -0.115299895405769, -0.11754634231329, -0.0318448171019554, -0.00082233443390577999, 0.064001239836216001, 0.12740932404995001, -0.049412809312343597, -0.074464745819568606, 0.19851659238338501, -0.088878795504569993, 0.135610401630402, -0.13736967742443101, -0.101879440248013, 0.174509197473526, 0.085864908993244199, -0.00527230883017182, 0.00029381641070358499, 0.0072919740341603799, -0.00026352208806201799, 0.0068838247098028703, 0.0047294595278799499, 0.11728889495134399, -0.062190085649490398, -0.15969239175319699, -0.10645916312933, 0.112943567335606, -0.0015487050404772199, -0.045737843960523598, 0.13196058571338701, 0.038660705089569099, -0.050266433507204097, -0.13745591044426, -0.037740666419267703, 0.11781705915927899, 0.044309988617896999, 0.215069741010666, -0.30874741077423101, -0.37815779447555498, 0.218254208564758, 0.378139227628708, -0.21839827299118, -0.30879178643226601, 0.025539221242070202, 0.112043909728527, 0.037676706910133403, 0.025347150862216901, -0.19991625845432301, 0.14163509011268599, 0.13326919078826899, 0.213842943310738, 0.131471157073975, 0.14626120030879999, -0.0118067460134625, 0.093547157943248693, 0.22439678013324699, -0.015649655833840401, 0.011034868657589, -1.52881526947021, -1.9163775505148799e-05, -2.1623263359069802, -1.5540548702119799e-05, 2.16284847259521, -1.52953028678894, -0.26402890682220498, 0.024944067001342801, -0.37140902876853898, -0.044071510434150703, -0.29843890666961698, -0.043834518641233403, -0.177285701036453, -0.058585006743669503, -0.0183692276477814, -0.026075478643178902, -0.15623773634433699, -0.011319605633616401, 0.054057534784078598, -0.00089795450912788499, 0.0060455044731497799, -0.067333526909351293, 0.00043256155913695698, 0.095035225152969402, -0.00043632232700474598, 0.094624765217304202, 0.066768139600753798, 0.14448715746402699, -0.12090456485748299, -0.203432962298393, -0.21472349762916601, 0.154288679361343, 0.017385326325893399, -0.071518480777740506, 0.27918133139610302, 0.080454014241695404, -0.0361245274543762, 0.0566458702087402, -0.030157707631588, 0.25150132179260298, -0.0114687560126185, -0.066873423755168901, -2.2779068946838401, -2.3393469746224599e-05, -3.2218937873840301, -1.7939200915861899e-05, 3.22278833389282, -2.2791447639465301, -0.47945570945739702, 0.045197278261184699, -0.64478874206543002, -0.080234304070472703, -0.47389149665832497, 0.024908870458602898, -0.23455648124218001, -0.10695217549800901, -0.033054649829864502, -0.195749551057816, -0.040076982229948002, 0.094596333801746396, 0.098401188850402804, -0.000334969052346423, 0.0023343940265476699, -0.0025739683769643298, -0.0020174346864223502, -0.0024467820767313199, -0.00204527890309691, 0.00013988610589876801, 0.00073755031917244196, 0.17339251935482, 0.099962860345840496, 0.14622049033641801, -0.088966980576515198, 0.040044382214546197, -0.14437201619148299, -1.5416861060657501e-05, -0.081159926950931494, -0.11909016221761699, 0.170753449201584, 0.031124599277973199, -0.052109345793724102, 0.121445283293724, 0.020809086039662399, -0.039619617164135, -0.0043673445470631097, -0.00014539182302542001, 0.0061219600029289696, 0.000189467871678062, 0.0060523007996380303, 0.00425094133242965, 0.0905451029539108, 0.058540407568216303, -0.15191346406936601, 0.11511342227459, 0.15016922354698201, 0.10024280846118901, -0.11066282540559801, -0.15089659392833699, -0.021995550021529201, 0.10021472722291901, 0.019932290539145501, 0.070281185209751101, -0.120184391736984, 0.0072
double Slater_inv2[441] = {-0.064091478343647895, 713.06177859755496, -52.018444001735702, 0.22294142549212401, 6739.1399551272798, 409.47332011602299, -4258.7085329747897, 2990.6198857571399, 2842.82174111627, -16.4362076124787, -19.6479024764561, -450.95698589104001, 1499.21617070053, 129.354086480939, -32.5623649949015, -1533.16729952379, -801.75588044433505, -1302.1490458327801, -1589.22747455529, -3895.5139631694801, 69.692448899038894, 0.078717515330272506, -890.24945003874404, 61.925371835460197, -1.59301536151858, -8399.1221177311709, -504.62121099099301, 5307.6897279898503, -3725.8058949477299, -3517.1067799943698, 20.0821069369808, 23.8371292623433, 563.21704392130198, -1868.99615781019, -160.79773282728499, 40.209679175758602, 1910.8342353503599, 1000.30652221449, 1610.10305651065, 1967.2714135880001, 4853.9412977750599, -86.938972564264105, 0.035772543187836597, 520.38606836760005, -37.3576026553948, 0.93824771966655895, 4913.8375731575898, 297.02078853525001, -3105.3054726097698, 2179.8383098375198, 2064.03863882426, -11.5897837400758, -13.7454041271756, -329.662381563542, 1093.32256576446, 94.164922311323394, -23.370538453071202, -1117.8490893585799, -585.07554862150903, -945.21841098621996, -1153.8443875258399, -2840.0954494666398, 50.509996472385303, -0.062171023850398299, -883.91684171449401, 65.956744158981707, -0.30209264289181997, -8354.6527651622891, -510.24292797392201, 5279.5127573481996, -3708.5126037463401, -3534.1651244869099, 20.118054506627601, 24.253115792914802, 558.90980701110504, -1858.4878939499599, -160.30158683160701, 40.074632460370502, 1900.40832678587, 993.28424539653497, 1619.3228924922901, 1974.7898591375399, 4829.6374496633298, -85.861801422967005, -0.045296585245108302, 503.608089249118, -37.409231418514402, 0.15972971393766999, 4760.0409737735199, 289.86000580815897, -3007.9000554557201, 2112.4472365954498, 2009.46047743586, -11.620588884143499, -13.435110153350299, -318.37193513541001, 1059.1982996055399, 91.392671066014998, -23.013250839168599, -1083.0984474979, -566.25418012298701, -920.44766441625598, -1123.03275109381, -2751.24966532286, 49.134449946833101, -0.050574070850706303, -735.268788752521, 52.034324888924203, -1.3246973887217, -6942.4299970826296, -419.091569409214, 4387.0782989889804, -3079.7112127626101, -2915.54203952635, 16.3862338355292, 20.069350976326799, 465.60100478088498, -1545.0790688177201, -133.05534575663401, 33.017866331913801, 1579.60766417078, 826.67092177544703, 1335.1987409164999, 1630.19252491066, 4012.1274550981102, -71.423354704319806, -0.00110547248442243, 66.086472294581398, -8.0640571635796796, 0.120361096092389, 643.91485283159295, 47.038233657028698, -406.863770451839, 288.91657765213398, 305.14483180816802, 0.028924090527909001, -2.4114842655012301, -39.3897517070086, 140.96244082771301, 13.290701732299899, -2.82866188450281, -144.284041440354, -74.462522726873601, -141.49789272768501, -167.03129270601499, -373.94726486427999, 5.3217330886365399, 0.00169058806401361, -100.681182204012, 9.9438606966790992, -0.16071917417936801, -947.22789644116904, -62.1819544991477, 598.507402482923, -421.456732120239, -417.90099350750501, -0.071547956877298299, 3.03910354158063, 58.825827638175099, -207.18138234350499, -16.8210784946065, 3.4953362885853601, 212.690544720802, 112.99993227712, 192.20097004017899, 232.57136070342199, 549.86422579734096, -7.9115450643297303, 0.00031167490112610501, 16.5702800434129, -0.96427444933400797, 0.017530807484877599, 153.70562899881901, 7.4253775224466096, -97.164535381519997, 68.391115404813903, 58.880541272935403, -0.00395431753418613, -0.27710621824447701, -6.5580129619132999, 32.685848366303901, 3.4769501442949502, -1.3237352902339601, -33.496400026660503, -18.678955342199799, -26.814344384641199, -34.430890540207301, -88.9250347843996, 1.51194685980572, 0.00045214452900354401, 37.241126324245002, 0.56653388165982099, -0.0153845901652373, 340.71285233348101, 15.613889705904301, -215.25964333041301, 151.26187968817001, 126.42681477968, 0.012659113756634199, -0.62711870198964903, -16.4637553538752, 71.9864648940
const uint64_t N_updates3 = 3;
const uint64_t Updates_index3[3] = {12, 13, 21};
const double Updates3[63] = {-0.32320992089807998, -0.1768221333622936, -0.044189138920046375, -0.0083658993244170032, -0.13345118239521958, -0.088903807103633908, -0.25946535170078289, 0.14435659933042233, -0.21090563386678701, 0.019084170460701003, 0.073579406365752206, -0.013528384268283893, 0.031103432178496981, -0.25866857916116798, 0.022464506328106093, 0.032419085502625011, -0.083496315404772786, -0.392626002430916, -0.057606873102486092, 0.1326765250414613, 0.034384071826935009, 0.25532682985067379, -0.025051936507225009, 0.17769842967391061, 0.080573752522469011, 0.11870069429278349, 0.35069981217384349, 0.127604305744171, -0.081144510089870836, -0.04023376852273898, 0.0045312587171792984, 0.017237693071365301, 0.0386847481131554, -0.21935135498642899, 0.040930040180683996, -0.057743255048990395, -0.289656192064286, 0.26965582184493592, 0.257760990411043, -0.1886596158146856, -0.036238497123122194, -0.1386065091937782, 0.0015524076297879202, -0.1212832294404507, 0.028310360386967659, -0.005844349274411801, 0.0046449899673459971, 0.044258039444684996, 0.0025386139750480999, 0.0354412645101548, 0.064816890284419101, 0.1017611473798752, -0.25162110477685901, -0.0285851191729307, -0.024413086473941803, -0.22420015186071351, 0.059393912553786996, 0.04153373837471, 0.01215143408626318, 0.203658737242222, -0.0020150262862444011, -0.037645858246833121, 0.0714646056294439};
const double Slater3[441] = {-2.8945870399475102, 3.5455725193023699, 2.0470812320709202, -3.5464441776275599, -2.0474903583526598, -2.89596366882324, -0.61329728364944502, 0.70991641283035301, 0.45664468407630898, 0.59523195028305098, 0.26079276204109197, -0.35093706846237199, -0.095610238611698206, -0.130077064037323, 0.10946778208017301, 0.021800471469759899, 0.048480678349733401, -0.092234551906585693, -0.0160505045205355, -0.028107110410928698, 0.0080896317958831804, -0.78021097183227495, -0.95558542013168302, -0.55174458026885997, -0.95579683780670199, -0.55179446935653698, 0.78044348955154397, -0.115299895405769, -0.11754634231329, -0.0318448171019554, -0.00082233443390577999, 0.064001239836216001, -0.049412809312343597, -0.074464745819568606, 0.19851659238338501, -0.088878795504569993, 0.135610401630402, -0.13736967742443101, -0.101879440248013, 0.102390937507153, 0.174509197473526, -0.0599150508642197, -0.00527230883017182, 0.00029381641070358499, 0.0072919740341603799, -0.00026352208806201799, 0.0068838247098028703, 0.0047294595278799499, 0.11728889495134399, -0.062190085649490398, -0.15969239175319699, -0.10645916312933, 0.112943567335606, -0.045737843960523598, 0.13196058571338701, 0.038660705089569099, -0.050266433507204097, -0.13745591044426, -0.037740666419267703, 0.11781705915927899, 0.19417031109332999, 0.044309988617896999, 0.0096954777836799604, -0.30874741077423101, -0.37815779447555498, 0.218254208564758, 0.378139227628708, -0.21839827299118, -0.30879178643226601, 0.025539221242070202, 0.112043909728527, 0.037676706910133403, 0.025347150862216901, -0.19991625845432301, 0.13326919078826899, 0.213842943310738, 0.131471157073975, 0.14626120030879999, -0.0118067460134625, 0.093547157943248693, 0.22439678013324699, 0.0082425400614738499, -0.015649655833840401, -0.00653020711615682, -1.52881526947021, -1.9163775505148799e-05, -2.1623263359069802, -1.5540548702119799e-05, 2.16284847259521, -1.52953028678894, -0.26402890682220498, 0.024944067001342801, -0.37140902876853898, -0.044071510434150703, -0.29843890666961698, -0.177285701036453, -0.058585006743669503, -0.0183692276477814, -0.026075478643178902, -0.15623773634433699, -0.011319605633616401, 0.054057534784078598, -0.23151709139347099, -0.00089795450912788499, 0.10536465793848, -0.067333526909351293, 0.00043256155913695698, 0.095035225152969402, -0.00043632232700474598, 0.094624765217304202, 0.066768139600753798, 0.14448715746402699, -0.12090456485748299, -0.203432962298393, -0.21472349762916601, 0.154288679361343, -0.071518480777740506, 0.27918133139610302, 0.080454014241695404, -0.0361245274543762, 0.0566458702087402, -0.030157707631588, 0.25150132179260298, -0.082445412874221802, -0.0114687560126185, 0.027113301679492, -2.2779068946838401, -2.3393469746224599e-05, -3.2218937873840301, -1.7939200915861899e-05, 3.22278833389282, -2.2791447639465301, -0.47945570945739702, 0.045197278261184699, -0.64478874206543002, -0.080234304070472703, -0.47389149665832497, -0.23455648124218001, -0.10695217549800901, -0.033054649829864502, -0.195749551057816, -0.040076982229948002, 0.094596333801746396, 0.098401188850402804, -0.059563819319009802, -0.000334969052346423, 0.030707944184541699, -0.0025739683769643298, -0.0020174346864223502, -0.0024467820767313199, -0.00204527890309691, 0.00013988610589876801, 0.00073755031917244196, 0.17339251935482, 0.099962860345840496, 0.14622049033641801, -0.088966980576515198, 0.040044382214546197, -1.5416861060657501e-05, -0.081159926950931494, -0.11909016221761699, 0.170753449201584, 0.031124599277973199, -0.052109345793724102, 0.121445283293724, 0.0289684906601906, 0.020809086039662399, 0.0544226132333279, -0.0043673445470631097, -0.00014539182302542001, 0.0061219600029289696, 0.000189467871678062, 0.0060523007996380303, 0.00425094133242965, 0.0905451029539108, 0.058540407568216303, -0.15191346406936601, 0.11511342227459, 0.15016922354698201, -0.11066282540559801, -0.15089659392833699, -0.021995550021529201, 0.10021472722291901, 0.019932290539145501, 0.070281185209751101, -0.120184391736984, -0.030551897361874601, 0.00720479153
// WB3
double Slater_inv3_1[441] = {-0.056658742514349103, -9.5006825962797503, -18.307026635935799, -0.36744251865365501, 31.610841135834299, 39.084671241308698, -19.703711132730799, 18.666497646833299, 121.964198363159, -17.598444551917201, -2.9520245875548698, 24.935160344122199, 13.992695020707, 4.6612706754314903, -3.6234738869455199, -28.604054847833801, 9.5739197944338503, -62.128443913882499, -67.232459151277695, -11.803301377693099, -3.8508731471217099, 0.069411805979426894, 10.8007994955779, 22.458394727678201, -0.89282708659453403, -39.476164388726097, -47.5420201141931, 24.536662733275399, -23.351558477080602, -146.60044409871099, 21.549928185153298, 3.32891560872004, -30.5549979706051, -17.384760133954, -5.4401084934617003, 4.3730845094487396, 35.2040348414587, -11.414817838980399, 74.926509290741393, 80.145978783950099, 14.6842509146363, 5.7893824160342398, 0.041199702875412497, -6.4779318938824302, -13.2989719358296, 0.51508213587065199, 23.9434304832055, 27.923633944040301, -14.994539048994801, 13.545190256632999, 84.658750799602501, -12.440996315879, -1.6347815496642299, 17.4053340515119, 10.4563880132908, 3.2716511123291698, -2.32038524814726, -20.901080430435499, 6.5092932119827198, -43.310144838525801, -46.266317618221201, -9.0209731319320596, -3.3217327746699499, -0.071364036152653701, 11.529909516191401, 22.9175425489932, 0.44722818276140403, -39.897585842875301, -48.876557451246697, 24.801523764393401, -23.690956189203401, -151.28132599138999, 21.549304148995901, 3.4094704881876998, -30.687597458642099, -17.659848105023201, -5.7069215817131997, 4.0886157921745996, 35.564950307785701, -12.1803912964186, 77.135406641174299, 83.318689013925706, 14.819424258558801, 4.7782939163062803, -0.040051226770134098, -6.7134150132482402, -13.312257348225099, -0.26126985985289403, 22.204441530888801, 27.729627986407898, -13.7012421876807, 13.0503428483308, 85.297842345116607, -12.439369240025201, -1.60848061970372, 17.698309764178902, 10.1798985121263, 3.3110252774145601, -2.5466287558690199, -20.410023365852201, 6.7633439222022202, -43.4149505982517, -46.886928763683102, -7.88950712988734, -2.6927901019362199, -0.058248999176268398, 9.0508747193022892, 18.587941141225699, -0.73447131422595202, -35.308936282117699, -39.940782969048001, 21.9119312895712, -20.076960036338601, -123.935102576568, 17.592669318251101, 3.0259816420954402, -24.782453059479, -15.379073099754899, -4.67601503080067, 3.33242427220405, 30.0344780145867, -9.0772191462285008, 63.391820359557002, 67.760837774760006, 13.3701295015636, 4.8413052345891403, -0.00048534604109837199, -1.8238859168577699, 0.56645976558050504, -0.011652978445607001, 3.3511865791586, 2.0262681174419299, -2.14462816857628, 1.8811114733258301, 1.72516289469276, -0.0410444803602956, -0.17992011162600399, 4.6474809764852196, 0.32638302042170703, 1.28068103014527, 0.42381953652975901, -1.62158383882094, 1.84081037263479, -1.2636267813034701, -0.97726719185308597, -0.84319863660797401, 0.57245121780013697, 0.00072965273109526496, -0.74868100456805098, 0.61375144696985096, -0.013670796650539699, -10.880212250529199, -2.2389072194253701, 6.8446360871752603, -3.8338819099152301, -0.83708262317409499, 0.055638693619984997, 0.16395892268443801, -6.4020540512225601, -0.87619791103918299, 0.67273193621704896, -0.96217527481654197, 3.5322304037025498, 0.74780382069687101, 0.457651060635482, 2.4127372990709501, 5.8207812126355103, 0.41153697505715597, 0.00047300381075295099, -0.24273841428794099, 0.62931133864105804, -0.0075433586410305498, -3.8733549821913198, -2.7042122832831499, 2.4066006453807001, -1.90501007795962, -11.4961539239149, -0.025181201641384401, 0.20951155742128399, 4.41315756942062, -2.02827911285326, 0.53245648205106699, -0.57146100385840204, 1.6986603412216601, 0.20698230514139301, 5.5487259477006203, 4.3923373342920504, 2.6422104099671002, 0.12108516249854701, 0.00082691786224778699, 0.113024111733347, 2.7966391201096901, -0.052696331353419397, -4.8718908471754601, -4.3479524721016301, 3.1326652654405498, -2.1520186098021998, -17.7293360378548, -0.043482158450766098, 0.29
// SM+spl
double Slater_inv3_2[441] = {-0.056658742514349103, -9.5006825962797503, -18.307026635935799, -0.36744251865365501, 31.610841135834299, 39.084671241308698, -19.703711132730799, 18.666497646833299, 121.964198363159, -17.598444551917201, -2.9520245875548698, 24.935160344122199, 13.992695020707, 4.6612706754314903, -3.6234738869455199, -28.604054847833801, 9.5739197944338503, -62.128443913882499, -67.232459151277695, -11.803301377693099, -3.8508731471217099, 0.069411805979426894, 10.8007994955779, 22.458394727678201, -0.89282708659453403, -39.476164388726097, -47.5420201141931, 24.536662733275399, -23.351558477080602, -146.60044409871099, 21.549928185153298, 3.32891560872004, -30.5549979706051, -17.384760133954, -5.4401084934617003, 4.3730845094487396, 35.2040348414587, -11.414817838980399, 74.926509290741393, 80.145978783950099, 14.6842509146363, 5.7893824160342398, 0.041199702875412497, -6.4779318938824302, -13.2989719358296, 0.51508213587065199, 23.9434304832055, 27.923633944040301, -14.994539048994801, 13.545190256632999, 84.658750799602501, -12.440996315879, -1.6347815496642299, 17.4053340515119, 10.4563880132908, 3.2716511123291698, -2.32038524814726, -20.901080430435499, 6.5092932119827198, -43.310144838525801, -46.266317618221201, -9.0209731319320596, -3.3217327746699499, -0.071364036152653701, 11.529909516191401, 22.9175425489932, 0.44722818276140403, -39.897585842875301, -48.876557451246697, 24.801523764393401, -23.690956189203401, -151.28132599138999, 21.549304148995901, 3.4094704881876998, -30.687597458642099, -17.659848105023201, -5.7069215817131997, 4.0886157921745996, 35.564950307785701, -12.1803912964186, 77.135406641174299, 83.318689013925706, 14.819424258558801, 4.7782939163062803, -0.040051226770134098, -6.7134150132482402, -13.312257348225099, -0.26126985985289403, 22.204441530888801, 27.729627986407898, -13.7012421876807, 13.0503428483308, 85.297842345116607, -12.439369240025201, -1.60848061970372, 17.698309764178902, 10.1798985121263, 3.3110252774145601, -2.5466287558690199, -20.410023365852201, 6.7633439222022202, -43.4149505982517, -46.886928763683102, -7.88950712988734, -2.6927901019362199, -0.058248999176268398, 9.0508747193022892, 18.587941141225699, -0.73447131422595202, -35.308936282117699, -39.940782969048001, 21.9119312895712, -20.076960036338601, -123.935102576568, 17.592669318251101, 3.0259816420954402, -24.782453059479, -15.379073099754899, -4.67601503080067, 3.33242427220405, 30.0344780145867, -9.0772191462285008, 63.391820359557002, 67.760837774760006, 13.3701295015636, 4.8413052345891403, -0.00048534604109837199, -1.8238859168577699, 0.56645976558050504, -0.011652978445607001, 3.3511865791586, 2.0262681174419299, -2.14462816857628, 1.8811114733258301, 1.72516289469276, -0.0410444803602956, -0.17992011162600399, 4.6474809764852196, 0.32638302042170703, 1.28068103014527, 0.42381953652975901, -1.62158383882094, 1.84081037263479, -1.2636267813034701, -0.97726719185308597, -0.84319863660797401, 0.57245121780013697, 0.00072965273109526496, -0.74868100456805098, 0.61375144696985096, -0.013670796650539699, -10.880212250529199, -2.2389072194253701, 6.8446360871752603, -3.8338819099152301, -0.83708262317409499, 0.055638693619984997, 0.16395892268443801, -6.4020540512225601, -0.87619791103918299, 0.67273193621704896, -0.96217527481654197, 3.5322304037025498, 0.74780382069687101, 0.457651060635482, 2.4127372990709501, 5.8207812126355103, 0.41153697505715597, 0.00047300381075295099, -0.24273841428794099, 0.62931133864105804, -0.0075433586410305498, -3.8733549821913198, -2.7042122832831499, 2.4066006453807001, -1.90501007795962, -11.4961539239149, -0.025181201641384401, 0.20951155742128399, 4.41315756942062, -2.02827911285326, 0.53245648205106699, -0.57146100385840204, 1.6986603412216601, 0.20698230514139301, 5.5487259477006203, 4.3923373342920504, 2.6422104099671002, 0.12108516249854701, 0.00082691786224778699, 0.113024111733347, 2.7966391201096901, -0.052696331353419397, -4.8718908471754601, -4.3479524721016301, 3.1326652654405498, -2.1520186098021998, -17.7293360378548, -0.043482158450766098, 0.29
const uint64_t N_updates5 = 5;
const uint64_t Updates_index5[5] = {17, 18, 19, 20, 21};
const double Updates5[105] = {-0.026680206879973502, 0.27298007905483301, -0.099715244024992294, -0.1053539039567112, -0.1449181307107206, 0.086803577840328203, -0.13467331603169441, 0.083233945071697304, -0.050348894670605604, 0.113249283283949, -0.18058512732386639, -0.013084722682833599, -0.2101329043507576, 0.29668186604976698, 0.011996015906333896, -0.25581711530685397, -0.18468007631599939, -0.19795016199350388, 0.13754389435052872, 0.066371031105519007, 0.1506568714976311, 0.1407152302563191, -0.035490237176418013, -0.15555772557854669, -0.13084962218999829, -0.065377140417695004, -0.28165902942419097, -0.0038048550486564081, -0.17355462908744812, 0.1904655769467351, 0.1015159860253334, -0.082969114184379605, 0.1009396240115165, 0.25867408514022822, -0.090622015297412997, -0.0261126160621643, -0.0070611685514450073, 0.34575527906417902, 0.1856994293630127, 0.1505507901310916, -0.24748291820287749, -0.1930690407752986, -0.10824421048164369, -0.2003080397844319, -0.050228431820869016, 0.2341227037832137, -0.1714876033365724, 0.32784904539585158, 0.040038149803876807, 0.12709141941741103, -0.088602993637323102, -0.1617441214621064, -0.237361330538988, -0.029239896684885004, 0.039443209767341697, -0.15546871721744498, -0.0127334967255592, -0.093294814229011994, -0.16617536172270789, 0.081297293305397006, -0.11561803519725759, -0.088157065212726496, 0.043941155076026403, -0.043334499001503005, -0.072291933000086989, 0.25947313755750623, -0.00654759816825385, -0.23022028384730242, -0.0992842186242342, -0.059074921417050084, 0.009816635400056898, -0.018885548226535299, 0.28729220479726802, 0.30690228054299934, 0.058271216228604303, 0.050839929841458791, -0.091462507843971003, 0.21059621125459649, -0.19497598148882408, 0.012118805199861502, -0.040971436537802185, -0.0676504261791709, -0.16150601254776081, 0.076545581221580811, -0.44388696737587502, 0.026131823658942982, -0.11532356310635784, 0.038358195684850299, 0.16635001881513747, 0.035624274984002099, 0.039094586856663179, 0.12423532153479717, 0.13070272840559499, 0.027002811431885002, -0.22131750034168385, 0.1172939799726007, -0.034513717517256695, 0.19796614628285178, 0.087915631011128009, 0.1338309701532128, -0.086308442056179102, -0.014136113226413699, 0.06877816375345, 0.11613245680928279, -0.16100704949349121};
const double Slater5[441] = {-2.8945870399475102, 3.5455725193023699, 2.0470812320709202, -3.5464441776275599, -2.0474903583526598, -2.89596366882324, -0.61329728364944502, 0.70991641283035301, 0.45664468407630898, 0.59523195028305098, 0.26079276204109197, -0.027727147564291999, -0.35093706846237199, -0.095610238611698206, -0.130077064037323, 0.10946778208017301, 0.021800471469759899, 0.048480678349733401, -0.092234551906585693, -0.0160505045205355, -0.42026260495185902, -0.78021097183227495, -0.95558542013168302, -0.55174458026885997, -0.95579683780670199, -0.55179446935653698, 0.78044348955154397, -0.115299895405769, -0.11754634231329, -0.0318448171019554, -0.00082233443390577999, 0.064001239836216001, 0.12740932404995001, -0.049412809312343597, -0.074464745819568606, 0.19851659238338501, -0.088878795504569993, 0.135610401630402, -0.13736967742443101, -0.101879440248013, 0.102390937507153, -0.12476679682731601, -0.00527230883017182, 0.00029381641070358499, 0.0072919740341603799, -0.00026352208806201799, 0.0068838247098028703, 0.0047294595278799499, 0.11728889495134399, -0.062190085649490398, -0.15969239175319699, -0.10645916312933, 0.112943567335606, -0.0015487050404772199, -0.045737843960523598, 0.13196058571338701, 0.038660705089569099, -0.050266433507204097, -0.13745591044426, -0.037740666419267703, 0.11781705915927899, 0.19417031109332999, -0.0079971449449658394, -0.30874741077423101, -0.37815779447555498, 0.218254208564758, 0.378139227628708, -0.21839827299118, -0.30879178643226601, 0.025539221242070202, 0.112043909728527, 0.037676706910133403, 0.025347150862216901, -0.19991625845432301, 0.14163509011268599, 0.13326919078826899, 0.213842943310738, 0.131471157073975, 0.14626120030879999, -0.0118067460134625, 0.093547157943248693, 0.22439678013324699, 0.0082425400614738499, 0.0250354297459126, -1.52881526947021, -1.9163775505148799e-05, -2.1623263359069802, -1.5540548702119799e-05, 2.16284847259521, -1.52953028678894, -0.26402890682220498, 0.024944067001342801, -0.37140902876853898, -0.044071510434150703, -0.29843890666961698, -0.043834518641233403, -0.177285701036453, -0.058585006743669503, -0.0183692276477814, -0.026075478643178902, -0.15623773634433699, -0.011319605633616401, 0.054057534784078598, -0.23151709139347099, -0.0016308756312355399, -0.067333526909351293, 0.00043256155913695698, 0.095035225152969402, -0.00043632232700474598, 0.094624765217304202, 0.066768139600753798, 0.14448715746402699, -0.12090456485748299, -0.203432962298393, -0.21472349762916601, 0.154288679361343, 0.017385326325893399, -0.071518480777740506, 0.27918133139610302, 0.080454014241695404, -0.0361245274543762, 0.0566458702087402, -0.030157707631588, 0.25150132179260298, -0.082445412874221802, -0.017600754275918, -2.2779068946838401, -2.3393469746224599e-05, -3.2218937873840301, -1.7939200915861899e-05, 3.22278833389282, -2.2791447639465301, -0.47945570945739702, 0.045197278261184699, -0.64478874206543002, -0.080234304070472703, -0.47389149665832497, 0.024908870458602898, -0.23455648124218001, -0.10695217549800901, -0.033054649829864502, -0.195749551057816, -0.040076982229948002, 0.094596333801746396, 0.098401188850402804, -0.059563819319009802, -0.0045592403039336196, -0.0025739683769643298, -0.0020174346864223502, -0.0024467820767313199, -0.00204527890309691, 0.00013988610589876801, 0.00073755031917244196, 0.17339251935482, 0.099962860345840496, 0.14622049033641801, -0.088966980576515198, 0.040044382214546197, -0.14437201619148299, -1.5416861060657501e-05, -0.081159926950931494, -0.11909016221761699, 0.170753449201584, 0.031124599277973199, -0.052109345793724102, 0.121445283293724, 0.0289684906601906, 0.120766848325729, -0.0043673445470631097, -0.00014539182302542001, 0.0061219600029289696, 0.000189467871678062, 0.0060523007996380303, 0.00425094133242965, 0.0905451029539108, 0.058540407568216303, -0.15191346406936601, 0.11511342227459, 0.15016922354698201, 0.10024280846118901, -0.11066282540559801, -0.15089659392833699, -0.021995550021529201, 0.10021472722291901, 0.019932290539145501, 0.070281185209751101, -0.120184391736984, -0.03055
2021-09-18 18:07:05 +02:00
double Slater_inv5[441] = {-0.054189244668834902, -105.426713929607, -88.458496476283003, 1.5333775291907901, -306.17152423250297, 211.834723659242, 189.348181800731, -164.85602878397, 1001.02895803198, -19.086531171244101, -14.1862411100869, 93.929906236367501, -177.880045125896, -6.7382862272631598, -14.511503830974201, 198.86948709278099, 116.24375034946701, -509.187787693936, -474.82187536023201, 185.49468275501101, -68.704359475869197, 0.066396729468773799, 128.61396390514599, 107.279152808508, -3.2214077352586501, 377.39638500462598, -258.32727230254102, -233.50899542599601, 202.73053181330999, -1222.6711957145401, 23.374193247961198, 17.0231928902612, -115.31992914359, 218.34781344598201, 8.5868800406738099, 17.706454078785399, -244.13052330624399, -142.42769283244499, 622.15335091370196, 579.15535219316803, -228.72024346964599, 85.262861528059901, 0.042977214694503003, -75.818887039648899, -63.4429496117761, 1.89037149893251, -221.440571714557, 152.54192851849001, 136.92110635753301, -119.073716594494, 720.98982202338402, -13.516656591400601, -9.7291496308540406, 67.372278162472895, -128.19425354399101, -4.9638472213270699, -10.200417492010599, 143.31769471617801, 83.615852083439094, -366.86654819753397, -341.426964056759, 133.98241248588599, -50.230799667840699, -0.074440153582489205, 131.003888225851, 110.340015781723, -1.92534262087524, 381.482058581253, -264.51148724341499, -236.035513173862, 204.43864191351099, -1249.3967405711001, 23.404380739562701, 17.427737318703301, -116.737852869633, 221.189913802663, 8.4677080090305399, 17.684126766935599, -247.44394443997601, -145.031902348194, 635.53348480709997, 592.58164589603302, -230.87850923206801, 85.694308069456397, -0.038299262589364599, -74.702005824131604, -63.171412849899497, 1.0859859506457199, -217.148397050222, 150.25957665606299, 134.431598757773, -117.013550651343, 708.34335491195498, -13.4940069422526, -9.5788103819957797, 66.593490547156307, -125.80391712960601, -4.7671091217963504, -10.2632044954116, 140.809682960885, 82.366527572926103, -360.28206877713399, -335.74815599116999, 131.925760239261, -48.660685665615503, -0.060749768867194, 106.92546569584999, 88.790343162838894, -2.6692288353747702, 312.16097467403699, -214.903931604847, -193.215823966269, 167.781122174079, -1018.10641886182, 19.107519450141002, 14.388719322765001, -95.267373183568594, 180.40857560667001, 6.9728512293157801, 14.4102075774752, -201.83015329787199, -117.917199027973, 518.09179195913805, 482.46342773412402, -189.25708821193001, 70.869624357767506, -0.000187411047696029, -11.9846117540694, -9.5726690972884096, 0.23970985157893901, -31.525734707294099, 26.424001662617101, 19.629541977294, -12.8111861333484, 126.186955595279, -0.21994779537689901, -1.75278408990554, 12.699562562401301, -19.050321348645301, 0.28466729039388, -1.07057307471916, 20.675591706275501, 13.1048598424623, -64.121267975873096, -59.394985304548896, 16.918592992744401, -7.70362909966763, 0.00036065684724530402, 12.9068309168011, 12.002712393435401, -0.33194601220397202, 40.326484701927001, -31.878253347415502, -25.170560077575001, 17.9905011861023, -156.05779622500401, 0.28522581276215098, 2.06069522309385, -17.151275363286, 25.397351917205899, 2.08581286823873, 0.91118057859336898, -26.515090998756101, -14.4035822159035, 78.890023885294596, 75.240667711773398, -20.769502289340899, 10.973365712218101, 0.00054214253230406803, -2.29504518628212, -2.1312244512866698, 0.047562978801466302, -10.6318457557031, 3.0591772010668401, 6.6375862982055702, -4.4363733110097296, 16.291945106274301, -0.062929125223047597, -0.16835581689436099, 6.0820271877694996, -5.8351950377758204, 0.35401259185447298, -0.90070453596893996, 6.0382685793059503, 2.4785619931575602, -8.4908203096171402, -8.6197966065449307, 5.8804527059182696, -1.67757738078734, 0.00096800508899916803, -2.4934502646428101, -4.5195305250493698, 0.052900907147883501, -6.5403147654399101, 8.5375829900159896, 4.1350660479474399, -1.63386715634978, 39.5635438458416, -0.11054656518749099, -0.56905231469759199, 10.2311793759502, -8.8014973731182806, -1.4221
2021-09-05 11:22:41 +02:00
2021-09-18 18:07:05 +02:00
assert(Updates1 != NULL);
assert(Updates_index1 != NULL);
assert(Slater_inv1 != NULL);
2021-10-13 10:31:24 +02:00
// original determinant of Slater1 (before applying updates)
double det = 3.407025646103221e-10;
rc = qmckl_sherman_morrison(context, Dim, N_updates1, Updates1, Updates_index1, breakdown, Slater_inv1, &det);
// Check that the determinant is updated properly
2021-10-13 10:45:54 +02:00
assert(fabs(det + 4.120398385068217e-10) < 1e-15);
2021-10-13 10:31:24 +02:00
2021-09-05 11:22:41 +02:00
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 * Dim + 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;
}
}
}
2021-07-20 19:34:51 +02:00
assert(rc == QMCKL_SUCCESS);
#+end_src
2021-07-21 17:30:12 +02:00
* Woodbury 2x2
2021-07-21 17:42:48 +02:00
2021-07-21 17:30:12 +02:00
** ~qmckl_woodbury_2~
:PROPERTIES:
:Name: qmckl_woodbury_2
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-09-07 11:22:54 +02:00
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
2021-09-01 16:06:51 +02:00
\[
(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
2021-10-06 23:44:06 +02:00
2021-10-13 10:45:54 +02:00
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.
2021-07-21 17:30:12 +02:00
#+NAME: qmckl_woodbury_2_args
2021-07-21 17:42:48 +02:00
| 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 |
2021-07-26 17:41:21 +02:00
| double | breakdown | in | Break-down parameter on which to fail or not |
2021-07-21 17:42:48 +02:00
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
2021-10-13 10:45:54 +02:00
| double* | determinant | inout | Determinant of Slater-matrix |
2021-07-21 17:30:12 +02:00
*** Requirements
2021-09-01 16:06:51 +02:00
- ~context~ is not ~qmckl_null_context~
2021-09-07 11:22:54 +02:00
- ~Dim >= 2~
- ~Updates~ is allocated with $2 \times Dim$ elements
- ~Updates_index~ is allocated with $2$ elements
2021-09-01 16:06:51 +02:00
- ~breakdown~ is a small number such that $0 < breakdown << 1$
2021-09-07 11:22:54 +02:00
- ~Slater_inv~ is allocated with $Dim \times Dim$ elements
2021-07-21 17:30:12 +02:00
*** C header
#+CALL: generate_c_header(table=qmckl_woodbury_2_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
2021-09-16 13:44:12 +02:00
qmckl_exit_code qmckl_woodbury_2(
2021-07-21 17:30:12 +02:00
const qmckl_context context,
2021-10-09 22:23:12 +02:00
const uint64_t Dim,
2021-07-21 17:30:12 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 10:45:54 +02:00
double* Slater_inv,
2021-10-13 12:03:18 +02:00
double* determinant);
2021-07-21 17:30:12 +02:00
#+end_src
2021-10-09 22:23:12 +02:00
*** C source
2021-07-21 17:30:12 +02:00
#+begin_src c :tangle (eval c) :comments org
#include <stdbool.h >
2021-09-10 17:25:56 +02:00
#include <math.h >
2021-07-21 17:30:12 +02:00
#include "qmckl.h"
2021-10-09 22:23:12 +02:00
qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
const uint64_t Dim,
2021-07-21 17:30:12 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 10:45:54 +02:00
double* Slater_inv,
double* determinant) {
2021-07-21 17:30:12 +02:00
/*
C := S^{-1} * U, dim x 2
B := 1 + V * C, 2 x 2
D := V * S^{-1}, 2 x dim
2021-09-10 17:25:56 +02:00
,*/
2021-09-16 13:44:12 +02:00
2021-10-06 23:44:06 +02:00
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
2021-07-22 09:59:02 +02:00
const uint64_t row1 = (Updates_index[0] - 1);
const uint64_t row2 = (Updates_index[1] - 1);
2021-07-21 17:30:12 +02:00
// Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE
// OF LAYOUT OF 'Updates' !!
double C[2 * Dim];
2021-07-22 09:59:02 +02:00
for (uint64_t i = 0; i < Dim; i++) {
for (uint64_t j = 0; j < 2; j++) {
2021-07-21 17:30:12 +02:00
C[i * 2 + j] = 0;
2021-07-22 09:59:02 +02:00
for (uint64_t k = 0; k < Dim; k++) {
2021-07-21 17:30:12 +02:00
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;
2021-07-26 17:41:21 +02:00
if (fabs(det) < breakdown) {
2021-07-21 17:30:12 +02:00
return QMCKL_FAILURE;
}
2021-10-13 10:45:54 +02:00
// Update det(S) when passed
if (determinant != NULL)
*determinant *= det;
2021-07-21 17:30:12 +02:00
// 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];
2021-07-22 09:59:02 +02:00
for (uint64_t i = 0; i < 2; i++) {
for (uint64_t j = 0; j < Dim; j++) {
2021-07-21 17:30:12 +02:00
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
2021-07-22 09:59:02 +02:00
for (uint64_t i = 0; i < Dim; i++) {
for (uint64_t j = 0; j < Dim; j++) {
2021-07-21 17:30:12 +02:00
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;
}
#+end_src
*** Performance
2021-09-07 12:22:39 +02:00
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.
2021-09-10 17:25:56 +02:00
2021-10-09 22:23:12 +02:00
** Fortran interface :noexport:
2021-09-14 12:13:17 +02:00
:PROPERTIES:
:Name: qmckl_woodbury_2
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-07-21 17:30:12 +02:00
#+CALL: generate_f_interface(table=qmckl_woodbury_2_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
2021-09-16 13:44:12 +02:00
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_woodbury_2 &
2021-10-13 10:45:54 +02:00
(context, Dim, Updates, Updates_index, breakdown, Slater_inv, determinant) &
2021-09-16 13:44:12 +02:00
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
2021-10-09 22:23:12 +02:00
integer (c_int64_t) , intent(in) , value :: Dim
2021-09-16 13:44:12 +02:00
real (c_double ) , intent(in) :: Updates(2*Dim)
integer (c_int64_t) , intent(in) :: Updates_index(2)
2021-10-09 22:23:12 +02:00
real (c_double ) , intent(in) , value :: breakdown
2021-09-16 13:44:12 +02:00
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
2021-10-13 10:45:54 +02:00
real (c_double ) , intent(inout) :: determinant
2021-09-16 13:44:12 +02:00
end function qmckl_woodbury_2
end interface
#+end_src
2021-07-21 17:30:12 +02:00
*** Test :noexport:
#+begin_src c :tangle (eval c_test)
2021-09-18 18:07:05 +02:00
assert(Updates2 != NULL);
assert(Updates_index2 != NULL);
assert(Slater_inv2 != NULL);
2021-10-13 10:45:54 +02:00
det = -1.4432116661319376e-11;
rc = qmckl_woodbury_2(context, Dim, Updates2, Updates_index2, breakdown, Slater_inv2, &det);
assert(fabs(det-2.367058141251457e-10) < 1e-15);
2021-09-05 11:22:41 +02:00
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] + = Slater2[i * Dim + k] * Slater_inv2[k * Dim + 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;
}
}
}
2021-07-21 17:30:12 +02:00
assert(rc == QMCKL_SUCCESS);
#+end_src
2021-07-22 11:38:50 +02:00
* Woodbury 3x3
2021-10-06 23:44:06 +02:00
2021-07-22 11:38:50 +02:00
** ~qmckl_woodbury_3~
:PROPERTIES:
:Name: qmckl_woodbury_3
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-09-01 16:06:51 +02:00
The 3x3 version of the Woodbury 2x2 kernel. It is used to apply three
2021-09-02 17:56:42 +02:00
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:
2021-09-01 16:06:51 +02:00
$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
2021-07-22 11:38:50 +02:00
2021-10-13 10:51:42 +02:00
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.
2021-07-22 11:38:50 +02:00
#+NAME: qmckl_woodbury_3_args
| 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 |
2021-07-26 17:41:21 +02:00
| double | breakdown | in | Break-down parameter on which to fail or not |
2021-07-22 11:38:50 +02:00
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
2021-10-13 10:51:42 +02:00
| double* | determinant | inout | Determinant of Slater-matrix |
2021-07-22 11:38:50 +02:00
*** Requirements
2021-09-01 16:10:59 +02:00
- ~context~ is not ~qmckl_null_context~
2021-09-07 11:22:54 +02:00
- ~Dim >= 2~
- ~Updates~ is allocated with $3 \times Dim$ elements
- ~Updates_index~ is allocated with $3$ elements
2021-09-01 16:10:59 +02:00
- ~breakdown~ is a small number such that $0 < breakdown << 1$
2021-09-07 11:22:54 +02:00
- ~Slater_inv~ is allocated with $Dim \times Dim$ elements
2021-07-22 11:38:50 +02:00
*** C header
#+CALL: generate_c_header(table=qmckl_woodbury_3_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
2021-09-16 13:44:12 +02:00
qmckl_exit_code qmckl_woodbury_3(
2021-07-22 11:38:50 +02:00
const qmckl_context context,
2021-10-09 22:23:12 +02:00
const uint64_t Dim,
2021-07-22 11:38:50 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 10:51:42 +02:00
double* Slater_inv,
double* determinant);
2021-07-22 11:38:50 +02:00
#+end_src
2021-10-09 22:23:12 +02:00
*** C source
2021-07-22 11:38:50 +02:00
#+begin_src c :tangle (eval c) :comments org
#include <stdbool.h >
2021-09-10 17:25:56 +02:00
#include <math.h >
2021-07-22 11:38:50 +02:00
#include "qmckl.h"
2021-10-09 22:23:12 +02:00
qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
const uint64_t Dim,
2021-07-22 11:38:50 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 10:51:42 +02:00
double* Slater_inv,
double* determinant) {
2021-07-22 11:38:50 +02:00
/*
C := S^{-1} * U, dim x 3
B := 1 + V * C, 3 x 3
D := V * S^{-1}, 3 x dim
,*/
2021-09-16 13:44:12 +02:00
2021-10-06 23:44:06 +02:00
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
2021-07-22 11:38:50 +02:00
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);
2021-07-26 17:41:21 +02:00
if (fabs(det) < breakdown) {
2021-07-22 11:38:50 +02:00
return QMCKL_FAILURE;
}
2021-10-13 10:51:42 +02:00
// Update det(Slater) if passed
if (determinant != NULL)
*determinant *= det;
2021-07-22 11:38:50 +02:00
// 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;
}
#+end_src
*** Performance...
2021-10-06 23:44:06 +02:00
2021-09-07 12:22:39 +02:00
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.
2021-09-01 16:06:51 +02:00
2021-10-09 22:23:12 +02:00
** Fortran interface :noexport:
2021-09-14 12:13:17 +02:00
:PROPERTIES:
:Name: qmckl_woodbury_3
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-07-22 11:38:50 +02:00
#+CALL: generate_f_interface(table=qmckl_woodbury_3_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
2021-09-16 13:44:12 +02:00
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_woodbury_3 &
2021-10-13 10:51:42 +02:00
(context, Dim, Updates, Updates_index, breakdown, Slater_inv, determinant) &
2021-09-16 13:44:12 +02:00
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
2021-10-09 22:23:12 +02:00
integer (c_int64_t) , intent(in) , value :: Dim
2021-09-16 13:44:12 +02:00
real (c_double ) , intent(in) :: Updates(3*Dim)
integer (c_int64_t) , intent(in) :: Updates_index(3)
2021-10-09 22:23:12 +02:00
real (c_double ) , intent(in) , value :: breakdown
2021-09-16 13:44:12 +02:00
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
2021-10-13 10:51:42 +02:00
real (c_double ) , intent(inout) :: determinant
2021-09-16 13:44:12 +02:00
end function qmckl_woodbury_3
end interface
#+end_src
2021-07-22 11:38:50 +02:00
*** Test :noexport:
#+begin_src c :tangle (eval c_test)
2021-09-18 18:07:05 +02:00
assert(Updates3 != NULL);
assert(Updates_index3 != NULL);
assert(Slater_inv3_1 != NULL);
2021-10-13 10:51:42 +02:00
det = -1.23743195512859e-09;
rc = qmckl_woodbury_3(context, Dim, Updates3, Updates_index3, breakdown, Slater_inv3_1, &det);
assert(fabs(det - 1.602708950725074e-10) < 1e-15);
2021-09-05 11:22:41 +02:00
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] + = Slater3[i * Dim + k] * Slater_inv3_1[k * Dim + 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;
}
}
}
2021-07-22 11:38:50 +02:00
assert(rc == QMCKL_SUCCESS);
#+end_src
2021-07-22 18:05:39 +02:00
* Sherman-Morrison with update splitting
2021-10-06 23:44:06 +02:00
2021-07-22 18:05:39 +02:00
** ~qmckl_sherman_morrison_splitting~
:PROPERTIES:
:Name: qmckl_sherman_morrison_splitting
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-09-02 17:56:42 +02:00
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:
2021-09-07 12:22:39 +02:00
$u_j \rightarrow \frac{1}{2} u_j$. One half is applied immediately --necessarily increasing the value of the
2021-09-02 17:56:42 +02:00
denominator because of the split-- while the other halve is put in a queue that will be applied when all the
2021-09-07 12:22:39 +02:00
remaining updates have been treated.
The kernel is executed recursively until the queue is eiter empty and all
2021-09-02 17:56:42 +02:00
updates are applied successfully, or the size of the queue equals the number of initial updates. In the last
2021-09-07 12:22:39 +02:00
case the Slater-matrix that would have resulted from applying the updates is singular and therefore the
2021-09-02 17:56:42 +02:00
kernel exits with an exit code.
2021-07-22 18:05:39 +02:00
2021-10-13 11:55:20 +02:00
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.
2021-07-22 18:05:39 +02:00
#+NAME: qmckl_sherman_morrison_splitting_args
2021-07-22 18:20:20 +02:00
| 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 |
2021-07-26 17:41:21 +02:00
| double | breakdown | in | Break-down parameter on which to fail or not |
2021-07-22 18:20:20 +02:00
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
2021-10-13 11:55:20 +02:00
| 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.
2021-07-22 18:05:39 +02:00
*** Requirements
2021-09-07 11:22:54 +02:00
* ~context~ is not ~QMCKL_NULL_CONTEXT~
* ~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
2021-10-06 23:44:06 +02:00
2021-07-22 18:05:39 +02:00
*** C header
#+CALL: generate_c_header(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
2021-07-22 18:20:20 +02:00
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
2021-09-16 13:44:12 +02:00
qmckl_exit_code qmckl_sherman_morrison_splitting(
2021-07-22 18:20:20 +02:00
const qmckl_context context,
2021-10-09 22:23:12 +02:00
const uint64_t Dim,
const uint64_t N_updates,
2021-07-22 18:20:20 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 11:55:20 +02:00
double* Slater_inv,
2021-10-13 12:03:18 +02:00
double* determinant);
2021-07-22 18:20:20 +02:00
#+end_src
2021-10-09 22:23:12 +02:00
*** C source
2021-07-22 18:05:39 +02:00
#+begin_src c :tangle (eval c) :comments org
#include <stdbool.h >
#include "qmckl.h"
2021-10-09 22:23:12 +02:00
qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
2021-07-22 18:05:39 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 11:55:20 +02:00
double* Slater_inv,
double* determinant) {
2021-10-06 23:44:06 +02:00
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
2021-10-09 22:23:12 +02:00
double later_updates[Dim * N_updates];
uint64_t later_index[N_updates];
2021-07-22 18:05:39 +02:00
uint64_t later = 0;
2021-10-09 22:23:12 +02:00
(void) qmckl_slagel_splitting(Dim, N_updates, Updates, Updates_index,
2021-10-13 11:55:20 +02:00
breakdown, Slater_inv, later_updates, later_index, &later, determinant);
2021-10-06 23:44:06 +02:00
2021-07-22 18:05:39 +02:00
if (later > 0) {
2021-10-09 22:23:12 +02:00
(void) qmckl_sherman_morrison_splitting(context, Dim, later,
2021-10-13 11:55:20 +02:00
later_updates, later_index, breakdown, Slater_inv, determinant);
2021-07-22 18:05:39 +02:00
}
return QMCKL_SUCCESS;
}
#+end_src
*** Performance...
2021-09-07 12:22:39 +02:00
This kernel performs best when there are 2 or more rank-1 update cycles and fail-rate is high.
2021-10-06 23:44:06 +02:00
2021-10-09 22:23:12 +02:00
** Fortran interface :noexport:
2021-09-14 12:13:17 +02:00
:PROPERTIES:
:Name: qmckl_sherman_morrison_splitting
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-07-22 18:05:39 +02:00
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
2021-09-16 17:24:19 +02:00
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison_splitting &
2021-10-13 11:55:20 +02:00
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
2021-09-16 17:24:19 +02:00
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
2021-10-09 22:23:12 +02:00
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
2021-09-16 17:24:19 +02:00
real (c_double ) , intent(in) :: Updates(N_updates*Dim)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
2021-10-09 22:23:12 +02:00
real (c_double ) , intent(in) , value :: breakdown
2021-09-16 17:24:19 +02:00
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
2021-10-13 11:55:20 +02:00
real (c_double ) , intent(inout) :: determinant
2021-09-16 17:24:19 +02:00
end function qmckl_sherman_morrison_splitting
end interface
#+end_src
2021-07-22 18:05:39 +02:00
*** Test :noexport:
#+begin_src c :tangle (eval c_test)
2021-09-18 18:07:05 +02:00
assert(Updates3 != NULL);
assert(Updates_index3 != NULL);
assert(Slater_inv3_2 != NULL);
2021-10-13 11:55:20 +02:00
det = -1.23743195512859e-09;
rc = qmckl_sherman_morrison_splitting(context, Dim, N_updates3, Updates3, Updates_index3, breakdown, Slater_inv3_2, &det);
assert(fabs(det - 1.602708950725074e-10) < 1e-15);
2021-09-05 11:22:41 +02:00
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] + = Slater3[i * Dim + k] * Slater_inv3_2[k * Dim + 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;
}
}
}
2021-07-22 18:05:39 +02:00
assert(rc == QMCKL_SUCCESS);
#+end_src
2021-07-26 17:41:21 +02:00
* Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting
2021-10-06 23:44:06 +02:00
2021-07-23 16:50:51 +02:00
** ~qmckl_sherman_morrison_smw32s~
:PROPERTIES:
:Name: qmckl_sherman_morrison_smw32s
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-09-02 17:56:42 +02:00
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.
2021-07-23 16:50:51 +02:00
2021-10-13 11:55:20 +02:00
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.
2021-07-23 16:50:51 +02:00
#+NAME: qmckl_sherman_morrison_smw32s_args
| 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 |
2021-07-26 17:41:21 +02:00
| double | breakdown | in | Break-down parameter on which to fail or not |
2021-07-23 16:50:51 +02:00
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
2021-10-13 11:55:20 +02:00
| double* | determinant | inout | Determinant of the Slater-matrix |
2021-07-23 16:50:51 +02:00
*** Requirements
2021-09-07 11:22:54 +02:00
* ~context~ is not ~QMCKL_NULL_CONTEXT~
* ~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
2021-07-23 16:50:51 +02:00
*** C header
#+CALL: generate_c_header(table=qmckl_sherman_morrison_smw32s_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
2021-07-23 19:58:24 +02:00
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
2021-09-16 17:24:19 +02:00
qmckl_exit_code qmckl_sherman_morrison_smw32s(
2021-07-23 19:58:24 +02:00
const qmckl_context context,
2021-10-09 22:23:12 +02:00
const uint64_t Dim,
const uint64_t N_updates,
2021-07-23 19:58:24 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 11:55:20 +02:00
double* Slater_inv,
2021-10-13 12:03:18 +02:00
double* determinant);
2021-07-23 19:58:24 +02:00
#+end_src
2021-10-09 22:23:12 +02:00
*** C source
2021-07-23 16:50:51 +02:00
#+begin_src c :tangle (eval c) :comments org
#include <stdbool.h >
#include "qmckl.h"
2021-10-09 22:23:12 +02:00
qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
2021-07-23 16:50:51 +02:00
const double* Updates,
const uint64_t* Updates_index,
2021-10-09 22:23:12 +02:00
const double breakdown,
2021-10-13 11:55:20 +02:00
double* Slater_inv,
double* determinant) {
2021-07-23 16:50:51 +02:00
2021-10-06 23:44:06 +02:00
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
2021-07-23 16:50:51 +02:00
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++) {
2021-07-27 08:48:28 +02:00
const double *Updates_3block = &Updates[i * length_3block];
const uint64_t *Updates_index_3block = &Updates_index[i * 3];
2021-10-13 11:55:20 +02:00
rc = qmckl_woodbury_3(context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, determinant);
2021-07-23 16:50:51 +02:00
if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0;
2021-10-09 22:23:12 +02:00
rc = qmckl_slagel_splitting(Dim, 3, Updates_3block, Updates_index_3block,
2021-10-13 11:55:20 +02:00
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
2021-07-23 16:50:51 +02:00
later = later + l;
}
}
}
2021-09-16 17:24:19 +02:00
// Apply last remaining block of 2 updates with Woodbury 2x2 kernel
if (remainder == 2) {
2021-07-27 08:48:28 +02:00
const double *Updates_2block = &Updates[n_of_3blocks * length_3block];
const uint64_t *Updates_index_2block = &Updates_index[3 * n_of_3blocks];
2021-10-13 11:55:20 +02:00
rc = qmckl_woodbury_2(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv, determinant);
2021-07-23 16:50:51 +02:00
if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0;
2021-10-09 22:23:12 +02:00
(void) qmckl_slagel_splitting(Dim, 2, Updates_2block, Updates_index_2block,
2021-10-13 11:55:20 +02:00
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
2021-07-23 16:50:51 +02:00
later = later + l;
}
}
2021-09-16 17:24:19 +02:00
// Apply last remaining update with slagel_splitting
else if (remainder == 1) {
2021-07-27 08:48:28 +02:00
const double *Updates_1block = &Updates[n_of_3blocks * length_3block];
const uint64_t *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
2021-07-23 16:50:51 +02:00
uint64_t l = 0;
2021-10-09 22:23:12 +02:00
(void) qmckl_slagel_splitting(Dim, 1, Updates_1block, Updates_index_1block,
2021-10-13 11:55:20 +02:00
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
2021-07-23 16:50:51 +02:00
later = later + l;
}
if (later > 0) {
2021-10-13 11:55:20 +02:00
(void) qmckl_sherman_morrison_splitting(context, Dim, later, later_updates, later_index, breakdown, Slater_inv, determinant);
2021-07-23 16:50:51 +02:00
}
2021-07-27 06:59:44 +02:00
return QMCKL_SUCCESS;
2021-07-23 16:50:51 +02:00
}
#+end_src
*** Performance...
2021-09-07 12:22:39 +02:00
This kernel performs best for update cycles with 2 or more rank-1 updates and the fail-rate is low.
2021-09-02 17:56:42 +02:00
2021-10-09 22:23:12 +02:00
** Fortran interface :noexport:
2021-09-14 12:13:17 +02:00
:PROPERTIES:
:Name: qmckl_sherman_morrison_smw32s
:CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code
:END:
2021-07-23 16:50:51 +02:00
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_smw32s_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
2021-09-16 17:24:19 +02:00
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sherman_morrison_smw32s &
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
2021-10-09 22:23:12 +02:00
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
2021-09-16 17:24:19 +02:00
real (c_double ) , intent(in) :: Updates(N_updates*Dim)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
2021-10-09 22:23:12 +02:00
real (c_double ) , intent(in) , value :: breakdown
2021-09-16 17:24:19 +02:00
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
end function qmckl_sherman_morrison_smw32s
end interface
#+end_src
2021-07-23 16:50:51 +02:00
*** Test :noexport:
#+begin_src c :tangle (eval c_test)
2021-09-18 18:07:05 +02:00
assert(Updates5 != NULL);
assert(Updates_index5 != NULL);
assert(Slater_inv5 != NULL);
2021-10-13 11:55:20 +02:00
det = -3.186005284713128e-10;
rc = qmckl_sherman_morrison_smw32s(context, Dim, N_updates5, Updates5, Updates_index5, breakdown, Slater_inv5, &det);
assert(fabs(det + 5.260200118412903e-10) < 1e-15);
2021-09-05 11:22:41 +02:00
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++) {
2021-09-18 18:07:05 +02:00
res[i * Dim + j] + = Slater5[i * Dim + k] * Slater_inv5[k * Dim + j];
2021-09-05 11:22:41 +02:00
}
}
}
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;
}
}
}
2021-07-23 16:50:51 +02:00
assert(rc == QMCKL_SUCCESS);
#+end_src
2021-10-06 23:44:06 +02:00
2021-09-07 09:27:22 +02:00
* Helper Functions
2021-10-06 23:44:06 +02:00
Private helper-functions that are used by the Sherman-Morrison-Woodbury kernels.
2021-09-07 11:22:54 +02:00
These functions can only be used internally by the kernels in this module.
2021-09-07 09:27:22 +02:00
** ~qmckl_slagel_splitting~
:PROPERTIES:
:Name: qmckl_slagel_splitting
:CRetType: double
:FRetType: double precision
:END:
2021-09-07 11:22:54 +02:00
~qmckl_slagel_splitting~ is the non-recursive, inner part of the 'Sherman-Morrison with update splitting'-kernel.
2021-09-07 12:22:39 +02:00
It is used internally to apply a collection of $N$ rank-1 updates to the inverse Slater-matrix $S^{-1}$ and
2021-09-07 11:22:54 +02:00
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.
2021-10-06 23:44:06 +02:00
2021-09-07 11:22:54 +02:00
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}.
2021-10-06 23:44:06 +02:00
2021-10-13 11:55:20 +02:00
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.
2021-09-07 09:27:22 +02:00
#+NAME: qmckl_slagel_splitting_args
| 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 |
2021-10-13 11:55:20 +02:00
| double* | determinant | inout | Determinant of the Slater-matrix |
2021-10-06 23:44:06 +02:00
2021-09-07 09:27:22 +02:00
*** Requirements
- ~Dim >= 2~
- ~N_updates >= 1~
2021-09-07 11:22:54 +02:00
- ~Updates~ is allocated with $N_updates \times Dim$ elements
- ~Updates_index~ is allocated with $N_updates$ elements
2021-09-07 09:27:22 +02:00
- ~breakdown~ is a small number such that $0 < breakdown << 1$
2021-09-07 11:22:54 +02:00
- ~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
2021-10-06 23:44:06 +02:00
- ~later >= 0~
2021-09-07 09:27:22 +02:00
*** C header
#+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
2021-10-09 22:23:12 +02:00
qmckl_exit_code qmckl_slagel_splitting (
2021-09-07 09:27:22 +02:00
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,
2021-10-13 11:55:20 +02:00
uint64_t* later,
2021-10-13 12:03:18 +02:00
double* determinant);
2021-09-07 09:27:22 +02:00
#+end_src
2021-10-09 22:23:12 +02:00
*** C source
2021-09-07 09:27:22 +02:00
#+begin_src c :tangle (eval c) :comments org
#include <stdbool.h >
#include <math.h >
#include "qmckl.h"
2021-10-09 22:23:12 +02:00
qmckl_exit_code qmckl_slagel_splitting(uint64_t Dim,
2021-09-07 09:27:22 +02:00
uint64_t N_updates,
2021-10-06 23:44:06 +02:00
const double *Updates,
2021-09-07 09:27:22 +02:00
const uint64_t *Updates_index,
const double breakdown,
double *Slater_inv,
2021-10-06 23:44:06 +02:00
double *later_updates,
2021-09-07 09:27:22 +02:00
uint64_t *later_index,
2021-10-13 11:55:20 +02:00
uint64_t *later,
double *determinant) {
2021-09-07 09:27:22 +02:00
// #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];
2021-09-07 12:22:39 +02:00
if (fabs(den) < breakdown) { // Here is decided to split the update, or not.
2021-09-07 09:27:22 +02:00
2021-09-07 12:22:39 +02:00
// U_l = U_l / 2: split the update in 2 equal halves and save the second halve in later_updates
2021-09-07 09:27:22 +02:00
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];
2021-09-07 12:22:39 +02:00
} // From here onwards we continue with applying the first havel of the update to Slater_inv
2021-09-07 09:27:22 +02:00
double iden = 1 / den;
2021-10-13 11:55:20 +02:00
if (determinant != NULL)
*determinant *= den;
2021-09-07 09:27:22 +02:00
// 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;
}
#+end_src
2021-10-06 23:44:06 +02:00
2021-09-07 09:27:22 +02:00
*** Performance
2021-09-07 12:22:39 +02:00
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.
2021-09-07 09:27:22 +02:00
2021-07-19 12:01:07 +02:00
* End of files
#+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0;
2021-10-06 23:44:06 +02:00
2021-07-19 12:01:07 +02:00
}
#+end_src
# -*- mode: org -*-
# vim: syntax=c
