Possibility to compute only values

org/qmckl_mo.org

@@ -92,10 +92,12 @@ int main() {
 
   Computed data:
 
-  |---------------+--------------------------+-------------------------------------------------------------------------------------|
+  |-----------------+--------------------------+-------------------------------------------------------------------------------------|
-  | ~mo_vgl~      | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions                           |
+  | ~mo_value~      | ~[point_num][mo_num]~    | Value of the MOs at point positions                                                 |
-  | ~mo_vgl_date~ | ~uint64_t~               | Late modification date of Value, gradients, Laplacian of the MOs at point positions |
+  | ~mo_value_date~ | ~uint64_t~               | Late modification date of the value of the MOs at point positions                   |
-  |---------------+--------------------------+-------------------------------------------------------------------------------------|
+  | ~mo_vgl~        | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions                           |
+  | ~mo_vgl_date~   | ~uint64_t~               | Late modification date of Value, gradients, Laplacian of the MOs at point positions |
+  |-----------------+--------------------------+-------------------------------------------------------------------------------------|
 
 ** Data structure
 
@@ -106,7 +108,9 @@ typedef struct qmckl_mo_basis_struct {
   double * restrict coefficient_t;
 
   double * restrict mo_vgl;
+  double * restrict mo_value;
   uint64_t  mo_vgl_date;
+  uint64_t  mo_value_date;
 
   int32_t   uninitialized;
   bool      provided;
@@ -418,7 +422,464 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) {
 
 * Computation
 
-** Computation of MOs
+** Computation of MOs: values only
+
+*** Get
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value(qmckl_context context,
+                            double* const mo_value,
+                            const int64_t size_max);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value(qmckl_context context,
+                            double* const mo_value,
+                            const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_exit_code rc;
+
+  rc = qmckl_provide_ao_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  rc = qmckl_provide_mo_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_mo_basis_mo_value",
+                           "input array too small");
+  }
+  memcpy(mo_value, ctx->mo_basis.mo_value, sze * sizeof(double));
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+  interface
+     integer(c_int32_t) function qmckl_get_mo_basis_mo_value (context, &
+          mo_value, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+
+       integer (c_int64_t) , intent(in)  , value :: context
+       double precision,     intent(out)         :: mo_value(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_mo_basis_mo_value
+  end interface
+    #+end_src
+
+    Uses the given array to compute the values.
+
+    #+begin_src c :comments org :tangle (eval h_func) :noweb yes
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
+                                     double* const mo_value,
+                                     const int64_t size_max);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code
+qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
+                                     double* const mo_value,
+                                     const int64_t size_max)
+{
+
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_CONTEXT,
+                           "qmckl_get_mo_basis_mo_value",
+                           NULL);
+  }
+
+  qmckl_exit_code rc;
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num;
+  if (size_max < sze) {
+    return qmckl_failwith( context,
+                           QMCKL_INVALID_ARG_3,
+                           "qmckl_get_mo_basis_mo_value",
+                           "input array too small");
+  }
+
+  rc = qmckl_context_touch(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  double* old_array = ctx->mo_basis.mo_value;
+
+  ctx->mo_basis.mo_value = mo_value;
+
+  rc = qmckl_provide_mo_value(context);
+  if (rc != QMCKL_SUCCESS) return rc;
+
+  ctx->mo_basis.mo_value = old_array;
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+    #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
+  interface
+     integer(c_int32_t) function qmckl_get_mo_basis_mo_value_inplace (context, &
+          mo_value, size_max) bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+       integer (c_int64_t) , intent(in)  , value :: context
+       double precision,     intent(out)         :: mo_value(*)
+       integer (c_int64_t) , intent(in)  , value :: size_max
+     end function qmckl_get_mo_basis_mo_value_inplace
+  end interface
+    #+end_src
+
+*** Provide
+
+    #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none
+qmckl_exit_code qmckl_provide_mo_value(qmckl_context context);
+    #+end_src
+
+    #+begin_src c :comments org :tangle (eval c) :noweb yes  :exports none
+qmckl_exit_code qmckl_provide_mo_value(qmckl_context context)
+{
+
+  qmckl_exit_code rc;
+  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
+    return QMCKL_NULL_CONTEXT;
+  }
+
+  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
+  assert (ctx != NULL);
+
+  if (!ctx->ao_basis.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_ao_basis",
+                           NULL);
+  }
+
+  rc = qmckl_provide_ao_value(context);
+  if (rc != QMCKL_SUCCESS) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_ao_value",
+                           NULL);
+  }
+
+  if (!ctx->mo_basis.provided) {
+    return qmckl_failwith( context,
+                           QMCKL_NOT_PROVIDED,
+                           "qmckl_mo_basis",
+                           NULL);
+  }
+
+ /* Compute if necessary */
+  if (ctx->point.date > ctx->mo_basis.mo_value_date) {
+
+ /* Allocate array */
+    if (ctx->mo_basis.mo_value == NULL) {
+
+      qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
+      mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double);
+      double* mo_value = (double*) qmckl_malloc(context, mem_info);
+
+      if (mo_value == NULL) {
+        return qmckl_failwith( context,
+                               QMCKL_ALLOCATION_FAILED,
+                               "qmckl_mo_basis_mo_value",
+                               NULL);
+      }
+      ctx->mo_basis.mo_value = mo_value;
+    }
+
+    if (ctx->mo_basis.mo_vgl_date == ctx->point.date) {
+
+      // mo_vgl has been computed at this step: Just copy the data.
+      
+      double * v = &(ctx->mo_basis.mo_value[0]);
+      double * vgl = &(ctx->mo_basis.mo_vgl[0]);
+      for (int i=0 ; i<ctx->point.num ; ++i) {
+        for (int k=0 ; k<ctx->mo_basis.mo_num ; ++k) {
+          v[k] = vgl[k];
+        }
+        v   += ctx->mo_basis.mo_num;
+        vgl += ctx->mo_basis.mo_num * 5;
+      }
+
+    } else {
+
+      rc = qmckl_compute_mo_basis_mo_value(context,
+                                           ctx->ao_basis.ao_num,
+                                           ctx->mo_basis.mo_num,
+                                           ctx->point.num,
+                                           ctx->mo_basis.coefficient_t,
+                                           ctx->ao_basis.ao_value,
+                                           ctx->mo_basis.mo_value);
+
+      if (rc != QMCKL_SUCCESS) {
+        return rc;
+      }
+
+    }
+
+    ctx->mo_basis.mo_value_date = ctx->date;
+  }
+
+  return QMCKL_SUCCESS;
+}
+    #+end_src
+
+*** Compute
+   :PROPERTIES:
+   :Name:     qmckl_compute_mo_basis_mo_value
+   :CRetType: qmckl_exit_code
+   :FRetType: qmckl_exit_code
+   :END:
+
+    #+NAME: qmckl_mo_basis_mo_value_args
+    | Variable            | Type                        | In/Out | Description                                     |
+    |---------------------+-----------------------------+--------+-------------------------------------------------|
+    | ~context~           | ~qmckl_context~             | in     | Global state                                    |
+    | ~ao_num~            | ~int64_t~                   | in     | Number of AOs                                   |
+    | ~mo_num~            | ~int64_t~                   | in     | Number of MOs                                   |
+    | ~point_num~         | ~int64_t~                   | in     | Number of points                                |
+    | ~coef_normalized_t~ | ~double[mo_num][ao_num]~    | in     | Transpose of the AO to MO transformation matrix |
+    | ~ao_value~          | ~double[point_num][ao_num]~ | in     | Value of the AOs                                |
+    | ~mo_value~          | ~double[point_num][mo_num]~ | out    | Value of the MOs                                |
+
+
+    The matrix of AO values is very sparse, so we use a sparse-dense
+    matrix multiplication instead of a dgemm, as exposed in
+    https://dx.doi.org/10.1007/978-3-642-38718-0_14.
+
+
+
+    #+begin_src f90 :comments org :tangle (eval f) :noweb yes
+integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
+     ao_num, mo_num, point_num, &
+     coef_normalized_t, ao_value, mo_value) &
+     result(info)
+  use qmckl
+  implicit none
+  integer(qmckl_context), intent(in)  :: context
+  integer*8             , intent(in)  :: ao_num, mo_num
+  integer*8             , intent(in)  :: point_num
+  double precision      , intent(in)  :: ao_value(ao_num,point_num)
+  double precision      , intent(in)  :: coef_normalized_t(mo_num,ao_num)
+  double precision      , intent(out) :: mo_value(mo_num,point_num)
+  integer*8 :: i,j,k
+  double precision :: c1, c2, c3, c4, c5
+
+  integer*8 :: LDA, LDB, LDC
+
+  info = QMCKL_SUCCESS
+  if (.True.)  then    ! fast algorithm
+     do j=1,point_num
+        mo_value(:,j) = 0.d0
+        do k=1,ao_num
+           if (ao_value(k,j) /= 0.d0) then
+              c1 = ao_value(k,j)
+              do i=1,mo_num
+                 mo_value(i,j) = mo_value(i,j) + coef_normalized_t(i,k) * c1
+              end do
+           end if
+        end do
+     end do
+     
+  else ! dgemm
+
+    LDA = size(coef_normalized_t,1)
+    LDB = size(ao_value,1) 
+    LDC = size(mo_value,1)
+
+    info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0,     &
+                                    coef_normalized_t, LDA, ao_value, LDB, &
+                                    0.d0, mo_value, LDC)
+
+  end if
+
+end function qmckl_compute_mo_basis_mo_value_doc_f
+    #+end_src
+
+    #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_value (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const double* coef_normalized_t,
+         const double* ao_value,
+         double* const mo_value );
+   #+end_src
+
+   #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
+
+   #+RESULTS:
+   #+begin_src c :tangle (eval h_func) :comments org
+   qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
+         const qmckl_context context,
+         const int64_t ao_num,
+         const int64_t mo_num,
+         const int64_t point_num,
+         const double* coef_normalized_t,
+         const double* ao_value,
+         double* const mo_value );
+   #+end_src
+
+   #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc"))
+
+    #+RESULTS:
+    #+begin_src f90 :tangle (eval f) :comments org :exports none
+    integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc &
+        (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value) &
+        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 :: ao_num
+      integer (c_int64_t) , intent(in)  , value :: mo_num
+      integer (c_int64_t) , intent(in)  , value :: point_num
+      real    (c_double ) , intent(in)          :: coef_normalized_t(ao_num,mo_num)
+      real    (c_double ) , intent(in)          :: ao_value(ao_num,point_num)
+      real    (c_double ) , intent(out)         :: mo_value(mo_num,point_num)
+
+      integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f
+      info = qmckl_compute_mo_basis_mo_value_doc_f &
+             (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value)
+
+    end function qmckl_compute_mo_basis_mo_value_doc
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value (const qmckl_context context,
+                                 const int64_t ao_num,
+                                 const int64_t mo_num,
+                                 const int64_t point_num,
+                                 const double* coef_normalized_t,
+                                 const double* ao_value,
+                                 double* const mo_value )
+{
+#ifdef HAVE_HPC
+  return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value);
+#else
+  return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value);
+#endif
+}
+    #+end_src
+
+*** HPC version
+
+
+    #+begin_src c :tangle (eval h_func) :comments org
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
+                                     const int64_t ao_num,
+                                     const int64_t mo_num,
+                                     const int64_t point_num,
+                                     const double* coef_normalized_t,
+                                     const double* ao_value,
+                                     double* const mo_value );
+#endif
+    #+end_src
+
+    #+begin_src c :tangle (eval c) :comments org
+#ifdef HAVE_HPC
+qmckl_exit_code
+qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
+                                     const int64_t ao_num,
+                                     const int64_t mo_num,
+                                     const int64_t point_num,
+                                     const double* restrict coef_normalized_t,
+                                     const double* restrict ao_value,
+                                     double* restrict const mo_value )
+{
+  assert (context != QMCKL_NULL_CONTEXT);
+
+#ifdef HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
+    double* restrict const vgl1  = &(mo_value[ipoint*mo_num]);
+    const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);
+
+    for (int64_t i=0 ; i<mo_num ; ++i) {
+      vgl1[i] = 0.;
+    }
+
+    int64_t nidx=0;
+    int64_t idx[ao_num];
+    double  av1[ao_num];
+    for (int64_t k=0 ; k<ao_num ; ++k) {
+      if (avgl1[k] != 0.) {
+        idx[nidx] = k;
+        av1[nidx] = avgl1[k];
+        ++nidx;
+      }
+    }
+
+    int64_t n;
+
+    for (n=0 ; n < nidx-4 ; n+=4) {
+      const double* restrict ck1 = coef_normalized_t + idx[n  ]*mo_num;
+      const double* restrict ck2 = coef_normalized_t + idx[n+1]*mo_num;
+      const double* restrict ck3 = coef_normalized_t + idx[n+2]*mo_num;
+      const double* restrict ck4 = coef_normalized_t + idx[n+3]*mo_num;
+
+      const double a11 = av1[n  ];
+      const double a21 = av1[n+1];
+      const double a31 = av1[n+2];
+      const double a41 = av1[n+3];
+
+#ifdef HAVE_OPENMP
+#pragma omp simd
+#endif
+      for (int64_t i=0 ; i<mo_num ; ++i) {
+        vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
+      }
+    }
+
+    const int64_t n0 = n < 0 ? 0 : n;
+    for (int64_t m=n0 ; m < nidx ; m+=1) {
+      const double* restrict ck = coef_normalized_t + idx[m]*mo_num;
+      const double a1 = av1[m];
+
+#ifdef HAVE_OPENMP
+  #pragma omp simd
+#endif
+      for (int64_t i=0 ; i<mo_num ; ++i) {
+        vgl1[i] += ck[i] * a1;
+      }
+    }
+  }
+  return QMCKL_SUCCESS;
+}
+#endif
+    #+end_src
+
+** Computation of MOs: values, gradient, Laplacian
 
 *** Get
 
