From 5060bde30fd224285b86a7f33d25edd1181b33f7 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 30 Jan 2024 23:46:06 +0100
Subject: [PATCH 1/2] Moved ivdep after omp simd

---
 org/qmckl_ao.org | 16 ++++++++--------
 org/qmckl_mo.org | 32 ++++++++++++++++----------------
 2 files changed, 24 insertions(+), 24 deletions(-)

diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index db8a9c6..03e46f4 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -5944,28 +5944,28 @@ IVDEP
               ao_value_1[0] = s1 * f[0];
               break;
             case 3:
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
               for (int il=0 ; il<3 ; ++il) {
                 ao_value_1[il] =  poly_vgl_1[il] * s1 * f[il];
               }
               break;
             case(6):
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
               for (int il=0 ; il<6 ; ++il) {
                 ao_value_1[il] =  poly_vgl_1[il] * s1 * f[il];
               }
               break;
             default:
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
               for (int il=0 ; il<n ; ++il) {
                 ao_value_1[il] =  poly_vgl_1[il] * s1 * f[il];
               }
@@ -6713,10 +6713,10 @@ qmckl_compute_ao_vgl_hpc_gaussian (
 
         //        if (do_sparse) {
         for (int i=0 ; i<nucleus_shell_num[inucl] ; ++i) {
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
           for (int j=0 ; j<8 ; ++j) {
             ce_mat[i][j] = 0.;
           }
@@ -6725,10 +6725,10 @@ IVDEP
           for (int l=1 ; l<idx[0]; ++l) {
             const int k = idx[l];
             if (k >= nidx) break;
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
             for (int j=0 ; j<8 ; ++j) {
               ce_mat[i][j] = ce_mat[i][j] + v[l] * exp_mat[k][j];
             }
@@ -6814,10 +6814,10 @@ IVDEP
               ao_vgl_5[0] = s5 * f[0];
               break;
             case 3:
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
               for (int il=0 ; il<3 ; ++il) {
                 ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
                 ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
@@ -6830,10 +6830,10 @@ IVDEP
               }
               break;
             case 6:
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
               for (int il=0 ; il<6 ; ++il) {
                 ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
                 ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
@@ -6846,10 +6846,10 @@ IVDEP
               }
               break;
             default:
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
               for (int il=0 ; il<n ; ++il) {
                 ao_vgl_1[il] =  poly_vgl_1[il] * s1 * f[il];
                 ao_vgl_2[il] = (poly_vgl_2[il] * s1 + poly_vgl_1[il] * s2) * f[il];
diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org
index b5028d3..efb2179 100644
--- a/org/qmckl_mo.org
+++ b/org/qmckl_mo.org
@@ -1301,10 +1301,10 @@ qmckl_compute_mo_basis_mo_value_hpc_sp (const qmckl_context context,
                            "coefficient_t_sp");
   };
 
-  IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+  IVDEP
     for (int64_t i=0 ; i<mo_num*ao_num ; ++i) {
       coefficient_t_sp[i] = (float) coefficient_t[i];
     }
@@ -1358,10 +1358,10 @@ qmckl_compute_mo_basis_mo_value_hpc_sp (const qmckl_context context,
         const float a71 = av1[n+6];
         const float a81 = av1[n+7];
 
-        IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+        IVDEP
           for (int64_t i=0 ; i<mo_num ; ++i) {
             vgl_sp[i] = vgl_sp[i] +
                          ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41 +
@@ -1373,19 +1373,19 @@ qmckl_compute_mo_basis_mo_value_hpc_sp (const qmckl_context context,
         const float* restrict ck = coefficient_t_sp + idx[m]*mo_num;
         const float a1 = av1[m];
 
-        IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+        IVDEP
           for (int64_t i=0 ; i<mo_num ; ++i) {
             vgl_sp[i] = vgl_sp[i] + ck[i] * a1;
           }
       }
 
-      IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+      IVDEP
       for (int64_t i=0 ; i<mo_num ; ++i) {
         vgl1[i] = (double) vgl_sp[i];
       }
@@ -1482,10 +1482,10 @@ qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
         const double a71 = av1[n+6];
         const double a81 = av1[n+7];
 
-        IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+        IVDEP
           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 +
@@ -1497,10 +1497,10 @@ qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
         const double* restrict ck = coefficient_t + idx[m]*mo_num;
         const double a1 = av1[m];
 
-        IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+        IVDEP
           for (int64_t i=0 ; i<mo_num ; ++i) {
             vgl1[i] = vgl1[i] + ck[i] * a1;
           }
@@ -2038,10 +2038,10 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
       const double a35 = av5[n+2];
       const double a45 = av5[n+3];
 
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
       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;
         vgl2[i] = vgl2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
@@ -2059,10 +2059,10 @@ IVDEP
       const double a4 = av4[m];
       const double a5 = av5[m];
 
-IVDEP
 #ifdef HAVE_OPENMP
   #pragma omp simd
 #endif
+IVDEP
       for (int64_t i=0 ; i<mo_num ; ++i) {
         vgl1[i] = vgl1[i] + ck[i] * a1;
         vgl2[i] = vgl2[i] + ck[i] * a2;
@@ -2120,10 +2120,10 @@ qmckl_compute_mo_basis_mo_vgl_hpc_sp (const qmckl_context context,
                            "coefficient_t_sp");
   };
 
-  IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+  IVDEP
     for (int64_t i=0 ; i<mo_num*ao_num ; ++i) {
       coefficient_t_sp[i] = (float) coefficient_t[i];
     }
@@ -2219,10 +2219,10 @@ qmckl_compute_mo_basis_mo_vgl_hpc_sp (const qmckl_context context,
         const float a35 = av5[n+2];
         const float a45 = av5[n+3];
 
-        IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+        IVDEP
           for (int64_t i=0 ; i<mo_num ; ++i) {
             vgl_sp1[i] = vgl_sp1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
             vgl_sp2[i] = vgl_sp2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
@@ -2240,10 +2240,10 @@ qmckl_compute_mo_basis_mo_vgl_hpc_sp (const qmckl_context context,
         const float a4 = av4[m];
         const float a5 = av5[m];
 
-        IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+        IVDEP
           for (int64_t i=0 ; i<mo_num ; ++i) {
             vgl_sp1[i] = vgl_sp1[i] + ck[i] * a1;
             vgl_sp2[i] = vgl_sp2[i] + ck[i] * a2;
@@ -2252,10 +2252,10 @@ qmckl_compute_mo_basis_mo_vgl_hpc_sp (const qmckl_context context,
             vgl_sp5[i] = vgl_sp5[i] + ck[i] * a5;
           }
       }
-      IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+      IVDEP
       for (int64_t i=0 ; i<mo_num ; ++i) {
         vgl1[i] = (double) vgl_sp1[i];
         vgl2[i] = (double) vgl_sp2[i];
@@ -2583,10 +2583,10 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context,
       const double a31 = av1[n+2];
       const double a41 = av1[n+3];
 
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
       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;
       }
@@ -2596,10 +2596,10 @@ IVDEP
       const double* restrict ck = coefficient_t + idx[m]*mo_num;
       const double a1 = av1[m];
 
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
       for (int64_t i=0 ; i<mo_num ; ++i) {
         vgl1[i] = vgl1[i] + ck[i] * a1;
       }
@@ -3008,10 +3008,10 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context,
       const double a35 = av5[n+2];
       const double a45 = av5[n+3];
 
-IVDEP
 #ifdef HAVE_OPENMP
 #pragma omp simd
 #endif
+IVDEP
       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;
         vgl2[i] = vgl2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
@@ -3029,10 +3029,10 @@ IVDEP
       const double a4 = av4[m];
       const double a5 = av5[m];
 
-IVDEP
 #ifdef HAVE_OPENMP
   #pragma omp simd
 #endif
+IVDEP
       for (int64_t i=0 ; i<mo_num ; ++i) {
         vgl1[i] = vgl1[i] + ck[i] * a1;
         vgl2[i] = vgl2[i] + ck[i] * a2;

From a0e1843963746ffa7c9a39ad502146dacacfe49f Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 6 Feb 2024 22:27:17 +0100
Subject: [PATCH 2/2] Fix memset

---
 org/qmckl_ao.org | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org
index db8a9c6..f09faff 100644
--- a/org/qmckl_ao.org
+++ b/org/qmckl_ao.org
@@ -5759,8 +5759,6 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
 {
   <<ao_value_constants>>
 
-    memset(ao_value, 0, ao_num*point_num*sizeof(double));
-
 #ifdef HAVE_OPENMP
 #pragma omp parallel if (point_num > 16)
 #endif
@@ -5800,6 +5798,9 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context,
 #pragma omp for
 #endif
     for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
+
+      memset(&ao_value[ipoint*ao_num], 0, ao_num*sizeof(double));
+
       const double e_coord[3]  __attribute__((aligned(64))) =
         { coord[ipoint],
           coord[ipoint + point_num],
@@ -6537,8 +6538,6 @@ qmckl_compute_ao_vgl_hpc_gaussian (
 {
   <<ao_value_constants>>
 
-    memset(ao_vgl, 0, 5*ao_num*point_num*sizeof(double));
-
 #ifdef HAVE_OPENMP
 #pragma omp parallel if (point_num > 16)
 #endif
@@ -6590,6 +6589,9 @@ qmckl_compute_ao_vgl_hpc_gaussian (
 #pragma omp for
 #endif
     for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
+
+      memset(&ao_vgl[ipoint*ao_num*5], 0, 5*ao_num*sizeof(double));
+
       const double e_coord[3]  __attribute__((aligned(64))) =
         { coord[ipoint],
           coord[ipoint + point_num],