2021-09-10 17:38:14 +02:00
<?xml version="1.0" encoding="utf-8"?>
< !DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
< html xmlns = "http://www.w3.org/1999/xhtml" lang = "en" xml:lang = "en" >
< head >
2022-03-07 23:23:02 +01:00
<!-- 2022 - 03 - 07 Mon 22:22 -->
2021-09-10 17:38:14 +02:00
< meta http-equiv = "Content-Type" content = "text/html;charset=utf-8" / >
< meta name = "viewport" content = "width=device-width, initial-scale=1" / >
< title > Sherman-Morrison-Woodbury< / title >
< meta name = "generator" content = "Org mode" / >
< meta name = "author" content = "TREX CoE" / >
< style type = "text/css" >
<!-- /* --> <![CDATA[/*> <!-- */
.title { text-align: center;
margin-bottom: .2em; }
.subtitle { text-align: center;
font-size: medium;
font-weight: bold;
margin-top:0; }
.todo { font-family: monospace; color: red; }
.done { font-family: monospace; color: green; }
.priority { font-family: monospace; color: orange; }
.tag { background-color: #eee; font-family: monospace;
padding: 2px; font-size: 80%; font-weight: normal; }
.timestamp { color: #bebebe; }
.timestamp-kwd { color: #5f9ea0; }
.org-right { margin-left: auto; margin-right: 0px; text-align: right; }
.org-left { margin-left: 0px; margin-right: auto; text-align: left; }
.org-center { margin-left: auto; margin-right: auto; text-align: center; }
.underline { text-decoration: underline; }
#postamble p, #preamble p { font-size: 90%; margin: .2em; }
p.verse { margin-left: 3%; }
pre {
border: 1px solid #ccc;
box-shadow: 3px 3px 3px #eee;
padding: 8pt;
font-family: monospace;
overflow: auto;
margin: 1.2em;
}
pre.src {
position: relative;
overflow: visible;
padding-top: 1.2em;
}
pre.src:before {
display: none;
position: absolute;
background-color: white;
top: -10px;
right: 10px;
padding: 3px;
border: 1px solid black;
}
pre.src:hover:before { display: inline;}
/* Languages per Org manual */
pre.src-asymptote:before { content: 'Asymptote'; }
pre.src-awk:before { content: 'Awk'; }
pre.src-C:before { content: 'C'; }
/* pre.src-C++ doesn't work in CSS */
pre.src-clojure:before { content: 'Clojure'; }
pre.src-css:before { content: 'CSS'; }
pre.src-D:before { content: 'D'; }
pre.src-ditaa:before { content: 'ditaa'; }
pre.src-dot:before { content: 'Graphviz'; }
pre.src-calc:before { content: 'Emacs Calc'; }
pre.src-emacs-lisp:before { content: 'Emacs Lisp'; }
pre.src-fortran:before { content: 'Fortran'; }
pre.src-gnuplot:before { content: 'gnuplot'; }
pre.src-haskell:before { content: 'Haskell'; }
pre.src-hledger:before { content: 'hledger'; }
pre.src-java:before { content: 'Java'; }
pre.src-js:before { content: 'Javascript'; }
pre.src-latex:before { content: 'LaTeX'; }
pre.src-ledger:before { content: 'Ledger'; }
pre.src-lisp:before { content: 'Lisp'; }
pre.src-lilypond:before { content: 'Lilypond'; }
pre.src-lua:before { content: 'Lua'; }
pre.src-matlab:before { content: 'MATLAB'; }
pre.src-mscgen:before { content: 'Mscgen'; }
pre.src-ocaml:before { content: 'Objective Caml'; }
pre.src-octave:before { content: 'Octave'; }
pre.src-org:before { content: 'Org mode'; }
pre.src-oz:before { content: 'OZ'; }
pre.src-plantuml:before { content: 'Plantuml'; }
pre.src-processing:before { content: 'Processing.js'; }
pre.src-python:before { content: 'Python'; }
pre.src-R:before { content: 'R'; }
pre.src-ruby:before { content: 'Ruby'; }
pre.src-sass:before { content: 'Sass'; }
pre.src-scheme:before { content: 'Scheme'; }
pre.src-screen:before { content: 'Gnu Screen'; }
pre.src-sed:before { content: 'Sed'; }
pre.src-sh:before { content: 'shell'; }
pre.src-sql:before { content: 'SQL'; }
pre.src-sqlite:before { content: 'SQLite'; }
/* additional languages in org.el's org-babel-load-languages alist */
pre.src-forth:before { content: 'Forth'; }
pre.src-io:before { content: 'IO'; }
pre.src-J:before { content: 'J'; }
pre.src-makefile:before { content: 'Makefile'; }
pre.src-maxima:before { content: 'Maxima'; }
pre.src-perl:before { content: 'Perl'; }
pre.src-picolisp:before { content: 'Pico Lisp'; }
pre.src-scala:before { content: 'Scala'; }
pre.src-shell:before { content: 'Shell Script'; }
pre.src-ebnf2ps:before { content: 'ebfn2ps'; }
/* additional language identifiers per "defun org-babel-execute"
in ob-*.el */
pre.src-cpp:before { content: 'C++'; }
pre.src-abc:before { content: 'ABC'; }
pre.src-coq:before { content: 'Coq'; }
pre.src-groovy:before { content: 'Groovy'; }
/* additional language identifiers from org-babel-shell-names in
ob-shell.el: ob-shell is the only babel language using a lambda to put
the execution function name together. */
pre.src-bash:before { content: 'bash'; }
pre.src-csh:before { content: 'csh'; }
pre.src-ash:before { content: 'ash'; }
pre.src-dash:before { content: 'dash'; }
pre.src-ksh:before { content: 'ksh'; }
pre.src-mksh:before { content: 'mksh'; }
pre.src-posh:before { content: 'posh'; }
/* Additional Emacs modes also supported by the LaTeX listings package */
pre.src-ada:before { content: 'Ada'; }
pre.src-asm:before { content: 'Assembler'; }
pre.src-caml:before { content: 'Caml'; }
pre.src-delphi:before { content: 'Delphi'; }
pre.src-html:before { content: 'HTML'; }
pre.src-idl:before { content: 'IDL'; }
pre.src-mercury:before { content: 'Mercury'; }
pre.src-metapost:before { content: 'MetaPost'; }
pre.src-modula-2:before { content: 'Modula-2'; }
pre.src-pascal:before { content: 'Pascal'; }
pre.src-ps:before { content: 'PostScript'; }
pre.src-prolog:before { content: 'Prolog'; }
pre.src-simula:before { content: 'Simula'; }
pre.src-tcl:before { content: 'tcl'; }
pre.src-tex:before { content: 'TeX'; }
pre.src-plain-tex:before { content: 'Plain TeX'; }
pre.src-verilog:before { content: 'Verilog'; }
pre.src-vhdl:before { content: 'VHDL'; }
pre.src-xml:before { content: 'XML'; }
pre.src-nxml:before { content: 'XML'; }
/* add a generic configuration mode; LaTeX export needs an additional
(add-to-list 'org-latex-listings-langs '(conf " ")) in .emacs */
pre.src-conf:before { content: 'Configuration File'; }
table { border-collapse:collapse; }
caption.t-above { caption-side: top; }
caption.t-bottom { caption-side: bottom; }
td, th { vertical-align:top; }
th.org-right { text-align: center; }
th.org-left { text-align: center; }
th.org-center { text-align: center; }
td.org-right { text-align: right; }
td.org-left { text-align: left; }
td.org-center { text-align: center; }
dt { font-weight: bold; }
.footpara { display: inline; }
.footdef { margin-bottom: 1em; }
.figure { padding: 1em; }
.figure p { text-align: center; }
.inlinetask {
padding: 10px;
border: 2px solid gray;
margin: 10px;
background: #ffffcc;
}
#org-div-home-and-up
{ text-align: right; font-size: 70%; white-space: nowrap; }
textarea { overflow-x: auto; }
.linenr { font-size: smaller }
.code-highlighted { background-color: #ffff00; }
.org-info-js_info-navigation { border-style: none; }
#org-info-js_console-label
{ font-size: 10px; font-weight: bold; white-space: nowrap; }
.org-info-js_search-highlight
{ background-color: #ffff00; color: #000000; font-weight: bold; }
.org-svg { width: 90%; }
/*]]>*/-->
< / style >
< link rel = "stylesheet" title = "Standard" href = "qmckl.css" type = "text/css" / >
< script type = "text/javascript" src = "org-info.js" >
/**
*
* @source: org-info.js
*
* @licstart The following is the entire license notice for the
* JavaScript code in org-info.js.
*
* Copyright (C) 2012-2019 Free Software Foundation, Inc.
*
*
* The JavaScript code in this tag is free software: you can
* redistribute it and/or modify it under the terms of the GNU
* General Public License (GNU GPL) as published by the Free Software
* Foundation, either version 3 of the License, or (at your option)
* any later version. The code is distributed WITHOUT ANY WARRANTY;
* without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU GPL for more details.
*
* As additional permission under GNU GPL version 3 section 7, you
* may distribute non-source (e.g., minimized or compacted) forms of
* that code without the copy of the GNU GPL normally required by
* section 4, provided you include this license notice and a URL
* through which recipients can access the Corresponding Source.
*
* @licend The above is the entire license notice
* for the JavaScript code in org-info.js.
*
*/
< / script >
< script type = "text/javascript" >
/*
@licstart The following is the entire license notice for the
JavaScript code in this tag.
Copyright (C) 2012-2019 Free Software Foundation, Inc.
The JavaScript code in this tag is free software: you can
redistribute it and/or modify it under the terms of the GNU
General Public License (GNU GPL) as published by the Free Software
Foundation, either version 3 of the License, or (at your option)
any later version. The code is distributed WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU GPL for more details.
As additional permission under GNU GPL version 3 section 7, you
may distribute non-source (e.g., minimized or compacted) forms of
that code without the copy of the GNU GPL normally required by
section 4, provided you include this license notice and a URL
through which recipients can access the Corresponding Source.
@licend The above is the entire license notice
for the JavaScript code in this tag.
*/
<!-- /* --> <![CDATA[/*> <!-- */
org_html_manager.set("TOC_DEPTH", "4");
org_html_manager.set("LINK_HOME", "index.html");
org_html_manager.set("LINK_UP", "");
org_html_manager.set("LOCAL_TOC", "1");
org_html_manager.set("VIEW_BUTTONS", "0");
org_html_manager.set("MOUSE_HINT", "underline");
org_html_manager.set("FIXED_TOC", "0");
org_html_manager.set("TOC", "1");
org_html_manager.set("VIEW", "info");
org_html_manager.setup(); // activate after the parameters are set
/*]]>*///-->
< / script >
< script type = "text/javascript" >
/*
@licstart The following is the entire license notice for the
JavaScript code in this tag.
Copyright (C) 2012-2019 Free Software Foundation, Inc.
The JavaScript code in this tag is free software: you can
redistribute it and/or modify it under the terms of the GNU
General Public License (GNU GPL) as published by the Free Software
Foundation, either version 3 of the License, or (at your option)
any later version. The code is distributed WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU GPL for more details.
As additional permission under GNU GPL version 3 section 7, you
may distribute non-source (e.g., minimized or compacted) forms of
that code without the copy of the GNU GPL normally required by
section 4, provided you include this license notice and a URL
through which recipients can access the Corresponding Source.
@licend The above is the entire license notice
for the JavaScript code in this tag.
*/
<!-- /* --> <![CDATA[/*> <!-- */
function CodeHighlightOn(elem, id)
{
var target = document.getElementById(id);
if(null != target) {
elem.cacheClassElem = elem.className;
elem.cacheClassTarget = target.className;
target.className = "code-highlighted";
elem.className = "code-highlighted";
}
}
function CodeHighlightOff(elem, id)
{
var target = document.getElementById(id);
if(elem.cacheClassElem)
elem.className = elem.cacheClassElem;
if(elem.cacheClassTarget)
target.className = elem.cacheClassTarget;
}
/*]]>*///-->
< / script >
< script type = "text/x-mathjax-config" >
MathJax.Hub.Config({
displayAlign: "center",
displayIndent: "0em",
"HTML-CSS": { scale: 100,
linebreaks: { automatic: "false" },
webFont: "TeX"
},
SVG: {scale: 100,
linebreaks: { automatic: "false" },
font: "TeX"},
NativeMML: {scale: 100},
TeX: { equationNumbers: {autoNumber: "AMS"},
MultLineWidth: "85%",
TagSide: "right",
TagIndent: ".8em"
}
});
< / script >
< script type = "text/javascript"
src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-AMS_HTML">< / script >
< / head >
< body >
< div id = "org-div-home-and-up" >
< a accesskey = "h" href = "" > UP < / a >
|
< a accesskey = "H" href = "index.html" > HOME < / a >
< / div > < div id = "content" >
< h1 class = "title" > Sherman-Morrison-Woodbury< / h1 >
< div id = "table-of-contents" >
< h2 > Table of Contents< / h2 >
< div id = "text-table-of-contents" >
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org259fd5a" > 1. Headers< / a > < / li >
< li > < a href = "#org39c4233" > 2. Naïve Sherman-Morrison< / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#orga50c3cf" > 2.1. < code > qmckl_sherman_morrison< / code > < / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#orgd058335" > 2.1.1. Requirements< / a > < / li >
< li > < a href = "#org07bf513" > 2.1.2. C header< / a > < / li >
< li > < a href = "#orge3c8214" > 2.1.3. C source< / a > < / li >
< li > < a href = "#orgacf5066" > 2.1.4. Performance< / a > < / li >
2021-09-10 17:38:14 +02:00
< / ul >
< / li >
< / ul >
< / li >
2022-03-07 23:23:02 +01:00
< li > < a href = "#orgaba04f3" > 3. Woodbury 2x2< / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org250b6da" > 3.1. < code > qmckl_woodbury_2< / code > < / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org4ed5d98" > 3.1.1. Requirements< / a > < / li >
< li > < a href = "#org71e06c1" > 3.1.2. C header< / a > < / li >
< li > < a href = "#org7f06128" > 3.1.3. C source< / a > < / li >
< li > < a href = "#org6940539" > 3.1.4. Performance< / a > < / li >
2021-09-10 17:38:14 +02:00
< / ul >
< / li >
< / ul >
< / li >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org7f96a2b" > 4. Woodbury 3x3< / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org1d0575f" > 4.1. < code > qmckl_woodbury_3< / code > < / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org39555a3" > 4.1.1. Requirements< / a > < / li >
< li > < a href = "#org9c6f00a" > 4.1.2. C header< / a > < / li >
< li > < a href = "#org69ffc69" > 4.1.3. C source< / a > < / li >
< li > < a href = "#org63ca566" > 4.1.4. Performance… < / a > < / li >
2021-09-10 17:38:14 +02:00
< / ul >
< / li >
< / ul >
< / li >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org2f6a975" > 5. Sherman-Morrison with update splitting< / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org5b93e6a" > 5.1. < code > qmckl_sherman_morrison_splitting< / code > < / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#orgcf73557" > 5.1.1. Requirements< / a > < / li >
< li > < a href = "#org8df1b93" > 5.1.2. C header< / a > < / li >
< li > < a href = "#org29f060a" > 5.1.3. C source< / a > < / li >
< li > < a href = "#orgc04a259" > 5.1.4. Performance… < / a > < / li >
2021-09-10 17:38:14 +02:00
< / ul >
< / li >
< / ul >
< / li >
2022-03-07 23:23:02 +01:00
< li > < a href = "#orgbd555d7" > 6. Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting< / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#orgd8036c8" > 6.1. < code > qmckl_sherman_morrison_smw32s< / code > < / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org0bfcd81" > 6.1.1. Requirements< / a > < / li >
< li > < a href = "#org5db5db6" > 6.1.2. C header< / a > < / li >
< li > < a href = "#orgd559f69" > 6.1.3. C source< / a > < / li >
< li > < a href = "#org98eb5b1" > 6.1.4. Performance… < / a > < / li >
2021-09-10 17:38:14 +02:00
< / ul >
< / li >
< / ul >
< / li >
2022-03-07 23:23:02 +01:00
< li > < a href = "#orgb240864" > 7. Helper Functions< / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org5595a56" > 7.1. < code > qmckl_slagel_splitting< / code > < / a >
2021-09-10 17:38:14 +02:00
< ul >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org1d0fba5" > 7.1.1. Requirements< / a > < / li >
< li > < a href = "#org773b895" > 7.1.2. C header< / a > < / li >
< li > < a href = "#org54a7cc2" > 7.1.3. C source< / a > < / li >
< li > < a href = "#org6822c84" > 7.1.4. Performance< / a > < / li >
2021-09-10 17:38:14 +02:00
< / ul >
< / li >
< / ul >
< / li >
2022-03-07 23:23:02 +01:00
< li > < a href = "#org70b2565" > 8. End of files< / a > < / li >
2021-09-10 17:38:14 +02:00
< / ul >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org259fd5a" class = "outline-2" >
< h2 id = "org259fd5a" > < span class = "section-number-2" > 1< / span > Headers< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-1" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "qmckl.h"< / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "assert.h"< / span >
< span style = "color: #483d8b;" > #ifdef< / span > HAVE_CONFIG_H
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "config.h"< / span >
< span style = "color: #483d8b;" > #endif< / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < math.h> < / span >
< span style = "color: #228b22;" > int< / span > < span style = "color: #0000ff;" > main< / span > () {
< span style = "color: #228b22;" > qmckl_context< / span > < span style = "color: #a0522d;" > context< / span > ;
context = qmckl_context_create();
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #a0522d;" > rc< / span > ;
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org39c4233" class = "outline-2" >
< h2 id = "org39c4233" > < span class = "section-number-2" > 2< / span > Naïve Sherman-Morrison< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-2" >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orga50c3cf" class = "outline-3" >
< h3 id = "orga50c3cf" > < span class = "section-number-3" > 2.1< / span > < code > qmckl_sherman_morrison< / code > < / h3 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-3" id = "text-2-1" >
< p >
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.
< / p >
< p >
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}
\]
< / p >
< p >
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.
< / p >
< p >
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.
< / p >
< p >
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}.
< / p >
2021-10-13 12:18:34 +02:00
< p >
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.
< / p >
2022-03-07 23:23:02 +01:00
< table id = "org40a720d" border = "2" cellspacing = "0" cellpadding = "6" rules = "groups" frame = "hsides" >
2021-09-10 17:38:14 +02:00
< colgroup >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< / colgroup >
< tbody >
< tr >
< td class = "org-left" > qmckl< sub > context< / sub > < / td >
< td class = "org-left" > context< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Global state< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > LDS< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > in< / td >
< td class = "org-left" > Leading dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-11-17 14:14:09 +01:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Dim< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-09-10 17:38:14 +02:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > N< sub > updates< / sub > < / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Number of rank-1 updates to be applied to Slater< sub > inv< / sub > < / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > Updates[N< sub > updates< / sub > *Dim]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the updates< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Updates< sub > index< / sub > [N< sub > updates< / sub > ]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the rank-1 updates< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > breakdown< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Break-down parameter on which to fail or not< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > Slater< sub > inv< / sub > [LDS*Dim]< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the inverse of a Slater-matrix< / td >
< / tr >
2021-10-13 12:18:34 +02:00
< tr >
< td class = "org-left" > double*< / td >
< td class = "org-left" > determinant< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Determinant of the Slater-matrix< / td >
< / tr >
2021-09-10 17:38:14 +02:00
< / tbody >
< / table >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgd058335" class = "outline-4" >
< h4 id = "orgd058335" > < span class = "section-number-4" > 2.1.1< / span > Requirements< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-2-1-1" >
< ul class = "org-ul" >
< li > < code > context< / code > is not < code > QMCKL_NULL_CONTEXT< / code > < / li >
2021-11-17 14:14:09 +01:00
< li > < code > LDS > = 2< / code > < / li >
2021-09-10 17:38:14 +02:00
< li > < code > Dim > = 2< / code > < / li >
< li > < code > N_updates > = 1< / code > < / li >
< li > < code > Updates< / code > is allocated with \(N_updates \times Dim\) elements< / li >
< li > < code > Updates_index< / code > is allocated with \(N_updates\) elements< / li >
< li > < code > breakdown< / code > is a small number such that \(0 < breakdown < < 1 \ ) < / li >
< li > < code > Slater_inv< / code > is allocated with \(Dim \times Dim\) elements< / li >
< / ul >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org07bf513" class = "outline-4" >
< h4 id = "org07bf513" > < span class = "section-number-4" > 2.1.2< / span > C header< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-2-1-2" >
< div class = "org-src-container" >
2021-09-14 10:01:20 +02:00
< pre class = "src src-c" > < span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_sherman_morrison< / span > (
2022-01-06 02:30:32 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > context< / span > < span style = "color: #a0522d;" > qmckl_context< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > LDS< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Dim< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > N_updates< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates_index< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > breakdown< / span > < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > Slater_inv< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > determinant< / span > * < span style = "color: #228b22;" > double< / span > * );
2021-09-10 17:38:14 +02:00
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orge3c8214" class = "outline-4" >
< h4 id = "orge3c8214" > < span class = "section-number-4" > 2.1.3< / span > C source< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-2-1-3" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < stdbool.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < math.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "qmckl.h"< / span >
2021-10-10 16:26:18 +02:00
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_sherman_morrison< / span > (< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > qmckl_context< / span > < span style = "color: #a0522d;" > context< / span > ,
2021-11-17 14:14:09 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > LDS< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > Dim< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > N_updates< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Updates< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > * < span style = "color: #a0522d;" > Updates_index< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > breakdown< / span > ,
2021-10-13 12:18:34 +02:00
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Slater_inv< / span > ,
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > determinant< / span > ) {
2021-09-18 18:13:58 +02:00
2021-10-06 23:45:54 +02:00
< span style = "color: #a020f0;" > if< / span > (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
< span style = "color: #a020f0;" > return< / span > QMCKL_NULL_CONTEXT;
}
2021-09-10 17:38:14 +02:00
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > C< / span > [Dim];
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > D< / span > [Dim];
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > l< / span > = 0;
// < span style = "color: #b22222;" > For each update< / span >
< span style = "color: #a020f0;" > while< / span > (l < N_updates) {
// < span style = "color: #b22222;" > C = A^{-1} x U_l< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
C[i] = 0;
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
C[i] += Slater_inv[i * LDS + j] * Updates[l * Dim + j];
2021-09-10 17:38:14 +02:00
}
}
// < span style = "color: #b22222;" > Denominator< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > den< / span > = 1 + C[Updates_index[l] - 1];
2021-10-13 12:18:34 +02:00
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > if< / span > (fabs(den) < breakdown) {
< span style = "color: #a020f0;" > return< / span > QMCKL_FAILURE;
}
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > iden< / span > = 1 / den;
2021-10-13 12:18:34 +02:00
// < span style = "color: #b22222;" > Update det(A)< / span >
< span style = "color: #a020f0;" > if< / span > (determinant != < span style = "color: #008b8b;" > NULL< / span > )
*determinant *= den;
2021-09-10 17:38:14 +02:00
// < span style = "color: #b22222;" > D = v^T x A^{-1}< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
D[j] = Slater_inv[(Updates_index[l] - 1) * LDS + j];
2021-09-10 17:38:14 +02:00
}
// < span style = "color: #b22222;" > A^{-1} = A^{-1} - C x D / den< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > update< / span > = C[i] * D[j] * iden;
2021-11-17 14:14:09 +01:00
Slater_inv[i * LDS + j] -= update;
2021-09-10 17:38:14 +02:00
}
}
l += 1;
}
< span style = "color: #a020f0;" > return< / span > QMCKL_SUCCESS;
}
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgacf5066" class = "outline-4" >
< h4 id = "orgacf5066" > < span class = "section-number-4" > 2.1.4< / span > Performance< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-2-1-4" >
< p >
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.
< / p >
< / div >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgaba04f3" class = "outline-2" >
< h2 id = "orgaba04f3" > < span class = "section-number-2" > 3< / span > Woodbury 2x2< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-3" >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org250b6da" class = "outline-3" >
< h3 id = "org250b6da" > < span class = "section-number-3" > 3.1< / span > < code > qmckl_woodbury_2< / code > < / h3 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-3" id = "text-3-1" >
< p >
The Woodbury 2x2 kernel. It is used to apply two rank-1 updates at once. The formula used in
this algorithm is called the Woodbury Matrix Identity
\[
(S + U V)^{-1} = S^{-1} - C B^{-1} D
\]
where
\(S\) is the Slater-matrix
\(U\) and \(V\) are the matrices containing the updates and the canonical basis matrix
\(S^{-1}\) is the inverse of the Slater-matrix
\(C:= S^{-1}U\), a Dim \(\times 2\) matrix
\(B := 1 + VC\), the \(2 \times 2\) matrix that is going to be inverted
\(D := VS^{-1}\), a \(2 \times Dim\) matrix
< / p >
2021-10-13 12:18:34 +02:00
< p >
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.
< / p >
2021-09-10 17:38:14 +02:00
2021-10-13 12:18:34 +02:00
2022-03-07 23:23:02 +01:00
< table id = "org86fb162" border = "2" cellspacing = "0" cellpadding = "6" rules = "groups" frame = "hsides" >
2021-09-10 17:38:14 +02:00
< colgroup >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< / colgroup >
< tbody >
< tr >
< td class = "org-left" > qmckl< sub > context< / sub > < / td >
< td class = "org-left" > context< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Global state< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > LDS< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > in< / td >
< td class = "org-left" > Leading dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-11-17 14:14:09 +01:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Dim< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-09-10 17:38:14 +02:00
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > Updates[2*Dim]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the updates< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Updates< sub > index< / sub > [2]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the rank-1 updates< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > breakdown< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Break-down parameter on which to fail or not< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > Slater< sub > inv< / sub > [LDS*Dim]< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the inverse of a Slater-matrix< / td >
< / tr >
2021-10-13 12:18:34 +02:00
< tr >
< td class = "org-left" > double*< / td >
< td class = "org-left" > determinant< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Determinant of Slater-matrix< / td >
< / tr >
2021-09-10 17:38:14 +02:00
< / tbody >
< / table >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org4ed5d98" class = "outline-4" >
< h4 id = "org4ed5d98" > < span class = "section-number-4" > 3.1.1< / span > Requirements< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-3-1-1" >
< ul class = "org-ul" >
< li > < code > context< / code > is not < code > qmckl_null_context< / code > < / li >
2021-11-17 14:14:09 +01:00
< li > < code > LDS > = 2< / code > < / li >
2021-09-10 17:38:14 +02:00
< li > < code > Dim > = 2< / code > < / li >
< li > < code > Updates< / code > is allocated with \(2 \times Dim\) elements< / li >
< li > < code > Updates_index< / code > is allocated with \(2\) elements< / li >
< li > < code > breakdown< / code > is a small number such that \(0 < breakdown < < 1 \ ) < / li >
< li > < code > Slater_inv< / code > is allocated with \(Dim \times Dim\) elements< / li >
< / ul >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org71e06c1" class = "outline-4" >
< h4 id = "org71e06c1" > < span class = "section-number-4" > 3.1.2< / span > C header< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-3-1-2" >
< div class = "org-src-container" >
2021-09-14 10:01:20 +02:00
< pre class = "src src-c" > < span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_woodbury_2< / span > (
2022-01-06 02:30:32 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > context< / span > < span style = "color: #a0522d;" > qmckl_context< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > LDS< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Dim< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates_index< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > breakdown< / span > < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > Slater_inv< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > determinant< / span > * < span style = "color: #228b22;" > double< / span > * );
2021-09-10 17:38:14 +02:00
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org7f06128" class = "outline-4" >
< h4 id = "org7f06128" > < span class = "section-number-4" > 3.1.3< / span > C source< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-3-1-3" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < stdbool.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < math.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "qmckl.h"< / span >
2021-10-10 16:26:18 +02:00
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_woodbury_2< / span > (< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > qmckl_context< / span > < span style = "color: #a0522d;" > context< / span > ,
2021-11-17 14:14:09 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > LDS< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > Dim< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Updates< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > * < span style = "color: #a0522d;" > Updates_index< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > breakdown< / span > ,
2021-10-13 12:18:34 +02:00
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Slater_inv< / span > ,
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > determinant< / span > ) {
2021-09-10 17:38:14 +02:00
/*
< span style = "color: #b22222;" > C := S^{-1} * U, dim x 2< / span >
< span style = "color: #b22222;" > B := 1 + V * C, 2 x 2< / span >
< span style = "color: #b22222;" > D := V * S^{-1}, 2 x dim< / span >
*/
2021-09-18 18:13:58 +02:00
2021-10-06 23:45:54 +02:00
< span style = "color: #a020f0;" > if< / span > (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
< span style = "color: #a020f0;" > return< / span > QMCKL_NULL_CONTEXT;
}
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > row1< / span > = (Updates_index[0] - 1);
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > row2< / span > = (Updates_index[1] - 1);
// < span style = "color: #b22222;" > Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE< / span >
// < span style = "color: #b22222;" > OF LAYOUT OF 'Updates' !!< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > C< / span > [2 * Dim];
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < 2; j++) {
C[i * 2 + j] = 0;
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > k< / span > = 0; k < Dim; k++) {
2021-11-17 14:14:09 +01:00
C[i * 2 + j] += Slater_inv[i * LDS + k] * Updates[Dim * j + k];
2021-09-10 17:38:14 +02:00
}
}
}
// < span style = "color: #b22222;" > Compute B = 1 + V * C< / span >
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B0< / span > = C[row1 * 2] + 1;
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B1< / span > = C[row1 * 2 + 1];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B2< / span > = C[row2 * 2];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B3< / span > = C[row2 * 2 + 1] + 1;
// < span style = "color: #b22222;" > Check if determinant of inverted matrix is not zero< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > det< / span > = B0 * B3 - B1 * B2;
< span style = "color: #a020f0;" > if< / span > (fabs(det) < breakdown) {
< span style = "color: #a020f0;" > return< / span > QMCKL_FAILURE;
}
2021-10-13 12:18:34 +02:00
// < span style = "color: #b22222;" > Update det(S) when passed< / span >
< span style = "color: #a020f0;" > if< / span > (determinant != < span style = "color: #008b8b;" > NULL< / span > )
*determinant *= det;
2021-09-10 17:38:14 +02:00
// < span style = "color: #b22222;" > Compute B^{-1} with explicit formula for 2x2 inversion< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > Binv< / span > [4], < span style = "color: #a0522d;" > idet< / span > = 1.0 / det;
Binv[0] = idet * B3;
Binv[1] = -1.0 * idet * B1;
Binv[2] = -1.0 * idet * B2;
Binv[3] = idet * B0;
// < span style = "color: #b22222;" > Compute tmp = B^{-1} x (V.S^{-1})< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > tmp< / span > [2 * Dim];
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < 2; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1 * LDS + j];
tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2 * LDS + j];
2021-09-10 17:38:14 +02:00
}
}
// < span style = "color: #b22222;" > Compute (S + U V)^{-1} = S^{-1} - C x tmp< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
Slater_inv[i * LDS + j] -= C[i * 2] * tmp[j];
Slater_inv[i * LDS + j] -= C[i * 2 + 1] * tmp[Dim + j];
2021-09-10 17:38:14 +02:00
}
}
< span style = "color: #a020f0;" > return< / span > QMCKL_SUCCESS;
}
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org6940539" class = "outline-4" >
< h4 id = "org6940539" > < span class = "section-number-4" > 3.1.4< / span > Performance< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-3-1-4" >
< p >
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.
< / p >
< / div >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org7f96a2b" class = "outline-2" >
< h2 id = "org7f96a2b" > < span class = "section-number-2" > 4< / span > Woodbury 3x3< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-4" >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org1d0575f" class = "outline-3" >
< h3 id = "org1d0575f" > < span class = "section-number-3" > 4.1< / span > < code > qmckl_woodbury_3< / code > < / h3 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-3" id = "text-4-1" >
< p >
The 3x3 version of the Woodbury 2x2 kernel. It is used to apply three
rank-1 updates at once. The formula used in this kernel is the same as for Woodbury 2x2,
except for the sizes of the following matrices:
< / p >
< p >
\(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
< / p >
2021-10-13 12:18:34 +02:00
< p >
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.
< / p >
2022-03-07 23:23:02 +01:00
< table id = "org7a650e3" border = "2" cellspacing = "0" cellpadding = "6" rules = "groups" frame = "hsides" >
2021-09-10 17:38:14 +02:00
< colgroup >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< / colgroup >
< tbody >
< tr >
< td class = "org-left" > qmckl< sub > context< / sub > < / td >
< td class = "org-left" > context< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Global state< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > LDS< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > in< / td >
< td class = "org-left" > Leading dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-11-17 14:14:09 +01:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Dim< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-09-10 17:38:14 +02:00
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > Updates[3*Dim]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the updates< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Updates< sub > index< / sub > [3]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the rank-1 updates< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > breakdown< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Break-down parameter on which to fail or not< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > Slater< sub > inv< / sub > [LDS*Dim]< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the inverse of a Slater-matrix< / td >
< / tr >
2021-10-13 12:18:34 +02:00
< tr >
< td class = "org-left" > double*< / td >
< td class = "org-left" > determinant< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Determinant of Slater-matrix< / td >
< / tr >
2021-09-10 17:38:14 +02:00
< / tbody >
< / table >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org39555a3" class = "outline-4" >
< h4 id = "org39555a3" > < span class = "section-number-4" > 4.1.1< / span > Requirements< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-4-1-1" >
< ul class = "org-ul" >
< li > < code > context< / code > is not < code > qmckl_null_context< / code > < / li >
2021-11-17 14:14:09 +01:00
< li > < code > LDS > = 2< / code > < / li >
2021-09-10 17:38:14 +02:00
< li > < code > Dim > = 2< / code > < / li >
< li > < code > Updates< / code > is allocated with \(3 \times Dim\) elements< / li >
< li > < code > Updates_index< / code > is allocated with \(3\) elements< / li >
< li > < code > breakdown< / code > is a small number such that \(0 < breakdown < < 1 \ ) < / li >
< li > < code > Slater_inv< / code > is allocated with \(Dim \times Dim\) elements< / li >
< / ul >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org9c6f00a" class = "outline-4" >
< h4 id = "org9c6f00a" > < span class = "section-number-4" > 4.1.2< / span > C header< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-4-1-2" >
< div class = "org-src-container" >
2021-09-14 10:01:20 +02:00
< pre class = "src src-c" > < span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_woodbury_3< / span > (
2022-01-06 02:30:32 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > context< / span > < span style = "color: #a0522d;" > qmckl_context< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > LDS< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Dim< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates_index< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > breakdown< / span > < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > Slater_inv< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > determinant< / span > * < span style = "color: #228b22;" > double< / span > * );
2021-09-10 17:38:14 +02:00
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org69ffc69" class = "outline-4" >
< h4 id = "org69ffc69" > < span class = "section-number-4" > 4.1.3< / span > C source< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-4-1-3" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < stdbool.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < math.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "qmckl.h"< / span >
2021-10-10 16:26:18 +02:00
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_woodbury_3< / span > (< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > qmckl_context< / span > < span style = "color: #a0522d;" > context< / span > ,
2021-11-17 14:14:09 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > LDS< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > Dim< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Updates< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > * < span style = "color: #a0522d;" > Updates_index< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > breakdown< / span > ,
2021-10-13 12:18:34 +02:00
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Slater_inv< / span > ,
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > determinant< / span > ) {
2021-09-10 17:38:14 +02:00
/*
< span style = "color: #b22222;" > C := S^{-1} * U, dim x 3< / span >
< span style = "color: #b22222;" > B := 1 + V * C, 3 x 3< / span >
< span style = "color: #b22222;" > D := V * S^{-1}, 3 x dim< / span >
*/
2021-09-18 18:13:58 +02:00
2021-10-06 23:45:54 +02:00
< span style = "color: #a020f0;" > if< / span > (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
< span style = "color: #a020f0;" > return< / span > QMCKL_NULL_CONTEXT;
}
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > row1< / span > = (Updates_index[0] - 1);
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > row2< / span > = (Updates_index[1] - 1);
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > row3< / span > = (Updates_index[2] - 1);
// < span style = "color: #b22222;" > Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE< / span >
// < span style = "color: #b22222;" > OF LAYOUT OF 'Updates' !!< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > C< / span > [3 * Dim];
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < 3; j++) {
C[i * 3 + j] = 0;
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > k< / span > = 0; k < Dim; k++) {
2021-11-17 14:14:09 +01:00
C[i * 3 + j] += Slater_inv[i * LDS + k] * Updates[Dim * j + k];
2021-09-10 17:38:14 +02:00
}
}
}
// < span style = "color: #b22222;" > Compute B = 1 + V.C< / span >
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B0< / span > = C[row1 * 3] + 1;
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B1< / span > = C[row1 * 3 + 1];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B2< / span > = C[row1 * 3 + 2];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B3< / span > = C[row2 * 3];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B4< / span > = C[row2 * 3 + 1] + 1;
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B5< / span > = C[row2 * 3 + 2];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B6< / span > = C[row3 * 3];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B7< / span > = C[row3 * 3 + 1];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > B8< / span > = C[row3 * 3 + 2] + 1;
// < span style = "color: #b22222;" > Check if determinant of B is not too close to zero< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > det< / span > ;
det = B0 * (B4 * B8 - B5 * B7) - B1 * (B3 * B8 - B5 * B6) +
B2 * (B3 * B7 - B4 * B6);
< span style = "color: #a020f0;" > if< / span > (fabs(det) < breakdown) {
< span style = "color: #a020f0;" > return< / span > QMCKL_FAILURE;
}
2021-10-13 12:18:34 +02:00
// < span style = "color: #b22222;" > Update det(Slater) if passed< / span >
< span style = "color: #a020f0;" > if< / span > (determinant != < span style = "color: #008b8b;" > NULL< / span > )
*determinant *= det;
2021-09-10 17:38:14 +02:00
// < span style = "color: #b22222;" > Compute B^{-1} with explicit formula for 3x3 inversion< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > Binv< / span > [9], < span style = "color: #a0522d;" > idet< / span > = 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;
// < span style = "color: #b22222;" > Compute tmp = B^{-1} x (V.S^{-1})< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > tmp< / span > [3 * Dim];
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < 3; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
tmp[i * Dim + j] = Binv[i * 3] * Slater_inv[row1 * LDS + j];
tmp[i * Dim + j] += Binv[i * 3 + 1] * Slater_inv[row2 * LDS + j];
tmp[i * Dim + j] += Binv[i * 3 + 2] * Slater_inv[row3 * LDS + j];
2021-09-10 17:38:14 +02:00
}
}
// < span style = "color: #b22222;" > Compute (S + U V)^{-1} = S^{-1} - C x tmp< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
Slater_inv[i * LDS + j] -= C[i * 3] * tmp[j];
Slater_inv[i * LDS + j] -= C[i * 3 + 1] * tmp[Dim + j];
Slater_inv[i * LDS + j] -= C[i * 3 + 2] * tmp[2 * Dim + j];
2021-09-10 17:38:14 +02:00
}
}
< span style = "color: #a020f0;" > return< / span > QMCKL_SUCCESS;
}
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org63ca566" class = "outline-4" >
< h4 id = "org63ca566" > < span class = "section-number-4" > 4.1.4< / span > Performance… < / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-4-1-4" >
< p >
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.
< / p >
< / div >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org2f6a975" class = "outline-2" >
< h2 id = "org2f6a975" > < span class = "section-number-2" > 5< / span > Sherman-Morrison with update splitting< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-5" >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org5b93e6a" class = "outline-3" >
< h3 id = "org5b93e6a" > < span class = "section-number-3" > 5.1< / span > < code > qmckl_sherman_morrison_splitting< / code > < / h3 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-3" id = "text-5-1" >
< p >
This is a variation on the 'Naive' Sherman-Morrison kernel. Whenever the denominator \(1+v_j^T S^{-1} u_j\) in
the Sherman-Morrison formula is deemed to be too close to zero, the update \(u_j\) is split in half:
\(u_j \rightarrow \frac{1}{2} u_j\). One half is applied immediately – necessarily increasing the value of the
denominator because of the split– while the other halve is put in a queue that will be applied when all the
remaining updates have been treated.
< / p >
< p >
The kernel is executed recursively until the queue is eiter empty and all
updates are applied successfully, or the size of the queue equals the number of initial updates. In the last
case the Slater-matrix that would have resulted from applying the updates is singular and therefore the
kernel exits with an exit code.
< / p >
2021-10-13 12:18:34 +02:00
< p >
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.
< / p >
2022-03-07 23:23:02 +01:00
< table id = "orgc828cca" border = "2" cellspacing = "0" cellpadding = "6" rules = "groups" frame = "hsides" >
2021-09-10 17:38:14 +02:00
< colgroup >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< / colgroup >
< tbody >
< tr >
< td class = "org-left" > qmckl< sub > context< / sub > < / td >
< td class = "org-left" > context< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Global state< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > LDS< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > in< / td >
< td class = "org-left" > Leading dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-11-17 14:14:09 +01:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Dim< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-09-10 17:38:14 +02:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > N< sub > updates< / sub > < / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Number of rank-1 updates to be applied to Slater< sub > inv< / sub > < / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > Updates[N< sub > updates< / sub > *Dim]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the updates< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Updates< sub > index< / sub > [N< sub > updates< / sub > ]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the rank-1 updates< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > breakdown< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Break-down parameter on which to fail or not< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > Slater< sub > inv< / sub > [LDS*Dim]< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the inverse of a Slater-matrix< / td >
< / tr >
2021-10-13 12:18:34 +02:00
< tr >
< td class = "org-left" > double*< / td >
< td class = "org-left" > determinant< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Determinant of the Slater-matrix< / td >
< / tr >
2021-09-10 17:38:14 +02:00
< / tbody >
< / table >
2021-10-13 12:18:34 +02:00
< p >
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.
< / p >
2021-09-10 17:38:14 +02:00
< / div >
2021-10-13 12:18:34 +02:00
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgcf73557" class = "outline-4" >
< h4 id = "orgcf73557" > < span class = "section-number-4" > 5.1.1< / span > Requirements< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-5-1-1" >
< ul class = "org-ul" >
< li > < code > context< / code > is not < code > QMCKL_NULL_CONTEXT< / code > < / li >
2021-11-17 14:14:09 +01:00
< li > < code > LDS > = 2< / code > < / li >
2021-09-10 17:38:14 +02:00
< li > < code > Dim > = 2< / code > < / li >
< li > < code > N_updates > = 1< / code > < / li >
< li > < code > Updates< / code > is allocated with \(N_updates \times Dim\) elements< / li >
< li > < code > Updates_index< / code > is allocated with \(N_updates\) elements< / li >
< li > < code > breakdown< / code > is a small number such that \(0 < breakdown < < 1 \ ) < / li >
< li > < code > Slater_inv< / code > is allocated with \(Dim \times Dim\) elements< / li >
< / ul >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org8df1b93" class = "outline-4" >
< h4 id = "org8df1b93" > < span class = "section-number-4" > 5.1.2< / span > C header< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-5-1-2" >
< div class = "org-src-container" >
2021-09-14 10:01:20 +02:00
< pre class = "src src-c" > < span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_sherman_morrison_splitting< / span > (
2022-01-06 02:30:32 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > context< / span > < span style = "color: #a0522d;" > qmckl_context< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > LDS< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Dim< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > N_updates< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates_index< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > breakdown< / span > < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > Slater_inv< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > determinant< / span > * < span style = "color: #228b22;" > double< / span > * );
2021-09-10 17:38:14 +02:00
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org29f060a" class = "outline-4" >
< h4 id = "org29f060a" > < span class = "section-number-4" > 5.1.3< / span > C source< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-5-1-3" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < stdbool.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "qmckl.h"< / span >
2021-10-10 16:26:18 +02:00
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_sherman_morrison_splitting< / span > (< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > qmckl_context< / span > < span style = "color: #a0522d;" > context< / span > ,
2021-11-17 14:14:09 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > LDS< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > Dim< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > N_updates< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Updates< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > * < span style = "color: #a0522d;" > Updates_index< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > breakdown< / span > ,
2021-10-13 12:18:34 +02:00
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Slater_inv< / span > ,
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > determinant< / span > ) {
2021-09-10 17:38:14 +02:00
2021-10-06 23:45:54 +02:00
< span style = "color: #a020f0;" > if< / span > (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
< span style = "color: #a020f0;" > return< / span > QMCKL_NULL_CONTEXT;
}
2021-10-10 16:26:18 +02:00
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > later_updates< / span > [Dim * N_updates];
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > later_index< / span > [N_updates];
2021-09-10 17:38:14 +02:00
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > later< / span > = 0;
2021-11-17 14:14:09 +01:00
(< span style = "color: #228b22;" > void< / span > ) qmckl_slagel_splitting(LDS, Dim, N_updates, Updates, Updates_index,
2021-10-13 12:18:34 +02:00
breakdown, Slater_inv, later_updates, later_index, & later, determinant);
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > if< / span > (later > 0) {
2021-11-17 14:14:09 +01:00
(< span style = "color: #228b22;" > void< / span > ) qmckl_sherman_morrison_splitting(context, LDS, Dim, later,
2021-10-13 12:18:34 +02:00
later_updates, later_index, breakdown, Slater_inv, determinant);
2021-09-10 17:38:14 +02:00
}
< span style = "color: #a020f0;" > return< / span > QMCKL_SUCCESS;
}
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgc04a259" class = "outline-4" >
< h4 id = "orgc04a259" > < span class = "section-number-4" > 5.1.4< / span > Performance… < / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-5-1-4" >
< p >
This kernel performs best when there are 2 or more rank-1 update cycles and fail-rate is high.
< / p >
< / div >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgbd555d7" class = "outline-2" >
< h2 id = "orgbd555d7" > < span class = "section-number-2" > 6< / span > Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-6" >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgd8036c8" class = "outline-3" >
< h3 id = "orgd8036c8" > < span class = "section-number-3" > 6.1< / span > < code > qmckl_sherman_morrison_smw32s< / code > < / h3 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-3" id = "text-6-1" >
< p >
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.
< / p >
2021-10-13 12:18:34 +02:00
< p >
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.
< / p >
2022-03-07 23:23:02 +01:00
< table id = "orge0bbf91" border = "2" cellspacing = "0" cellpadding = "6" rules = "groups" frame = "hsides" >
2021-09-10 17:38:14 +02:00
< colgroup >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< / colgroup >
< tbody >
< tr >
< td class = "org-left" > qmckl< sub > context< / sub > < / td >
< td class = "org-left" > context< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Global state< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > LDS< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > in< / td >
< td class = "org-left" > Leading dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-11-17 14:14:09 +01:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Dim< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-09-10 17:38:14 +02:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > N< sub > updates< / sub > < / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Number of rank-1 updates to be applied to Slater< sub > inv< / sub > < / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > Updates[N< sub > updates< / sub > *Dim]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the updates< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Updates< sub > index< / sub > [N< sub > updates< / sub > ]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the rank-1 updates< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > breakdown< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Break-down parameter on which to fail or not< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > Slater< sub > inv< / sub > [LDS*Dim]< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the inverse of a Slater-matrix< / td >
< / tr >
2021-10-13 12:18:34 +02:00
< tr >
< td class = "org-left" > double*< / td >
< td class = "org-left" > determinant< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Determinant of the Slater-matrix< / td >
< / tr >
2021-09-10 17:38:14 +02:00
< / tbody >
< / table >
< / div >
2021-10-13 12:18:34 +02:00
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org0bfcd81" class = "outline-4" >
< h4 id = "org0bfcd81" > < span class = "section-number-4" > 6.1.1< / span > Requirements< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-6-1-1" >
< ul class = "org-ul" >
< li > < code > context< / code > is not < code > QMCKL_NULL_CONTEXT< / code > < / li >
2021-11-17 14:14:09 +01:00
< li > < code > LDS > = 2< / code > < / li >
2021-09-10 17:38:14 +02:00
< li > < code > Dim > = 2< / code > < / li >
< li > < code > N_updates > = 1< / code > < / li >
< li > < code > Updates< / code > is allocated with \(N_updates \times Dim\) elements< / li >
< li > < code > Updates_index< / code > is allocated with \(N_updates\) elements< / li >
< li > < code > breakdown< / code > is a small number such that \(0 < breakdown < < 1 \ ) < / li >
< li > < code > Slater_inv< / code > is allocated with \(Dim \times Dim\) elements< / li >
< / ul >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org5db5db6" class = "outline-4" >
< h4 id = "org5db5db6" > < span class = "section-number-4" > 6.1.2< / span > C header< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-6-1-2" >
< div class = "org-src-container" >
2021-09-14 10:01:20 +02:00
< pre class = "src src-c" > < span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_sherman_morrison_smw32s< / span > (
2022-01-06 02:30:32 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > context< / span > < span style = "color: #a0522d;" > qmckl_context< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > LDS< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Dim< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > N_updates< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates_index< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > breakdown< / span > < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > Slater_inv< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > determinant< / span > * < span style = "color: #228b22;" > double< / span > * );
2021-09-10 17:38:14 +02:00
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgd559f69" class = "outline-4" >
< h4 id = "orgd559f69" > < span class = "section-number-4" > 6.1.3< / span > C source< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-6-1-3" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < stdbool.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "qmckl.h"< / span >
2021-10-10 16:26:18 +02:00
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_sherman_morrison_smw32s< / span > (< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > qmckl_context< / span > < span style = "color: #a0522d;" > context< / span > ,
2021-11-17 14:14:09 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > LDS< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > Dim< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > N_updates< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Updates< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > * < span style = "color: #a0522d;" > Updates_index< / span > ,
2021-10-10 16:26:18 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > breakdown< / span > ,
2021-10-13 12:18:34 +02:00
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > Slater_inv< / span > ,
< span style = "color: #228b22;" > double< / span > * < span style = "color: #a0522d;" > determinant< / span > ) {
2021-09-18 18:13:58 +02:00
2021-10-06 23:45:54 +02:00
< span style = "color: #a020f0;" > if< / span > (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
< span style = "color: #a020f0;" > return< / span > QMCKL_NULL_CONTEXT;
}
2021-09-10 17:38:14 +02:00
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #a0522d;" > rc< / span > ;
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > n_of_3blocks< / span > = N_updates / 3;
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > remainder< / span > = N_updates % 3;
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > length_3block< / span > = 3 * Dim;
// < span style = "color: #b22222;" > Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with< / span >
// < span style = "color: #b22222;" > Woodbury 3x3 kernel< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > later_updates< / span > [Dim * N_updates];
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > later_index< / span > [N_updates];
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > later< / span > = 0;
< span style = "color: #a020f0;" > if< / span > (n_of_3blocks > 0) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < n_of_3blocks; i++) {
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > *< span style = "color: #a0522d;" > Updates_3block< / span > = & Updates[i * length_3block];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > *< span style = "color: #a0522d;" > Updates_index_3block< / span > = & Updates_index[i * 3];
2021-11-17 14:14:09 +01:00
rc = qmckl_woodbury_3(context, LDS, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, determinant);
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > if< / span > (rc != 0) { // < span style = "color: #b22222;" > Send the entire block to slagel_splitting< / span >
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > l< / span > = 0;
2021-11-17 14:14:09 +01:00
rc = qmckl_slagel_splitting(LDS, Dim, 3, Updates_3block, Updates_index_3block,
2021-10-13 12:18:34 +02:00
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, & l, determinant);
2021-09-10 17:38:14 +02:00
later = later + l;
}
}
}
2021-09-18 18:13:58 +02:00
// < span style = "color: #b22222;" > Apply last remaining block of 2 updates with Woodbury 2x2 kernel< / span >
< span style = "color: #a020f0;" > if< / span > (remainder == 2) {
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > *< span style = "color: #a0522d;" > Updates_2block< / span > = & Updates[n_of_3blocks * length_3block];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > *< span style = "color: #a0522d;" > Updates_index_2block< / span > = & Updates_index[3 * n_of_3blocks];
2021-11-17 14:14:09 +01:00
rc = qmckl_woodbury_2(context, LDS, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv, determinant);
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > if< / span > (rc != 0) { // < span style = "color: #b22222;" > Send the entire block to slagel_splitting< / span >
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > l< / span > = 0;
2021-11-17 14:14:09 +01:00
(< span style = "color: #228b22;" > void< / span > ) qmckl_slagel_splitting(LDS, Dim, 2, Updates_2block, Updates_index_2block,
2021-10-13 12:18:34 +02:00
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, & l, determinant);
2021-09-10 17:38:14 +02:00
later = later + l;
}
}
2021-09-18 18:13:58 +02:00
// < span style = "color: #b22222;" > Apply last remaining update with slagel_splitting< / span >
< span style = "color: #a020f0;" > else< / span > < span style = "color: #a020f0;" > if< / span > (remainder == 1) {
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > *< span style = "color: #a0522d;" > Updates_1block< / span > = & Updates[n_of_3blocks * length_3block];
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > *< span style = "color: #a0522d;" > Updates_index_1block< / span > = & Updates_index[3 * n_of_3blocks];
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > l< / span > = 0;
2021-11-17 14:14:09 +01:00
(< span style = "color: #228b22;" > void< / span > ) qmckl_slagel_splitting(LDS, Dim, 1, Updates_1block, Updates_index_1block,
2021-10-13 12:18:34 +02:00
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, & l, determinant);
2021-09-10 17:38:14 +02:00
later = later + l;
}
< span style = "color: #a020f0;" > if< / span > (later > 0) {
2021-11-17 14:14:09 +01:00
(< span style = "color: #228b22;" > void< / span > ) qmckl_sherman_morrison_splitting(context, LDS, Dim, later, later_updates, later_index, breakdown, Slater_inv, determinant);
2021-09-10 17:38:14 +02:00
}
< span style = "color: #a020f0;" > return< / span > QMCKL_SUCCESS;
}
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org98eb5b1" class = "outline-4" >
< h4 id = "org98eb5b1" > < span class = "section-number-4" > 6.1.4< / span > Performance… < / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-6-1-4" >
< p >
This kernel performs best for update cycles with 2 or more rank-1 updates and the fail-rate is low.
< / p >
< / div >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-orgb240864" class = "outline-2" >
< h2 id = "orgb240864" > < span class = "section-number-2" > 7< / span > Helper Functions< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-7" >
< p >
2021-10-06 23:45:54 +02:00
Private helper-functions that are used by the Sherman-Morrison-Woodbury kernels.
2021-09-10 17:38:14 +02:00
These functions can only be used internally by the kernels in this module.
< / p >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org5595a56" class = "outline-3" >
< h3 id = "org5595a56" > < span class = "section-number-3" > 7.1< / span > < code > qmckl_slagel_splitting< / code > < / h3 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-3" id = "text-7-1" >
< p >
< code > qmckl_slagel_splitting< / code > 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.
< / p >
< p >
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}.
< / p >
2021-10-13 12:18:34 +02:00
< p >
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.
< / p >
2022-03-07 23:23:02 +01:00
< table id = "org41cceb3" border = "2" cellspacing = "0" cellpadding = "6" rules = "groups" frame = "hsides" >
2021-09-10 17:38:14 +02:00
< colgroup >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< col class = "org-left" / >
< / colgroup >
< tbody >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > LDS< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > in< / td >
< td class = "org-left" > Leading dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-11-17 14:14:09 +01:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Dim< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Dimension of Slater< sub > inv< / sub > < / td >
< / tr >
2021-09-10 17:38:14 +02:00
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > N< sub > updates< / sub > < / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Number of rank-1 updates to be applied to Slater< sub > inv< / sub > < / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > Updates[N< sub > updates< / sub > *Dim]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing the rank-1 updates< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > Updates< sub > index< / sub > [N< sub > updates< / sub > ]< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Array containing positions of the rank-1 updates< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > breakdown< / td >
< td class = "org-left" > in< / td >
< td class = "org-left" > Break-down parameter on which to fail or not< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
2021-11-17 14:14:09 +01:00
< td class = "org-left" > Slater< sub > inv< / sub > [LDS*Dim]< / td >
2021-09-10 17:38:14 +02:00
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the inverse Slater-matrix< / td >
< / tr >
< tr >
< td class = "org-left" > double< / td >
< td class = "org-left" > later< sub > updates< / sub > [Dim * N< sub > updates< / sub > ]< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the split updates for later< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > later< sub > index< / sub > [N< sub > updates< / sub > ]< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Array containing the positions of the split updates for later< / td >
< / tr >
< tr >
< td class = "org-left" > uint64< sub > t< / sub > < / td >
< td class = "org-left" > later< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Number of split updates for later< / td >
< / tr >
2021-10-13 12:18:34 +02:00
< tr >
< td class = "org-left" > double*< / td >
< td class = "org-left" > determinant< / td >
< td class = "org-left" > inout< / td >
< td class = "org-left" > Determinant of the Slater-matrix< / td >
< / tr >
2021-09-10 17:38:14 +02:00
< / tbody >
< / table >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org1d0fba5" class = "outline-4" >
< h4 id = "org1d0fba5" > < span class = "section-number-4" > 7.1.1< / span > Requirements< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-7-1-1" >
< ul class = "org-ul" >
2021-11-17 14:14:09 +01:00
< li > < code > LDS > = 2< / code > < / li >
2021-09-10 17:38:14 +02:00
< li > < code > Dim > = 2< / code > < / li >
< li > < code > N_updates > = 1< / code > < / li >
< li > < code > Updates< / code > is allocated with \(N_updates \times Dim\) elements< / li >
< li > < code > Updates_index< / code > is allocated with \(N_updates\) elements< / li >
< li > < code > breakdown< / code > is a small number such that \(0 < breakdown < < 1 \ ) < / li >
< li > < code > Slater_inv< / code > is allocated with \(Dim \times Dim\) elements< / li >
< li > < code > later_updates< / code > is allocated with \(later \times Dim\) elements< / li >
< li > < code > later_index< / code > is allocated with \(N_updates\) elements< / li >
< li > < code > later > = 0< / code > < / li >
< / ul >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org773b895" class = "outline-4" >
< h4 id = "org773b895" > < span class = "section-number-4" > 7.1.2< / span > C header< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-7-1-2" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #228b22;" > double< / span > < span style = "color: #0000ff;" > qmckl_slagel_splitting< / span > (
2022-01-06 02:30:32 +01:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > LDS< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Dim< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > N_updates< / span > < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > Updates_index< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > breakdown< / span > < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > Slater_inv< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > later_updates< / span > * < span style = "color: #228b22;" > double< / span > ,
< span style = "color: #228b22;" > later_index< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #228b22;" > later< / span > * < span style = "color: #a0522d;" > uint64_t< / span > ,
< span style = "color: #228b22;" > determinant< / span > * < span style = "color: #228b22;" > double< / span > * );
2021-09-10 17:38:14 +02:00
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org54a7cc2" class = "outline-4" >
< h4 id = "org54a7cc2" > < span class = "section-number-4" > 7.1.3< / span > C source< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-7-1-3" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < stdbool.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > < math.h> < / span >
< span style = "color: #483d8b;" > #include< / span > < span style = "color: #8b2252;" > "qmckl.h"< / span >
2021-11-17 14:14:09 +01:00
< span style = "color: #228b22;" > qmckl_exit_code< / span > < span style = "color: #0000ff;" > qmckl_slagel_splitting< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > LDS< / span > ,
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > Dim< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > N_updates< / span > ,
2021-10-06 23:45:54 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > *< span style = "color: #a0522d;" > Updates< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > uint64_t< / span > *< span style = "color: #a0522d;" > Updates_index< / span > ,
< span style = "color: #a020f0;" > const< / span > < span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > breakdown< / span > ,
< span style = "color: #228b22;" > double< / span > *< span style = "color: #a0522d;" > Slater_inv< / span > ,
2021-10-06 23:45:54 +02:00
< span style = "color: #228b22;" > double< / span > *< span style = "color: #a0522d;" > later_updates< / span > ,
2021-09-10 17:38:14 +02:00
< span style = "color: #228b22;" > uint64_t< / span > *< span style = "color: #a0522d;" > later_index< / span > ,
2021-10-13 12:18:34 +02:00
< span style = "color: #228b22;" > uint64_t< / span > *< span style = "color: #a0522d;" > later< / span > ,
< span style = "color: #228b22;" > double< / span > *< span style = "color: #a0522d;" > determinant< / span > ) {
2021-09-10 17:38:14 +02:00
// < span style = "color: #b22222;" > #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl.< / span >
// < span style = "color: #b22222;" > std::cerr < < "Called slagel_splitting with " < < N_updates < < " updates" < < std::endl;< / span >
// < span style = "color: #b22222;" > #endif< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > C< / span > [Dim];
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > D< / span > [Dim];
< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > l< / span > = 0;
// < span style = "color: #b22222;" > For each update< / span >
< span style = "color: #a020f0;" > while< / span > (l < N_updates) {
// < span style = "color: #b22222;" > C = S^{-1} x U_l< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
C[i] = 0;
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
C[i] += Slater_inv[i * LDS + j] * Updates[l * Dim + j];
2021-09-10 17:38:14 +02:00
}
}
// < span style = "color: #b22222;" > Denominator< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > den< / span > = 1 + C[Updates_index[l] - 1];
< span style = "color: #a020f0;" > if< / span > (fabs(den) < breakdown) { // < span style = "color: #b22222;" > Here is decided to split the update, or not.< / span >
// < span style = "color: #b22222;" > U_l = U_l / 2: split the update in 2 equal halves and save the second halve in later_updates< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 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];
} // < span style = "color: #b22222;" > From here onwards we continue with applying the first havel of the update to Slater_inv< / span >
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > iden< / span > = 1 / den;
2021-10-13 12:18:34 +02:00
< span style = "color: #a020f0;" > if< / span > (determinant != < span style = "color: #008b8b;" > NULL< / span > )
*determinant *= den;
2021-09-10 17:38:14 +02:00
// < span style = "color: #b22222;" > D = v^T x S^{-1}< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
2021-11-17 14:14:09 +01:00
D[j] = Slater_inv[(Updates_index[l] - 1) * LDS + j];
2021-09-10 17:38:14 +02:00
}
// < span style = "color: #b22222;" > S^{-1} = S^{-1} - C x D / den< / span >
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > i< / span > = 0; i < Dim; i++) {
< span style = "color: #a020f0;" > for< / span > (< span style = "color: #228b22;" > uint64_t< / span > < span style = "color: #a0522d;" > j< / span > = 0; j < Dim; j++) {
< span style = "color: #228b22;" > double< / span > < span style = "color: #a0522d;" > update< / span > = C[i] * D[j] * iden;
2021-11-17 14:14:09 +01:00
Slater_inv[i * LDS + j] -= update;
2021-09-10 17:38:14 +02:00
}
}
l += 1;
}
< span style = "color: #a020f0;" > return< / span > QMCKL_SUCCESS;
}
< / pre >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org6822c84" class = "outline-4" >
< h4 id = "org6822c84" > < span class = "section-number-4" > 7.1.4< / span > Performance< / h4 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-4" id = "text-7-1-4" >
< p >
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.
< / p >
< / div >
< / div >
< / div >
< / div >
2022-03-07 23:23:02 +01:00
< div id = "outline-container-org70b2565" class = "outline-2" >
< h2 id = "org70b2565" > < span class = "section-number-2" > 8< / span > End of files< / h2 >
2021-09-10 17:38:14 +02:00
< div class = "outline-text-2" id = "text-8" >
< div class = "org-src-container" >
< pre class = "src src-c" > < span style = "color: #0000ff;" > assert< / span > (< span style = "color: #228b22;" > qmckl_context_destroy< / span > (< span style = "color: #a0522d;" > context< / span > ) == QMCKL_SUCCESS);
< span style = "color: #a020f0;" > return< / span > 0;
2021-09-18 18:13:58 +02:00
}
2021-09-10 17:38:14 +02:00
< / pre >
< / div >
< / div >
< / div >
< / div >
< div id = "postamble" class = "status" >
< p class = "author" > Author: TREX CoE< / p >
2022-03-07 23:23:02 +01:00
< p class = "date" > Created: 2022-03-07 Mon 22:22< / p >
2021-09-10 17:38:14 +02:00
< p class = "validation" > < a href = "http://validator.w3.org/check?uri=referer" > Validate< / a > < / p >
< / div >
< / body >
< / html >