@@ -774,7 +1235,6 @@ qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
 }
     #+end_src
 
-
 *** HPC version
 
 
@@ -918,9 +1378,9 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
 #endif
     #+end_src
 
-*** Test
+** Test
 
-    #+begin_src python :results output :exports none
+   #+begin_src python :results output :exports none
 import numpy as np
 
 def f(a,x,y):
@@ -980,9 +1440,9 @@ print ( "[2][1][15][14] : %25.15e"% df(a,x,y,2))
 print ( "[3][1][15][14] : %25.15e"% df(a,x,y,3))
 print ( "[4][1][15][14] : %25.15e"% lf(a,x,y))
 
-    #+end_src
+   #+end_src
 
-     #+begin_src c :tangle (eval c_test) :exports none
+    #+begin_src c :tangle (eval c_test) :exports none
 {
 #define walk_num chbrclf_walk_num
 #define elec_num chbrclf_elec_num
@@ -1107,10 +1567,36 @@ assert (rc == QMCKL_SUCCESS);
 
 assert(qmckl_mo_basis_provided(context));
 
+rc = qmckl_context_touch(context);
+assert (rc == QMCKL_SUCCESS);
+
+double mo_value[walk_num*elec_num][chbrclf_mo_num];
+rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
+assert (rc == QMCKL_SUCCESS);
+
 double mo_vgl[walk_num*elec_num][5][chbrclf_mo_num];
 rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), walk_num*elec_num*5*chbrclf_mo_num);
 assert (rc == QMCKL_SUCCESS);
 
+for (int i=0 ; i< walk_num*elec_num; ++i) {
+  for (int k=0 ; k< chbrclf_mo_num ; ++k) {
+    assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
+  }
+ }
+
+rc = qmckl_context_touch(context);
+assert (rc == QMCKL_SUCCESS);
+ 
+rc = qmckl_get_mo_basis_mo_value(context, &(mo_value[0][0]), walk_num*elec_num*chbrclf_mo_num);
+assert (rc == QMCKL_SUCCESS);
+
+for (int i=0 ; i< walk_num*elec_num; ++i) {
+  for (int k=0 ; k< chbrclf_mo_num ; ++k) {
+    assert(fabs(mo_vgl[i][0][k] - mo_value[i][k]) < 1.e-12) ;
+  }
+ }
+
+
 // Test overlap of MO
 //double point_x[10];
 //double point_y[10];
@@ -1168,7 +1654,7 @@ printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[2][1][3]);
 printf("\n");
 }
 
-     #+end_src
+    #+end_src
 
 * End of files                                                     :noexport: