From 6e047046f4f15f4e66e285440abc5a7adad2c0ea Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Mon, 19 Jul 2021 18:25:10 +0200 Subject: [PATCH] Get the sherman morrison to compile. Tests still dont compile. #25 --- org/qmckl_sherman_morrison_woodbury.org | 126 ++++++++++++++++++++---- 1 file changed, 106 insertions(+), 20 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 17ca931..255c0f4 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -21,6 +21,8 @@ Low- and high-level functions that use the Sherman-Morrison and Woodbury matrix #define THRESHOLD 1e-3 #endif +#include "qmckl.h" + int main() { qmckl_context context; context = qmckl_context_create(); @@ -28,7 +30,7 @@ int main() { #+end_src * Naïve Sherman-Morrison - + ** ~qmckl_sherman_morrison~ :PROPERTIES: :Name: qmckl_sherman_morrison @@ -37,6 +39,7 @@ int main() { :END: The Sherman-Morrison formula + \begin{align} S_k^{-1} &= (S_l + U_k)^-1 \\ &= S_l^{-1} - \frac{S_l^{-1}U_kS_l}{1+\underline{v}_k^tS_l^{-1}\underline{u}_k} @@ -44,12 +47,12 @@ int main() { #+NAME: qmckl_sherman_morrison_args - | qmckl_context | context | in | Global state | - | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | - | uint | Dim | in | Leading dimension of Slater_inv | - | uint | N_updates | in | Number of rank-1 updates to be applied to Slater_inv | - | double | Updates[N_updates*Dim] | in | Array containing the updates | - | double | Updates_index | in | Array containing the rank-1 updates | + | qmckl_context | context | in | Global state | + | uint64_t | Dim | in | Leading dimension of Slater_inv | + | uint64_t | N_updates | in | Number of rank-1 updates to be applied to Slater_inv | + | double | Updates[N_updates*Dim] | in | Array containing the updates | + | uint64_t | Updates_index[N_updates] | in | Array containing the rank-1 updates | + | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements @@ -61,20 +64,59 @@ int main() { #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org - qmckl_exit_code qmckl_sherman_morrison ( - const qmckl_context context, - double *Slater_inv, - unsigned int Dim, - unsigned int N_updates, - double *Updates, - unsigned int *Updates_index ); +// Sherman-Morrison-Woodbury break-down threshold +#ifndef THRESHOLD +#define THRESHOLD 1e-3 +#endif +static double threshold(); + + qmckl_exit_code qmckl_sherman_morrison_c ( + const qmckl_context context, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + double* Slater_inv ); #+end_src -*** Source +*** Source Fortran - #+begin_src c :tangle (eval c_func) :comments org -bool qmckl_sherman_morrison(double *Slater_inv, unsigned int Dim, unsigned int N_updates, - double *Updates, unsigned int *Updates_index) { + #+begin_src f90 :tangle (eval f) +integer function qmckl_sherman_morrison_f(context, Slater_inv, Dim, N_updates, & + Updates, Updates_index) result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + integer*8 , intent(in) :: Dim, N_updates + integer*8 , intent(in) :: Updates_index(N_updates) + real*8 , intent(in) :: Updates(N_updates*Dim) + real*8 , intent(inout) :: Slater_inv(Dim*Dim) + !logical, external :: qmckl_sherman_morrison_f + info = qmckl_sherman_morrison(context, Dim, N_updates, Updates, Updates_index, Slater_inv) +end function qmckl_sherman_morrison_f + #+end_src + +*** Source C + + #+begin_src c :tangle (eval c) :comments org +#include +#include "qmckl.h" + +// Sherman-Morrison-Woodbury break-down threshold +static double threshold() { + const double threshold = THRESHOLD; +#ifdef DEBUG + std::cerr << "Break-down threshold set to: " << threshold << std::endl; +#endif + return threshold; +} + +qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + double * Slater_inv) { #ifdef DEBUG std::cerr << "Called qmckl_sherman_morrison with " << N_updates << " updates" << std::endl; #endif @@ -95,7 +137,7 @@ bool qmckl_sherman_morrison(double *Slater_inv, unsigned int Dim, unsigned int N // Denominator double den = 1 + C[Updates_index[l] - 1]; - if (std::fabs(den) < threshold()) { + if (fabs(den) < threshold()) { return false; } double iden = 1 / den; @@ -115,7 +157,7 @@ bool qmckl_sherman_morrison(double *Slater_inv, unsigned int Dim, unsigned int N l += 1; } - return true; + return QMCKL_SUCCESS; } #+end_src @@ -128,8 +170,52 @@ bool qmckl_sherman_morrison(double *Slater_inv, unsigned int Dim, unsigned int N #+CALL: generate_c_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_sherman_morrison & + (context, Dim, N_updates, Updates, Updates_index, Slater_inv) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*Dim) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) + + integer(c_int32_t), external :: qmckl_sherman_morrison_c + info = qmckl_sherman_morrison_c & + (context, Dim, N_updates, Updates, Updates_index, Slater_inv) + + end function qmckl_sherman_morrison + #+end_src + #+CALL: generate_f_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_sherman_morrison & + (context, Dim, N_updates, Updates, Updates_index, Slater_inv) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*Dim) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) + + end function qmckl_sherman_morrison + end interface + #+end_src + *** Test :noexport: [TODO: FMJC] Write tests for the Sherman-Morrison part.